cern.colt.matrix.DoubleMatrix1D Java Examples

@Override
public double getEnergy() {
    if(curDOFVals==null){
        throw new RuntimeException("ERROR: Trying to evaluate an EPICEnergyFunction "
                + "before assigning it to a vector of DOF values");
    }
    
    double E = 0;
    for(int termNum=0; termNum<terms.size(); termNum++){
        EPoly term = terms.get(termNum);
        
        DoubleMatrix1D DOFValsForTerm = DoubleFactory1D.dense.make(term.numDOFs);
        for(int DOFCount=0; DOFCount<term.numDOFs; DOFCount++)
            DOFValsForTerm.set( DOFCount, curDOFVals.get(termDOFs.get(termNum).get(DOFCount)) );
        
        double termVal = term.evaluate(DOFValsForTerm, includeMinE, useSharedMolec);
        E += termVal;
    }
    
    return E;
}
 

/**
Assigns the result of a function to each cell; <tt>x[i] = function(x[i])</tt>.
(Iterates downwards from <tt>[size()-1]</tt> to <tt>[0]</tt>).
<p>
<b>Example:</b>
<pre>
// change each cell to its sine
matrix =   0.5      1.5      2.5       3.5 
matrix.assign(cern.jet.math.Functions.sin);
-->
matrix ==  0.479426 0.997495 0.598472 -0.350783
</pre>
For further examples, see the <a href="package-summary.html#FunctionObjects">package doc</a>.

@param function a function object taking as argument the current cell's value.
@return <tt>this</tt> (for convenience only).
@see cern.jet.math.Functions
*/
public DoubleMatrix1D assign(cern.colt.function.DoubleFunction function) {
	int s=stride;
	int i=index(0);
	double[] elems = this.elements;
	if (elems==null) throw new InternalError();

	// specialization for speed
	if (function instanceof cern.jet.math.Mult) { // x[i] = mult*x[i]
		double multiplicator = ((cern.jet.math.Mult)function).multiplicator;
		if (multiplicator==1) return this;
		for (int k=size; --k >= 0; ) {
			elems[i] *= multiplicator;
			i += s;
		}
	}
	else { // the general case x[i] = f(x[i])
		for (int k=size; --k >= 0; ) {
			elems[i] = function.apply(elems[i]);
			i += s;
		}
	}
	return this;
}
 

static void updateSeriesHessian(DoubleMatrix2D hess, DoubleMatrix1D x, double coeff, int...dofs ){
    //update the Hessian hess of a series at DOF values x
    //by adding in the Hessian of the term
    //coeff*prod_i x.get(dofs[i])
    int deg = dofs.length;
    
    for(int d1=0; d1<deg; d1++){
        for(int d2=0; d2<d1; d2++){//non diagonal Hessian contributions
            double dhess = coeff;
            
            for(int d3=0; d3<deg; d3++){
                if(d3!=d1 && d3!=d2)
                    dhess *= x.get(dofs[d3]);
            }
            
            hess.set(dofs[d1],dofs[d2], hess.get(dofs[d1],dofs[d2])+dhess);
            hess.set(dofs[d2],dofs[d1], hess.get(dofs[d2],dofs[d1])+dhess);
        }
    }
}
 

@Override
  public Minimizer.Result minimizeFrom(DoubleMatrix1D x) {

  	this.singleInitVal = x;

//First figure out the constraints and initial values
//(since the minimizer might be used for many rotameric states, etc.,
//this can't be done in the constructor)
DoubleMatrix1D constr[] = objFcn.getConstraints();
DOFmin = constr[0];
DOFmax = constr[1];

if(!compInitVals())//Initialize x
	return null;//No initial values found (this is currently only for IBEX)

      minimizeFromCurPoint();
      return new Minimizer.Result(this.x, objFcn.getValue(this.x));
  }
 

private static DoubleMatrix2D getGt(final DenseDoubleMatrix2D p, final DenseDoubleMatrix2D q, double lambda) {
    final int K = p.columns();

    DenseDoubleMatrix2D A1 = new DenseDoubleMatrix2D(K, K);
    q.zMult(q, A1, 1.0, 0.0, true, false);
    for (int k = 0; k < K; k++) {
        A1.setQuick(k, k, lambda + A1.getQuick(k, k));
    }

    EigenvalueDecomposition eig = new EigenvalueDecomposition(A1);
    DoubleMatrix1D d = eig.getRealEigenvalues();
    DoubleMatrix2D gt = eig.getV();
    for (int k = 0; k < K; k++) {
        double a = sqrt(d.get(k));
        gt.viewColumn(k).assign(x -> a * x);
    }

    return gt;
}
 

SampleNormalization(ArrayList<DoubleMatrix1D> fullSamples){
    center = DoubleFactory1D.dense.make(numDOFs);
    scaling = DoubleFactory1D.dense.make(numDOFs);
    jacDet = 1;
    
    for(int dofNum=0; dofNum<numDOFs; dofNum++){
        ArrayList<Double> vals = new ArrayList<>();
        int numSamples = fullSamples.size();
        double cen = 0;
        for(DoubleMatrix1D samp : fullSamples){
            vals.add(samp.get(dofNum));
            cen += samp.get(dofNum);
        }
        center.set(dofNum, cen/numSamples);
        Collections.sort(vals);
        //estimate spread using interquartile range
        double sc = vals.get(3*numSamples/4) - vals.get(numSamples/4);
        jacDet *= sc;
        scaling.set(dofNum, sc);
    }
}
 

static DoubleMatrix2D QuQt(DoubleMatrix2D Q, DoubleMatrix1D u){
    //Return matrix product of Q * diagonal version of u * Q^T
    //answer = \sum_i ( u_i * q_i * q_i')  is a (d+1)x(d+1) matrix

    int m = Q.rows();
    int n = Q.columns();
    if(u.size()!=n){
        throw new Error("Size mismatch in QuQt: "+n+" columns in Q, u length="+u.size());
    }

    DoubleMatrix2D ans = DoubleFactory2D.dense.make(m,m);

    for(int i=0; i<m; i++){
        for(int j=0; j<m; j++){
            double s = 0;
            for(int k=0; k<n; k++)
                s += Q.getQuick(i,k)*Q.getQuick(j,k)*u.get(k);
            ans.setQuick(i, j, s);
        }
    }

    return ans;
}
 

/**
 */
public static void doubleTest28(DoubleFactory2D f) {
	double[] data={1,2,3,4,5,6};
	double[][] arrMatrix = 
	{ 
		{ 1, 2, 3, 4, 5, 6},
		{ 2, 3, 4, 5, 6, 7}
	};
	
	DoubleMatrix1D vector = new DenseDoubleMatrix1D(data);
	DoubleMatrix2D matrix = f.make(arrMatrix);
	DoubleMatrix1D res = vector.like(matrix.rows());
	
	matrix.zMult(vector,res);

	System.out.println(res);
}
 

/**
 * Returns the row of the item matrix corresponding to the given item.
 *
 * @param i item
 * @return row of the item matrix
 */
public DoubleMatrix1D getItemVector(I i) {
    int iidx = item2iidx(i);
    if (iidx < 0) {
        return null;
    } else {
        return itemMatrix.viewRow(iidx);
    }
}
 

/**
 * Returns whether all cells of the given matrix <tt>A</tt> are equal to the given value.
 * The result is <tt>true</tt> if and only if <tt>A != null</tt> and
 * <tt>! (Math.abs(value - A[i]) > tolerance())</tt> holds for all coordinates.
 * @param   A   the first matrix to compare.
 * @param   value   the value to compare against.
 * @return  <tt>true</tt> if the matrix is equal to the value;
 *          <tt>false</tt> otherwise.
 */
public boolean equals(DoubleMatrix1D A, double value) {
	if (A==null) return false;
	double epsilon = tolerance();
	for (int i = A.size(); --i >= 0;) {
		//if (!(A.getQuick(i) == value)) return false;
		//if (Math.abs(value - A.getQuick(i)) > epsilon) return false;
		double x = A.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;
}
 

public VoxelSeriesChecker(List<Residue> residues, int numFreeDOFs, int numFullDOFs, 
        double[][] fullDOFPolys, PepPlaneLinModel[] pepPlanes, DoubleMatrix2D freeDOFMatrix,
        DoubleMatrix1D freeDOFCenter
) {
    //to be called in init conf!!
    numRes = residues.size();
    this.numFreeDOFs = numFreeDOFs;
    this.numFullDOFs = numFullDOFs;
    this.fullDOFPolys = fullDOFPolys;
    this.pepPlanes = pepPlanes;
    this.freeDOFMatrix = freeDOFMatrix;
    this.freeDOFCenter = freeDOFCenter;
    
    NCoord = new double[numRes][];
    CACoord = new double[numRes][];
    CCoord = new double[numRes][];
    
    for(int resNum=0; resNum<numRes; resNum++){
        NCoord[resNum] = residues.get(resNum).getCoordsByAtomName("N");
        CACoord[resNum] = residues.get(resNum).getCoordsByAtomName("CA");
        if(resNum==numRes-1)//use actual, instead of projected, C'
            CCoord[resNum] = residues.get(resNum).getCoordsByAtomName("C");
    }
    
    for(int resNum=0; resNum<numRes-1; resNum++){//make projected C' from N, CA
        CCoord[resNum]= pepPlanes[resNum].calcCCoords(CACoord[resNum], NCoord[resNum+1],
                CACoord[resNum+1], true);
    }
    
    targetConstraintVals = calcConstraintVals();
}
 

@Override
public double error(Factorization<U, I> factorization, FastPreferenceData<U, I> data) {
    DenseDoubleMatrix2D pu_z = factorization.getUserMatrix();
    DenseDoubleMatrix2D piz = factorization.getItemMatrix();

    return data.getUidxWithPreferences().parallel().mapToDouble(uidx -> {
        DoubleMatrix1D pU_z = pu_z.viewRow(uidx);
        DoubleMatrix1D pUi = piz.zMult(pU_z, null);
        return data.getUidxPreferences(uidx)
                .mapToDouble(iv -> -iv.v2 * pUi.getQuick(iv.v1))
                .sum();
    }).sum();

}
 

public EPolyPC( EPoly template, int fullOrder, int PCOrder, double PCFac) {
    //set up the principal component basis and the DOF information
    //coeffs will be added later
    //since we will be using this EPolyPC object when performing the fitting
    
    super(template.numDOFs,template.DOFmax,template.DOFmin,template.center,
            template.minE,null,fullOrder,template.DOFNames);
    
    //get the coordinate transformation into the eigenbasis of the template's Hessian
    //so we can use the high-eigenvalue directions as our principcal components
    DoubleMatrix2D hess = SeriesFitter.getHessian(template.coeffs,numDOFs,false);
    EigenvalueDecomposition edec = new EigenvalueDecomposition(hess);

    DoubleMatrix1D eigVals = edec.getRealEigenvalues();
    DoubleMatrix2D eigVecs = edec.getV();//The columns of eigVec are the eigenvectors

    invAxisCoeffs = eigVecs;//axisCoeffs' inverse is its transpose, since it's orthogonal
    axisCoeffs = Algebra.DEFAULT.transpose(eigVecs);
    

    this.fullOrder = fullOrder;
    this.PCOrder = PCOrder;


    //Now figure out what PC's to use
    //Let's try the criterion to use all PCs
    //whose eigenvalues' absolute values are within a factor of PCFac
    //of the biggest

    double maxEigVal = 0;
    for(int n=0; n<numDOFs; n++)
        maxEigVal = Math.max( Math.abs(eigVals.get(n)), maxEigVal );

    isPC = new boolean[numDOFs];

    for(int n=0; n<numDOFs; n++){
        if( Math.abs(eigVals.get(n)) >= maxEigVal*PCFac )
            isPC[n] = true;
    }
}
 

public int idamax(DoubleMatrix1D x) {
	int maxIndex = -1;
	double maxValue = Double.MIN_VALUE;
	for (int i=x.size(); --i >= 0; ) {
		double v = Math.abs(x.getQuick(i));
		if (v > maxValue) {
			maxValue = v;
			maxIndex = i;
		}
	}
	return maxIndex;
}
 

/**
 * Normalizes matrix of p(z|u) such that \forall_z: \sum_u p(z|u) = 1.
 *
 * @param pu_z normalized matrix of p(z|u)
 */
protected void normalizePuz(DoubleMatrix2D pu_z) {
    for (int z = 0; z < pu_z.columns(); z++) {
        final DoubleMatrix1D pu_Z = pu_z.viewColumn(z);
        pu_Z.assign(mult(1 / pu_Z.aggregate(plus, identity)));
    }
}
 

public boolean isOutOfRange(DoubleMatrix1D x){
    if(x.size()!=DOFs.size())
        throw new RuntimeException("ERROR: Trying to check range on "+DOFs.size()+" DOFs with "+x.size()+" values");
    
    for(int dof=0; dof<DOFs.size(); dof++){
        if( x.get(dof) < constraints[0].get(dof)-1e-6 )
            return true;
        if( x.get(dof) > constraints[1].get(dof)+1e-6 )
            return true;
    }
    
    return false;
}
 

public DoubleMatrix2D hessian(DoubleMatrix1D x) {
     
    DoubleMatrix1D z = toRelCoords(x);
    
    if(this instanceof EPolyPC)
        throw new RuntimeException("ERROR: Hessian for EPolyPC not currently supported");
    
    DoubleMatrix2D hess = seriesHessian(z);
             
    if(sapeTerm!=null)
        throw new RuntimeException("ERROR: SVE Hessian not currently supported");
    else 
        return hess;
}
 

private boolean checkStaysPositive(EPoly term, MoleculeModifierAndScorer mof){
    //Check, by minimization, that this EPIC term doesn't go too far negative
    //mof defines the voxel for the term conveniently
    
    ArrayList<EPoly> termList = new ArrayList<>();
    termList.add(term);
    EPICEnergyFunction epicEF = new EPICEnergyFunction(termList, false);
    
    ObjectiveFunction ofEPIC = new MoleculeModifierAndScorer( epicEF, mof.getConstraints(),
            mof.getMolec(), mof.getDOFs() );
    
    CCDMinimizer emin = new CCDMinimizer(ofEPIC, true);
    DoubleMatrix1D lowPoint = emin.minimize().dofValues;
    
    //energies will be relative to the voxel center energy (from the initial minimization)
    double lowPointTrueE = mof.getValue(lowPoint) - term.getMinE();
    double lowPointEPICE = ofEPIC.getValue(lowPoint);
    
    
    double tol = 0.1;//we'll only consider this a problem if the minimum
    //(where EPIC should be pretty accurate) is well below the actual minimum
    if(lowPointEPICE < lowPointTrueE - tol){
        System.out.println("Rejecting fit because of low minimum for EPIC term, "
                +lowPointEPICE+".  Corresponding true E: "+lowPointTrueE);
        return false;
    }
    
    //But let's warn if the minimum is low in absolute terms
    if(lowPointEPICE < -tol){
        System.out.println("Warning: low minimum for EPIC term, "
                +lowPointEPICE+".  Corresponding true E: "+lowPointTrueE);
    }
    
    epicEF.unassignSharedMolec();
    
    return true;
}
 

/**
 * Sets all cells to the state specified by <tt>value</tt>.
 * @param    value the value to be filled into the cells.
 * @return <tt>this</tt> (for convenience only).
 */
public DoubleMatrix1D assign(double value) {
	int index = index(0);
	int s = this.stride;
	double[] elems = this.elements;
	for (int i=size; --i >= 0; ) {
		elems[index] = value;
		index += s;
	}
	return this;
}
 

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);
}
 

boolean checkValidPt(DoubleMatrix1D y){
    //is y in the voxel and below the threshold?  (The distribution we're sampling is
    //uniform over the space of valid pts and 0 elsewhere)
    
    //first check voxel bounds...this is cheap
    for(int dof=0; dof<numDOFs; dof++){
        if(y.get(dof)<DOFmin.get(dof) || y.get(dof)>DOFmax.get(dof))
            return false;
    }
    
    if(of.getValue(y)>thresh)
        return false;
    
    return true;
}
 

public boolean isInBounds(DoubleMatrix1D x) {
	for (int d=0; d<size(); d++) {
		if (!isInBounds(d, x.get(d))) {
			return false;
		}
	}
	return true;
}
 

protected void scaleEigenvalues(DoubleMatrix1D eigenvalues) {
    double sum = 0.0;
    for (int i = 0; i < eigenvalues.cardinality(); ++i) {
        sum += eigenvalues.get(i);
    }

    double mean = -sum / eigenvalues.cardinality();

    for (int i = 0; i < eigenvalues.cardinality(); ++i) {
        eigenvalues.set(i, eigenvalues.get(i) / mean);
    }
}
 

public SAPE (MoleculeModifierAndScorer objFcn, double distCutoff, DoubleMatrix1D[] sampAbs){
    //We'll initialize the SAPE for standalone work, and also
    //initialize information sufficient to reinitialize with a shared molecule
    
    
    //The standalone molecule and energy function for SAPE will belong only to the SAPE object
    //(we'll deep-copy, deleting fields not needed for SAPE)
    EnergyFunction mainEF = objFcn.getEfunc();
    objFcn.setEfunc(null);
    mofStandalone = (MoleculeModifierAndScorer) ObjectIO.deepCopy(objFcn);
    objFcn.setEfunc(mainEF);
    
    
    Molecule standaloneMolec = mofStandalone.getMolec();
    //compressMolecule(standaloneMolec);//if needed...see comments below (compressSVE)
    
    ffParams = null;
    findFFParams(mainEF);//get forcefield params from some component of mainEF
    if(ffParams==null){
        throw new RuntimeException("ERROR: Trying to approximate energy function with SAPE but can't find any forcefield-type energies");
    }

    pickAtomPairs(mainEF,distCutoff,sampAbs);//select which atom-pairs to include in SAPE
    
    //now make a standalone energy function including only these pairs
    EnergyFunction standaloneEFcn = makeSparseEFcn( standaloneMolec );
    mofStandalone.setEfunc(standaloneEFcn);
}
 

/**
 * Returns a string <tt>s</tt> such that <tt>Object[] m = s</tt> is a legal Java statement.
 * @param matrix the matrix to format.
 */
public String toSourceCode(DoubleMatrix1D matrix) {
	Formatter copy = (Formatter) this.clone();
	copy.setPrintShape(false);
	copy.setColumnSeparator(", ");
	String lead  = "{";
	String trail = "};";
	return lead + copy.toString(matrix) + trail;
}
 

public void dger(double alpha, DoubleMatrix1D x, DoubleMatrix1D y, DoubleMatrix2D A) {
	cern.jet.math.PlusMult fun = cern.jet.math.PlusMult.plusMult(0);
	for (int i=A.rows(); --i >= 0; ) {
		fun.multiplicator = alpha * x.getQuick(i);
 		A.viewRow(i).assign(y,fun);
		
	}
}
 

private Minimizer.Result minimizeWithVdw(ParametricMolecule pmol, ResidueInteractions inters, DoubleMatrix1D x) {
	try (EnergyFunction efunc = cpuContext.efuncs.make(inters, pmol.mol)) {
		ResidueForcefieldEnergy.Vdw vdwEfunc = new ResidueForcefieldEnergy.Vdw((ResidueForcefieldEnergy)efunc);
		try (Minimizer minimizer = cpuContext.minimizers.make(new MoleculeObjectiveFunction(pmol, vdwEfunc))) {
			return minimizer.minimizeFrom(x);
		}
	}
}
 

/**
 * Returns whether all cells of the given matrix <tt>A</tt> are equal to the given value.
 * The result is <tt>true</tt> if and only if <tt>A != null</tt> and
 * <tt>! (Math.abs(value - A[i]) > tolerance())</tt> holds for all coordinates.
 * @param   A   the first matrix to compare.
 * @param   value   the value to compare against.
 * @return  <tt>true</tt> if the matrix is equal to the value;
 *          <tt>false</tt> otherwise.
 */
public boolean equals(DoubleMatrix1D A, double value) {
	if (A==null) return false;
	double epsilon = tolerance();
	for (int i = A.size(); --i >= 0;) {
		//if (!(A.getQuick(i) == value)) return false;
		//if (Math.abs(value - A.getQuick(i)) > epsilon) return false;
		double x = A.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;
}
 

/**
 * Normalizes matrix of p(z|u) such that \forall_u: \sum_z p(z|u) = 1.
 *
 * @param pu_z normalized matrix of p(z|u)
 */
@Override
protected void normalizePuz(DoubleMatrix2D pu_z) {
    for (int u = 0; u < pu_z.rows(); u++) {
        DoubleMatrix1D tmp = pu_z.viewRow(u);
        double norm = tmp.aggregate(Functions.plus, Functions.identity);
        if (norm != 0.0) {
            tmp.assign(Functions.mult(1 / norm));
        }
    }
}
 

public void dgemv(boolean transposeA, double alpha, DoubleMatrix2D A, DoubleMatrix1D x, double beta, DoubleMatrix1D y) {
	A.zMult(x,y,alpha,beta,transposeA);
}