Java Code Examples for cern.colt.matrix.DoubleMatrix1D#size()

public void dsymv(boolean isUpperTriangular, double alpha, DoubleMatrix2D A, DoubleMatrix1D x, double beta, DoubleMatrix1D y) {
	if (isUpperTriangular) A = A.viewDice();
	Property.DEFAULT.checkSquare(A);
	int size = A.rows();
	if (size != x.size() || size!=y.size()) {
		throw new IllegalArgumentException(A.toStringShort() + ", " + x.toStringShort() + ", " + y.toStringShort());
	}
	DoubleMatrix1D tmp = x.like();
	for (int i = 0; i < size; i++) {
		double sum = 0;
		for (int j = 0; j <= i; j++) {
			sum += A.getQuick(i,j) * x.getQuick(j);
		}
		for (int j = i + 1; j < size; j++) {
			sum += A.getQuick(j,i) * x.getQuick(j);
		}
		tmp.setQuick(i, alpha * sum + beta * y.getQuick(i));
	}
	y.assign(tmp);
}
 

/**
 * Returns whether both given matrices <tt>A</tt> and <tt>B</tt> are equal.
 * The result is <tt>true</tt> if <tt>A==B</tt>. 
 * Otherwise, the result is <tt>true</tt> if and only if both arguments are <tt>!= null</tt>, 
 * have the same size and 
 * <tt>! (Math.abs(A[i] - B[i]) > tolerance())</tt> holds for all indexes.
 * @param   A   the first matrix to compare.
 * @param   B   the second matrix to compare.
 * @return  <tt>true</tt> if both matrices are equal;
 *          <tt>false</tt> otherwise.
 */
public boolean equals(DoubleMatrix1D A, DoubleMatrix1D B) {
	if (A==B) return true;
	if (! (A != null && B != null)) return false;
	int size = A.size();
	if (size != B.size()) return false;

	double epsilon = tolerance();
	for (int i=size; --i >= 0;) {
		//if (!(getQuick(i) == B.getQuick(i))) return false;
		//if (Math.abs(A.getQuick(i) - B.getQuick(i)) > epsilon) return false;
		double x = A.getQuick(i);
		double value = B.getQuick(i);
		double diff = Math.abs(value - x);
		if ((diff!=diff) && ((value!=value && x!=x) || value==x)) diff = 0;
		if (!(diff <= epsilon)) return false;
	}
	return true;
}
 

/**
 * Linear algebraic matrix-vector multiplication; <tt>z = A * y</tt>.
 * <tt>z[i] = alpha*Sum(A[i,j] * y[j]) + beta*z[i], i=0..A.rows()-1, j=0..y.size()-1</tt>.
 * Where <tt>A == this</tt>.
 * @param y the source vector.
 * @param z the vector where results are to be stored.
 * 
 * @throws IllegalArgumentException if <tt>A.columns() != y.size() || A.rows() > z.size())</tt>.
 */
protected void zMult(final DoubleMatrix1D y, final DoubleMatrix1D z, cern.colt.list.IntArrayList nonZeroIndexes, DoubleMatrix1D[] allRows, final double alpha, final double beta) {
	if (columns != y.size() || rows > z.size())	
		throw new IllegalArgumentException("Incompatible args: "+toStringShort()+", "+y.toStringShort()+", "+z.toStringShort());

	z.assign(cern.jet.math.Functions.mult(beta/alpha));
	for (int i = indexes.length; --i >= 0; ) {
		if (indexes[i] != null) {
			for (int k = indexes[i].size(); --k >= 0; ) {
				int j = indexes[i].getQuick(k);
				double value = values[i].getQuick(k);
				z.setQuick(i,z.getQuick(i) + value * y.getQuick(j));
			}
		}
	}
		
	z.assign(cern.jet.math.Functions.mult(alpha));
}
 

/**
 * Linear algebraic matrix-vector multiplication; <tt>z = A * y</tt>.
 * <tt>z[i] = alpha*Sum(A[i,j] * y[j]) + beta*z[i], i=0..A.rows()-1, j=0..y.size()-1</tt>.
 * Where <tt>A == this</tt>.
 * @param y the source vector.
 * @param z the vector where results are to be stored.
 * 
 * @throws IllegalArgumentException if <tt>A.columns() != y.size() || A.rows() > z.size())</tt>.
 */
protected void zMult(final DoubleMatrix1D y, final DoubleMatrix1D z, cern.colt.list.IntArrayList nonZeroIndexes, DoubleMatrix1D[] allRows, final double alpha, final double beta) {
	if (columns != y.size() || rows > z.size())	
		throw new IllegalArgumentException("Incompatible args: "+toStringShort()+", "+y.toStringShort()+", "+z.toStringShort());

	z.assign(cern.jet.math.Functions.mult(beta/alpha));
	for (int i = indexes.length; --i >= 0; ) {
		if (indexes[i] != null) {
			for (int k = indexes[i].size(); --k >= 0; ) {
				int j = indexes[i].getQuick(k);
				double value = values[i].getQuick(k);
				z.setQuick(i,z.getQuick(i) + value * y.getQuick(j));
			}
		}
	}
		
	z.assign(cern.jet.math.Functions.mult(alpha));
}
 

/**
Constructs and returns a sampling view with a size of <tt>round(matrix.size() * fraction)</tt>.
Samples "without replacement" from the uniform distribution.
@param matrix any matrix.
@param rowFraction the percentage of rows to be included in the view.
@param columnFraction the percentage of columns to be included in the view.
@param randomGenerator a uniform random number generator; set this parameter to <tt>null</tt> to use a default generator seeded with the current time.
@return the sampling view.
@throws IllegalArgumentException if <tt>! (0 <= rowFraction <= 1 && 0 <= columnFraction <= 1)</tt>.
@see cern.jet.random.sampling.RandomSampler
*/
public static DoubleMatrix1D viewSample(DoubleMatrix1D matrix, double fraction, RandomEngine randomGenerator) {
	// check preconditions and allow for a little tolerance
	double epsilon = 1e-09;
	if (fraction < 0 - epsilon || fraction > 1 + epsilon) throw new IllegalArgumentException();
	if (fraction < 0) fraction = 0;
	if (fraction > 1) fraction = 1;

	// random generator seeded with current time
	if (randomGenerator==null) randomGenerator = new cern.jet.random.engine.MersenneTwister((int) System.currentTimeMillis());

	int ncols = (int) Math.round(matrix.size() * fraction);
	int max = ncols;
	long[] selected = new long[max]; // sampler works on long's, not int's

	// sample 
	int n = ncols;
	int N = matrix.size();
	cern.jet.random.sampling.RandomSampler.sample(n,N,n,0,selected,0,randomGenerator);
	int[] selectedCols = new int[n];
	for (int i=0; i<n; i++) selectedCols[i] = (int) selected[i];

	return matrix.viewSelection(selectedCols);
}
 

public void dtrmv(boolean isUpperTriangular, boolean transposeA, boolean isUnitTriangular, DoubleMatrix2D A, DoubleMatrix1D x) {
	if (transposeA) {
		A = A.viewDice();
		isUpperTriangular = !isUpperTriangular;
	}
	
	Property.DEFAULT.checkSquare(A);
	int size = A.rows();
	if (size != x.size()) {
		throw new IllegalArgumentException(A.toStringShort() + ", " + x.toStringShort());
	}
	    
	DoubleMatrix1D b = x.like();
	DoubleMatrix1D y = x.like();
	if (isUnitTriangular) {
		y.assign(1);
	}
	else {
		for (int i = 0; i < size; i++) {
			y.setQuick(i, A.getQuick(i,i));
		}
	}
	
	for (int i = 0; i < size; i++) {
		double sum = 0;
		if (!isUpperTriangular) {
			for (int j = 0; j < i; j++) {
				sum += A.getQuick(i,j) * x.getQuick(j);
			}
			sum += y.getQuick(i) * x.getQuick(i);
		}
		else {
			sum += y.getQuick(i) * x.getQuick(i);
			for (int j = i + 1; j < size; j++) {
				sum += A.getQuick(i,j) * x.getQuick(j);			}
		}
		b.setQuick(i,sum);
	}
	x.assign(b);
}
 

public FactorizationUserIntentModel(DoubleMatrix1D userVector) {
    Set<Integer> nonZeroFidx = new HashSet<>();
    for (int i = 0; i < userVector.size(); i++) {
        if (userVector.getQuick(i) > 0) {
            nonZeroFidx.add(i);
        }
    }
    this.userVector = userVector;
    this.nonZeroFactors = nonZeroFidx.stream()
            .map(featureData::fidx2feature)
            .collect(Collectors.toSet());
}
 

public GaussianLowEnergySampler(double thresh, ObjectiveFunction of, DoubleMatrix1D DOFmin, 
        DoubleMatrix1D DOFmax, DoubleMatrix1D center) {
    this.EPICThresh1 = thresh;
    this.of = of;
    this.DOFmin = DOFmin;
    this.DOFmax = DOFmax;
    this.center = center;
    numDOFs = DOFmin.size();
    
    //OK we need to figure out sigmas based on thresh
    
    
    double sigmaVal = chooseNumSigmas();
    
    //now go along axes to find sigma
    sigmas = DoubleFactory1D.dense.make(numDOFs);
    for(int dofNum=0; dofNum<numDOFs; dofNum++){
                    
        //bound the good region along dofNum axis
        double goodRegionLB = bisectForGoodRegionEnd(dofNum, DOFmin.get(dofNum));
        double goodRegionUB = bisectForGoodRegionEnd(dofNum, DOFmax.get(dofNum));
        
        //sigma will based on the farther of these from center
        double sigma;
        if( goodRegionUB - center.get(dofNum) >= center.get(dofNum) - goodRegionLB ){
            sigma = (goodRegionUB - center.get(dofNum)) / sigmaVal;
        }
        else
            sigma = (center.get(dofNum) - goodRegionLB) / sigmaVal;
        
        sigmas.set(dofNum, sigma);
    }
}
 

@Override
public void setDOFs(DoubleMatrix1D x) {
               curDOFVals.assign(x);
	for (int d=0; d<x.size(); d++) {
		pmol.dofs.get(d).apply(x.get(d));
                       //handleBlocksTogetherMaybe();//DEBUG!!!
	}
}
 

@Override
public double train(List<Minimizer.Result> trainingSet, List<Minimizer.Result> testSet) {

	// make sure all the samples have the right dimension
	for (Minimizer.Result sample : trainingSet) {
		if (sample.dofValues.size() != numDofs) {
			throw new IllegalArgumentException("samples have wrong number of dimensions");
		}
	}

	LinearSystem trainingSystem = new LinearSystem(trainingSet);
	LinearSystem testSystem = new LinearSystem(testSet);

	// solve Ax = b in least squares sense
	coefficients.assign(new QRDecomposition(trainingSystem.A).solve(trainingSystem.b).viewColumn(0));

	// calculate the residual (Ax - b) for the test set
	DoubleMatrix1D residual = new Algebra()
		.mult(testSystem.A, coefficients)
		.assign(testSystem.b.viewColumn(0), (ri, bi) -> Math.abs(ri - bi));

	// calculate the max error
	maxe = 0.0;
	for (int i=0; i<residual.size(); i++) {
		assert (residual.get(i) >= 0.0);
		maxe = Math.max(maxe, residual.get(i));
	}

	return maxe;
}
 

/**
 * Demonstrates how to use this class.
 */
public static void demo2() {
// parameters
double[] values = {
	//5, 0.0, -0.0, -Double.NaN, Double.NaN, 0.0/0.0, Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY, Double.MIN_VALUE, Double.MAX_VALUE
	5, 0.0, -0.0, -Double.NaN, Double.NaN, 0.0/0.0, Double.MIN_VALUE, Double.MAX_VALUE , Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY
	//Double.MIN_VALUE, Double.MAX_VALUE //, Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY
};
//String[] formats =         {"%G", "%1.10G", "%f", "%1.2f", "%0.2e"};
String[] formats =         {"%G", "%1.19G"};


// now the processing
int size = formats.length;
DoubleMatrix1D matrix = new DenseDoubleMatrix1D(values);

String[] strings = new String[size];
//String[] javaStrings = new String[size];

for (int i=0; i<size; i++) {
	String format = formats[i];
	strings[i] = new Formatter(format).toString(matrix);
	for (int j=0; j<matrix.size(); j++) {
		System.out.println(String.valueOf(matrix.get(j)));
	}
}

System.out.println("original:\n"+new Formatter().toString(matrix));

for (int i=0; i<size; i++) {
	System.out.println("\nstring("+formats[i]+"):\n"+strings[i]);
}

}
 

private static Matrix getMatrixSqrt(Matrix M, Boolean invert) {
    DoubleMatrix2D S = new DenseDoubleMatrix2D(M.toComponents());
    RobustEigenDecomposition eigenDecomp = new RobustEigenDecomposition(S, 100);
    DoubleMatrix1D eigenValues = eigenDecomp.getRealEigenvalues();
    int dim = eigenValues.size();
    for (int i = 0; i < dim; i++) {
        double value = sqrt(eigenValues.get(i));
        if (invert) {
            value = 1 / value;
        }
        eigenValues.set(i, value);
    }

    DoubleMatrix2D eigenVectors = eigenDecomp.getV();
    for (int i = 0; i < dim; i++) {
        for (int j = 0; j < dim; j++) {
            eigenVectors.set(i, j, eigenVectors.get(i, j) * eigenValues.get(j));

        }
    }
    DoubleMatrix2D storageMatrix = new DenseDoubleMatrix2D(dim, dim);
    eigenVectors.zMult(eigenDecomp.getV(), storageMatrix, 1, 0, false, true);


    return new Matrix(storageMatrix.toArray());

}
 

public void dtrmv(boolean isUpperTriangular, boolean transposeA, boolean isUnitTriangular, DoubleMatrix2D A, DoubleMatrix1D x) {
	if (transposeA) {
		A = A.viewDice();
		isUpperTriangular = !isUpperTriangular;
	}
	
	Property.DEFAULT.checkSquare(A);
	int size = A.rows();
	if (size != x.size()) {
		throw new IllegalArgumentException(A.toStringShort() + ", " + x.toStringShort());
	}
	    
	DoubleMatrix1D b = x.like();
	DoubleMatrix1D y = x.like();
	if (isUnitTriangular) {
		y.assign(1);
	}
	else {
		for (int i = 0; i < size; i++) {
			y.setQuick(i, A.getQuick(i,i));
		}
	}
	
	for (int i = 0; i < size; i++) {
		double sum = 0;
		if (!isUpperTriangular) {
			for (int j = 0; j < i; j++) {
				sum += A.getQuick(i,j) * x.getQuick(j);
			}
			sum += y.getQuick(i) * x.getQuick(i);
		}
		else {
			sum += y.getQuick(i) * x.getQuick(i);
			for (int j = i + 1; j < size; j++) {
				sum += A.getQuick(i,j) * x.getQuick(j);			}
		}
		b.setQuick(i,sum);
	}
	x.assign(b);
}
 

/**
Sorts the vector into ascending order, according to the <i>natural ordering</i>.
The returned view is backed by this matrix, so changes in the returned view are reflected in this matrix, and vice-versa.
To sort ranges use sub-ranging views. To sort descending, use flip views ...
<p>
<b>Example:</b> 
<table border="1" cellspacing="0">
  <tr nowrap> 
	<td valign="top"><tt> 7, 1, 3, 1<br>
	  </tt></td>
	<td valign="top"> 
	  <p><tt> ==&gt; 1, 1, 3, 7<br>
		The vector IS NOT SORTED.<br>
		The new VIEW IS SORTED.</tt></p>
	</td>
  </tr>
</table>

@param vector the vector to be sorted.
@return a new sorted vector (matrix) view. 
		<b>Note that the original matrix is left unaffected.</b>
*/
public DoubleMatrix1D sort(final DoubleMatrix1D vector) {
	int[] indexes = new int[vector.size()]; // row indexes to reorder instead of matrix itself
	for (int i=indexes.length; --i >= 0; ) indexes[i] = i;

	IntComparator comp = new IntComparator() {  
		public int compare(int a, int b) {
			double av = vector.getQuick(a);
			double bv = vector.getQuick(b);
			if (av!=av || bv!=bv) return compareNaN(av,bv); // swap NaNs to the end
			return av<bv ? -1 : (av==bv ? 0 : 1);
		}
	};

	runSort(indexes,0,indexes.length,comp);

	return vector.viewSelection(indexes);
}
 

/**
 * Returns the one-norm of vector <tt>x</tt>, which is <tt>Sum(abs(x[i]))</tt>.
 */
public double norm1(DoubleMatrix1D x) {
	if (x.size()==0) return 0;
	return x.aggregate(cern.jet.math.Functions.plus,cern.jet.math.Functions.abs);
}
 

/**
 * Returns the one-norm of vector <tt>x</tt>, which is <tt>Sum(abs(x[i]))</tt>.
 */
public double norm1(DoubleMatrix1D x) {
	if (x.size()==0) return 0;
	return x.aggregate(cern.jet.math.Functions.plus,cern.jet.math.Functions.abs);
}
 

/**
Modifies the given vector <tt>A</tt> such that it is permuted as specified; Useful for pivoting.
Cell <tt>A[i]</tt> will go into cell <tt>A[indexes[i]]</tt>.
<p>
<b>Example:</b>
<pre>
Reordering
[A,B,C,D,E] with indexes [0,4,2,3,1] yields 
[A,E,C,D,B]
In other words A[0]<--A[0], A[1]<--A[4], A[2]<--A[2], A[3]<--A[3], A[4]<--A[1].

Reordering
[A,B,C,D,E] with indexes [0,4,1,2,3] yields 
[A,E,B,C,D]
In other words A[0]<--A[0], A[1]<--A[4], A[2]<--A[1], A[3]<--A[2], A[4]<--A[3].
</pre>

@param   A   the vector to permute.
@param   indexes the permutation indexes, must satisfy <tt>indexes.length==A.size() && indexes[i] >= 0 && indexes[i] < A.size()</tt>;
@param   work the working storage, must satisfy <tt>work.length >= A.size()</tt>; set <tt>work==null</tt> if you don't care about performance.
@return the modified <tt>A</tt> (for convenience only).
@throws	IndexOutOfBoundsException if <tt>indexes.length != A.size()</tt>.
*/
public DoubleMatrix1D permute(DoubleMatrix1D A, int[] indexes, double[] work) {
	// check validity
	int size = A.size();
	if (indexes.length != size) throw new IndexOutOfBoundsException("invalid permutation");

	/*
	int i=size;
	int a;
	while (--i >= 0 && (a=indexes[i])==i) if (a < 0 || a >= size) throw new IndexOutOfBoundsException("invalid permutation");
	if (i<0) return; // nothing to permute
	*/

	if (work==null || size > work.length) {
		work = A.toArray();
	}
	else {
		A.toArray(work);
	}
	for (int i=size; --i >= 0; ) A.setQuick(i, work[indexes[i]]);
	return A;
}
 

/**
Sorts the vector into ascending order, according to the order induced by the specified comparator.
The returned view is backed by this matrix, so changes in the returned view are reflected in this matrix, and vice-versa.
The algorithm compares two cells at a time, determinining whether one is smaller, equal or larger than the other.
To sort ranges use sub-ranging views. To sort descending, use flip views ...
<p>
<b>Example:</b>
<pre>
// sort by sinus of cells
DoubleComparator comp = new DoubleComparator() {
&nbsp;&nbsp;&nbsp;public int compare(double a, double b) {
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;double as = Math.sin(a); double bs = Math.sin(b);
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;return as < bs ? -1 : as == bs ? 0 : 1;
&nbsp;&nbsp;&nbsp;}
};
sorted = quickSort(vector,comp);
</pre>

@param vector the vector to be sorted.
@param c the comparator to determine the order.
@return a new matrix view sorted as specified.
		<b>Note that the original vector (matrix) is left unaffected.</b>
*/
public DoubleMatrix1D sort(final DoubleMatrix1D vector, final cern.colt.function.DoubleComparator c) {
	int[] indexes = new int[vector.size()]; // row indexes to reorder instead of matrix itself
	for (int i=indexes.length; --i >= 0; ) indexes[i] = i;

	IntComparator comp = new IntComparator() {  
		public int compare(int a, int b) {
			return c.compare(vector.getQuick(a), vector.getQuick(b));
		}
	};

	runSort(indexes,0,indexes.length,comp);

	return vector.viewSelection(indexes);
}
 

/**
Sorts the vector into ascending order, according to the order induced by the specified comparator.
The returned view is backed by this matrix, so changes in the returned view are reflected in this matrix, and vice-versa.
The algorithm compares two cells at a time, determinining whether one is smaller, equal or larger than the other.
To sort ranges use sub-ranging views. To sort descending, use flip views ...
<p>
<b>Example:</b>
<pre>
// sort by sinus of cells
DoubleComparator comp = new DoubleComparator() {
&nbsp;&nbsp;&nbsp;public int compare(double a, double b) {
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;double as = Math.sin(a); double bs = Math.sin(b);
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;return as < bs ? -1 : as == bs ? 0 : 1;
&nbsp;&nbsp;&nbsp;}
};
sorted = quickSort(vector,comp);
</pre>

@param vector the vector to be sorted.
@param c the comparator to determine the order.
@return a new matrix view sorted as specified.
		<b>Note that the original vector (matrix) is left unaffected.</b>
*/
public DoubleMatrix1D sort(final DoubleMatrix1D vector, final cern.colt.function.DoubleComparator c) {
	int[] indexes = new int[vector.size()]; // row indexes to reorder instead of matrix itself
	for (int i=indexes.length; --i >= 0; ) indexes[i] = i;

	IntComparator comp = new IntComparator() {  
		public int compare(int a, int b) {
			return c.compare(vector.getQuick(a), vector.getQuick(b));
		}
	};

	runSort(indexes,0,indexes.length,comp);

	return vector.viewSelection(indexes);
}
 

/**
2-d OLAP cube operator; Fills all cells of the given vectors into the given histogram.
If you use hep.aida.ref.Converter.toString(histo) on the result, the OLAP cube of x-"column" vs. y-"column" , summing the weights "column" will be printed.
For example, aggregate sales by product by region.
<p>
Computes the distinct values of x and y, yielding histogram axes that capture one distinct value per bin.
Then fills the histogram.
<p>
Example output:
<table>
<td class="PRE"> 
<pre>
Cube:
&nbsp;&nbsp;&nbsp;Entries=5000, ExtraEntries=0
&nbsp;&nbsp;&nbsp;MeanX=4.9838, RmsX=NaN
&nbsp;&nbsp;&nbsp;MeanY=2.5304, RmsY=NaN
&nbsp;&nbsp;&nbsp;xAxis: Min=0, Max=10, Bins=11
&nbsp;&nbsp;&nbsp;yAxis: Min=0, Max=5, Bins=6
Heights:
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;| X
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;| 0   1   2   3   4   5   6   7   8   9   10  | Sum 
----------------------------------------------------------
Y 5   |  30  53  51  52  57  39  65  61  55  49  22 |  534
&nbsp;&nbsp;4   |  43 106 112  96  92  94 107  98  98 110  47 | 1003
&nbsp;&nbsp;3   |  39 134  87  93 102 103 110  90 114  98  51 | 1021
&nbsp;&nbsp;2   |  44  81 113  96 101  86 109  83 111  93  42 |  959
&nbsp;&nbsp;1   |  54  94 103  99 115  92  98  97 103  90  44 |  989
&nbsp;&nbsp;0   |  24  54  52  44  42  56  46  47  56  53  20 |  494
----------------------------------------------------------
&nbsp;&nbsp;Sum | 234 522 518 480 509 470 535 476 537 493 226 |     
</pre>
</td>
</table>
@return the histogram containing the cube.
@throws IllegalArgumentException if <tt>x.size() != y.size() || y.size() != weights.size()</tt>.
*/
public static hep.aida.IHistogram2D cube(DoubleMatrix1D x, DoubleMatrix1D y, DoubleMatrix1D weights) {
	if (x.size() != y.size() || y.size() != weights.size()) throw new IllegalArgumentException("vectors must have same size");
	
	double epsilon = 1.0E-9;
	cern.colt.list.DoubleArrayList distinct = new cern.colt.list.DoubleArrayList();
	double[] vals = new double[x.size()];
	cern.colt.list.DoubleArrayList sorted = new cern.colt.list.DoubleArrayList(vals);

	// compute distinct values of x
	x.toArray(vals); // copy x into vals
	sorted.sort();
	cern.jet.stat.Descriptive.frequencies(sorted, distinct, null);
	// since bins are right-open [from,to) we need an additional dummy bin so that the last distinct value does not fall into the overflow bin
	if (distinct.size()>0) distinct.add(distinct.get(distinct.size()-1) + epsilon);
	distinct.trimToSize();
	hep.aida.IAxis xaxis = new hep.aida.ref.VariableAxis(distinct.elements());

	// compute distinct values of y
	y.toArray(vals);
	sorted.sort();
	cern.jet.stat.Descriptive.frequencies(sorted, distinct, null);
	// since bins are right-open [from,to) we need an additional dummy bin so that the last distinct value does not fall into the overflow bin
	if (distinct.size()>0) distinct.add(distinct.get(distinct.size()-1) + epsilon);
	distinct.trimToSize();
	hep.aida.IAxis yaxis = new hep.aida.ref.VariableAxis(distinct.elements());

	hep.aida.IHistogram2D histo = new hep.aida.ref.Histogram2D("Cube",xaxis,yaxis);
	return histogram(histo,x,y,weights);
}