code: package cc.mallet.types;
  
  import java.io.*;
  
  import cc.mallet.util.Maths;
  
  public abstract class DenseMatrix implements Matrix, Serializable
  {
  	double[] values;
  	protected boolean hasInfinite;             // if true, at least one value =  -Inf or +Inf
  	
  	public abstract int getNumDimensions ();
  	public abstract int getDimensions (int[] sizes);
  	
  	public abstract double value (int[] indices);
  	public abstract void setValue (int[] indices, double value);
  	public abstract ConstantMatrix cloneMatrix ();
  
  	public abstract int singleIndex (int[] indices);
  	public abstract void singleToIndices (int i, int[] indices);
  	public double singleValue (int i) { return values[i]; }
  	public void setSingleValue (int i, double value) { values[i] = value; }
  	public void incrementSingleValue (int i, double delta) { values[i] += delta; }
  
    public void setValueAtLocation (int loc, double value)
    {
      // indices == locations
      setSingleValue (loc, value);
    }
  
    public int singleSize () { return values.length; }
  
  	public int numLocations () { return values.length; }
  	public int location (int index) { return index; }
  	public double valueAtLocation (int location) { return values[location]; }
  	// Returns a "singleIndex"
  	public int indexAtLocation (int location) { return location; }
  
  	
  	public void setAll (double v) { for (int i = 0; i < values.length; i++) values[i] = v; }
  	public void set (ConstantMatrix m) {
  		if (m instanceof DenseMatrix) {
  			assert (m.singleSize() == values.length);
  			System.arraycopy (((DenseMatrix)m).values, 0, values, 0, values.length);
  		} else
  			for (int i = m.numLocations()-1; i >= 0; i--)
  				values[m.indexAtLocation(i)] = m.valueAtLocation(i);
  	}
  	
  	public void setWithAddend (ConstantMatrix m, double addend) {
  		if (m instanceof DenseMatrix) {
  			assert (m.singleSize() == values.length);
  			for (int i = 0; i < values.length; i++)
  				values[i] = ((DenseMatrix)m).values[i] + addend;
  		} else 
  			for (int i = m.numLocations()-1; i >= 0; i--)
  				values[m.indexAtLocation(i)] = m.valueAtLocation(i) + addend;
  	}
  		
  	public void setWithFactor (ConstantMatrix m, double factor) {
  		if (m instanceof DenseMatrix) {
  			assert (m.singleSize() == values.length);
  			for (int i = 0; i < values.length; i++)
  				values[i] = ((DenseMatrix)m).values[i] * factor;
  		} else 
  			for (int i = m.numLocations()-1; i >= 0; i--)
  				values[m.indexAtLocation(i)] = m.valueAtLocation(i) * factor;
  	}
  
  	public void plusEquals (double v) {
  		for (int i = 0; i < values.length; i++)
  			values[i] += v;
  	}
  	
  	public void plusEquals (ConstantMatrix m) {
  		if (m instanceof DenseMatrix) {
  			assert (m.singleSize() == values.length);
  			for (int i = 0; i < values.length; i++) {
  				// added by Culotta - 12.10.0 to enforce INF - INF = 0
  				if(Double.isInfinite(values[i]) &&
  					 Double.isInfinite(m.valueAtLocation(i))) {
  					double newValue = m.valueAtLocation(i);
  					// make sure they're opposite signed 
  					if((newValue * values[i]) < 0) {
  						values[i] = 0.0; // inf - inf = 0
  					}
  					else
  						values[i] += newValue;
  				}
  				else 
  					values[i] += m.valueAtLocation(i);
  			}
  		}
  		else  
  			for (int i = m.numLocations()-1; i >= 0; i--) {
  				// added by Culotta - 12.10.02 to enforce INF - INF = 0
  				if(Double.isInfinite(values[m.indexAtLocation(i)]) &&
  					 Double.isInfinite(((DenseMatrix)m).values[i])) {
  					double newValue = m.valueAtLocation(i);
  					// make sure they're oppisite signed 
  					if((newValue * values[m.indexAtLocation(i)]) < 0) {
  						values[m.indexAtLocation(i)] = 0.0;
  					}
  					else
  						values[m.indexAtLocation(i)] += newValue;
  				}
  				else
  					values[m.indexAtLocation(i)] += m.valueAtLocation(i);
  			}
  	}
  
  	public void plusEquals (ConstantMatrix m, double factor) {
  		if (m instanceof DenseMatrix) {
  			assert (m.singleSize() == values.length);
  			for (int i = 0; i < values.length; i++) {
  				// added by Culotta - 12.10.0 to enforce INF - INF = 0
  				if(Double.isInfinite(values[i]) &&
  					 Double.isInfinite(m.valueAtLocation(i))) {
  					double newValue = factor*(m.valueAtLocation(i));
  					// make sure they're opposite signed 
  					if((newValue * values[i]) < 0) {
  						values[i] = 0.0; // inf - inf = 0
  					}
  					else
  						values[i] += newValue;
  				}
  				else
  					values[i] += (m.valueAtLocation(i) * factor);
  			}
  		}
  		else 
  			for (int i = m.numLocations()-1; i >= 0; i--){
  				// added by Culotta - 12.10.02 to enforce INF - INF = 0
  				if(Double.isInfinite(values[m.indexAtLocation(i)]) &&
  					 Double.isInfinite(m.valueAtLocation(i))) {
  					double newValue = factor*m.valueAtLocation(i);
  					// make sure they're oppisite signed 
  					if((newValue * values[m.indexAtLocation(i)]) < 0) {
  						values[m.indexAtLocation(i)] = 0.0;
  					}
  					else
  						values[m.indexAtLocation(i)] += newValue;
  				}
  				else
  					values[m.indexAtLocation(i)] += m.valueAtLocation(i) * factor;
  			}
  	}
  
  
  
  	public void equalsPlus (double factor, ConstantMatrix m) {
  		if (m instanceof DenseMatrix) {
  			assert (m.singleSize() == values.length);
  			for (int i = 0; i < values.length; i++) {
  				// added by Culotta - 12.10.0 to enforce INF - INF = 0
  				if(Double.isInfinite(values[i]) &&
  					 Double.isInfinite(((DenseMatrix)m).values[i])) {
  					double lhs = factor*values[i];
  					double rhs = ((DenseMatrix)m).values[i];
  					// make sure they're opposite signed 
  					if((lhs * rhs) < 0) {
  						values[i] = 0.0; // inf - inf = 0
  					}
  					else
  						values[i] = lhs + rhs;
  				}
  				else
  					values[i] = factor*values[i] + ((DenseMatrix)m).values[i];
  			}
  		}
  		else 
  			for (int i = m.numLocations()-1; i >= 0; i--) {				
  				// added by Culotta - 12.10.02 to enforce INF - INF = 0
  				if(Double.isInfinite(values[m.indexAtLocation(i)]) &&
  					 Double.isInfinite(((DenseMatrix)m).values[i])) {
  					double lhs = factor * values[m.indexAtLocation(i)];
  					double rhs = m.valueAtLocation(i);
  					// make sure they're oppisite signed 
  					if((lhs * rhs) < 0) {
  						values[m.indexAtLocation(i)] = 0.0;
  					}
  					else
  						values[m.indexAtLocation(i)] = lhs + rhs;
  				}
  				else
  					values[m.indexAtLocation(i)] = factor * values[m.indexAtLocation(i)] + m.valueAtLocation(i);
  			}
  	}
  	
  	public void timesEquals (double factor) {
  		for (int i = 0; i < values.length; i++)
  			values[i] *= factor;
  	}
  
  	public void elementwiseTimesEquals (ConstantMatrix m) {
  		if (m instanceof DenseMatrix) {
  			assert (m.singleSize() == values.length);
  			for (int i = 0; i < values.length; i++)
  				values[i] *= ((DenseMatrix)m).values[i];
  		} else 
  			for (int i = m.numLocations()-1; i >= 0; i--)
  				values[m.indexAtLocation(i)] *= m.valueAtLocation(i);
  	}
  	
  	public void elementwiseTimesEquals (ConstantMatrix m, double factor) {
  		if (m instanceof DenseMatrix) {
  			assert (m.singleSize() == values.length);
  			for (int i = 0; i < values.length; i++)
  				values[i] *= ((DenseMatrix)m).values[i] * factor;
  		} else 
  			for (int i = m.numLocations()-1; i >= 0; i--)
  				values[m.indexAtLocation(i)] *= m.valueAtLocation(i) * factor;
  	}
  	
  	public void divideEquals (double factor) {
  		for (int i = 0; i < values.length; i++)
  			values[i] /= factor;
  	}
  
  	public void elementwiseDivideEquals (ConstantMatrix m) {
  		if (m instanceof DenseMatrix) {
  			assert (m.singleSize() == values.length);
  			for (int i = 0; i < values.length; i++)
  				values[i] /= ((DenseMatrix)m).values[i];
  		} else 
  			for (int i = m.numLocations()-1; i >= 0; i--)
  				values[m.indexAtLocation(i)] /= m.valueAtLocation(i);
  	}
  
  	public void elementwiseDivideEquals (ConstantMatrix m, double factor) {
  		if (m instanceof DenseMatrix) {
  			assert (m.singleSize() == values.length);
  			for (int i = 0; i < values.length; i++)
  				values[i] /= ((DenseMatrix)m).values[i] * factor;
  		} else 
  			for (int i = m.numLocations()-1; i >= 0; i--)
  				values[m.indexAtLocation(i)] /= m.valueAtLocation(i) * factor;
  	}
  
  	// xxx Perhaps make a special efficient case for binary vectors
  	public double dotProduct (ConstantMatrix m) {
  		double ret = 0;
  		if (m instanceof DenseMatrix) {
  			assert (m.singleSize() == values.length);
  			for (int i = 0; i < values.length; i++)
  				ret += values[i] * ((DenseMatrix)m).values[i];
  		} else {
  			for (int i = m.numLocations()-1; i >= 0; i--)
  				if(m.indexAtLocation(i) < values.length)//fix problem
  					ret += values[m.indexAtLocation(i)] * m.valueAtLocation(i);
  				else{
  //					System.out.println(m.indexAtLocation(i) + ":" + values.length);
  //					throw new ArrayIndexOutOfBoundsException(m.indexAtLocation(i));
  				}
  		}
  
  		return ret;
  	}
  	
  	public double absNorm() {
  		double ret = 0;
  		for (int i = 0; i < values.length; i++)
  			ret += Math.abs(values[i]);
  		return ret;
  	}
  
  	public double oneNorm () {
  		double ret = 0;
  		for (int i = 0; i < values.length; i++)
  			ret += values[i];
  		return ret;
  	}
  	
  	public double twoNorm () {
  		double ret = 0;
  		for (int i = 0; i < values.length; i++)
  			ret += values[i] * values[i];
  		return Math.sqrt (ret);
  	}
  	
  	public double infinityNorm () {
  		double max = Double.NEGATIVE_INFINITY;
  		for (int i = 0; i < values.length; i++)
  			if (Math.abs(values[i]) > max)
  				max = Math.abs(values[i]);
  		return max;
  	}		
  
  	public double oneNormalize ()	{
  		double norm = oneNorm();
  		for (int i = 0; i < values.length; i++)
  			values[i] /= norm;
  		return norm;
  	}
  
  	public double twoNormalize ()	{
  		double norm = twoNorm();
  		for (int i = 0; i < values.length; i++)
  			values[i] /= norm;
  		return norm;
  	}
  
  	public double absNormalize ()	{
  		double norm = absNorm();
  		if (norm > 0)
  			for (int i = 0; i < values.length; i++)
  				values[i] /= norm;
  		return norm;
  	}
  
  	public double infinityNormalize () {
  		double norm = infinityNorm();
  		for (int i = 0; i < values.length; i++)
  			values[i] /= norm;
  		return norm;
  	}
  	
  	public void print() {
  		for (int i = 0; i < values.length; i++)
  			System.out.println ("DenseMatrix["+i+"] = "+values[i]);
  	}
  		
  	public boolean isNaN() {
  		for (int i = 0; i < values.length; i++)
  			if (Double.isNaN(values[i]))
  				return true;
  		return false;
  	}
  
  
  	public final void substitute (double oldValue, double newValue)
  	{
  		for (int i = values.length-1; i >= 0; i--)
  			if (values[i] == oldValue)
  				values[i] = newValue;
  	}
  	
  	// Serialization
  		
  	private static final long serialVersionUID = 1;
  	private static final int CURRENT_SERIAL_VERSION = 0;
  	private static final int NULL_INTEGER = -1;
  	
  	private void writeObject (ObjectOutputStream out) throws IOException {
  		int i, size;
  		out.writeInt (CURRENT_SERIAL_VERSION);
  		if (values != null) {
  			size = values.length;
  			out.writeInt(size);
  			for (i=0; i<size; i++) {
  				out.writeDouble(values[i]);
  			}
  		}
  		else {
  			out.writeInt(NULL_INTEGER);
  		}
  	}
  	
  	private void readObject (ObjectInputStream in) throws IOException, ClassNotFoundException {
  		int i, size;
  		this.hasInfinite = false;
  		int version = in.readInt ();
  		size = in.readInt();
  		if (size != NULL_INTEGER) {
  			values = new double[size];
  			for (i = 0; i<size; i++) {
  				values[i] = in.readDouble();
  				if (Double.isInfinite (values[i]))
  					this.hasInfinite = true;
  			}
  		}
  		else {
  			values = null;
  		}
  	}
  
  
  	public static void plusEquals (double[] accumulator, double[] addend)
  	{
  		assert (accumulator.length == addend.length);
  		for (int i = 0; i < addend.length; i++)
  			accumulator[i] += addend[i];
  	}
  
  	public static void plusEquals (double[] accumulator, double[] addend, double factor)
  	{
  		assert (accumulator.length == addend.length);
  		for (int i = 0; i < addend.length; i++)
  			accumulator[i] += factor * addend[i];
  	}
  	
  	public static void timesEquals (double[] accumulator, double[] product)
  	{
  		assert (accumulator.length == product.length);
  		for (int i = 0; i < product.length; i++)
  			accumulator[i] *= product[i];
  	}
  
  	public static double infinityNorm (double[] vector)
  	{
  		double max = Double.NEGATIVE_INFINITY;
  		for (int i = 0; i < vector.length; i++)
  			if (Math.abs(vector[i]) > max)
  				max = Math.abs(vector[i]);
  		return max;
  	}
  	
  	// This should probably be generalized.
  	public boolean almostEquals (ConstantMatrix m2) {
  		if (getNumDimensions () != m2.getNumDimensions ()) {
  			return false;
  		}
  		if (numLocations () != m2.numLocations ()) {
  			return false;
  		}
  		int[] dims1 = new int [getNumDimensions ()];
  		int[] dims2 = new int [getNumDimensions ()];
  		getDimensions (dims1);
  		m2.getDimensions (dims2);
  		for (int i = 0; i < dims1.length; i++) {
  			if (dims1 [i] != dims2 [i]) {
  				return false;
  			}
  		}
  		for (int i = 0; i < numLocations(); i++) {
  			if (!Maths.almostEquals (valueAtLocation (i), m2.valueAtLocation (i))) {
  				return false;
  			}
  		}
  		return true;
  	}
  		
  }
  