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

/** 
Generates and returns the (economy-sized) orthogonal factor <tt>Q</tt>.
@return <tt>Q</tt>
*/
public DoubleMatrix2D getQ () {
	cern.jet.math.Functions F = cern.jet.math.Functions.functions;
	DoubleMatrix2D Q = QR.like();
	//double[][] Q = X.getArray();
	for (int k = n-1; k >= 0; k--) {
		DoubleMatrix1D QRcolk = QR.viewColumn(k).viewPart(k,m-k);
		Q.setQuick(k,k, 1);
		for (int j = k; j < n; j++) {
			if (QR.getQuick(k,k) != 0) {
				DoubleMatrix1D Qcolj = Q.viewColumn(j).viewPart(k,m-k);
				double s = QRcolk.zDotProduct(Qcolj);
				s = -s / QR.getQuick(k,k);
				Qcolj.assign(QRcolk, F.plusMult(s));
			}
		}
	}
	return Q;
}
 

/** 
Generates and returns the (economy-sized) orthogonal factor <tt>Q</tt>.
@return <tt>Q</tt>
*/
public DoubleMatrix2D getQ () {
	cern.jet.math.Functions F = cern.jet.math.Functions.functions;
	DoubleMatrix2D Q = QR.like();
	//double[][] Q = X.getArray();
	for (int k = n-1; k >= 0; k--) {
		DoubleMatrix1D QRcolk = QR.viewColumn(k).viewPart(k,m-k);
		Q.setQuick(k,k, 1);
		for (int j = k; j < n; j++) {
			if (QR.getQuick(k,k) != 0) {
				DoubleMatrix1D Qcolj = Q.viewColumn(j).viewPart(k,m-k);
				double s = QRcolk.zDotProduct(Qcolj);
				s = -s / QR.getQuick(k,k);
				Qcolj.assign(QRcolk, F.plusMult(s));
			}
		}
	}
	return Q;
}
 

double getScaledDistFromCenter(DoubleMatrix1D x){
    //How far are we from the center of the ellipse, relative to the boundary?
    
    if(x==null)
        System.out.println("woah...");
    
    DoubleMatrix1D xrel = x.copy().assign(nc,Functions.plus);
    return xrel.zDotProduct( Algebra.DEFAULT.mult(A, xrel) );
}
 

public double ddot(DoubleMatrix1D x, DoubleMatrix1D y) {
	return x.zDotProduct(y);
}
 

public double evalAtFullDOFs(DoubleMatrix1D fullDOFVals){
    return fullDOFVals.zDotProduct(coeffs);
}
 

double[] sampleResid(double[][] freeDOFVoxel){
    //draw a sample from the voxel and measure its constr resid and free-DOF resid
    
    //draw free DOFs for sample...
    DoubleMatrix1D sampFreeDOFs = DoubleFactory1D.dense.make(numFreeDOFs);
    
    for(int freeDOF=0; freeDOF<numFreeDOFs; freeDOF++){
        double voxWidth = freeDOFVoxel[1][freeDOF] - freeDOFVoxel[0][freeDOF];
        sampFreeDOFs.set( freeDOF, freeDOFVoxel[0][freeDOF] + Math.random()*voxWidth );
    }
    
    //OK now do full DOFs, placing them in NCoord and CACoord
    int fullDOFCount = 0;
    DoubleMatrix1D fullDOFVals = DoubleFactory1D.dense.make(numFullDOFs);
    
    for(int resNum=1; resNum<numRes; resNum++){
        
        for(int dim=0; dim<3; dim++){
            NCoord[resNum][dim] = evalFullDOF(fullDOFCount,sampFreeDOFs);
            fullDOFVals.set(fullDOFCount, NCoord[resNum][dim]);
            fullDOFCount++;
        }
        
        if(resNum==numRes-1)//no CA variables
            break;
        
        for(int dim=0; dim<3; dim++){
            CACoord[resNum][dim] = evalFullDOF(fullDOFCount,sampFreeDOFs);
            fullDOFVals.set(fullDOFCount, CACoord[resNum][dim]);
            fullDOFCount++;
        }
        
        
    }
    
    //Once N and CA in place, can calc C'.  Use plane projection, to match constr in jac
    for(int resNum=0; resNum<numRes-1; resNum++){
        CCoord[resNum]= pepPlanes[resNum].calcCCoords(CACoord[resNum], NCoord[resNum+1],
                CACoord[resNum+1], true);
    }
    
    //OK now handle add up constraint resids!
    ArrayList<Double> sampConstraintVals = calcConstraintVals();
    
    double constrResid = 0;
    for(int c=0; c<numFullDOFs-numFreeDOFs; c++){
        double dev = sampConstraintVals.get(c) - targetConstraintVals.get(c);
        constrResid += dev*dev;
    }
    
    constrResid /= (numFullDOFs-numFreeDOFs);//normalize resid by # of constraints
    
    
    DoubleMatrix1D freeDOFsCheck = Algebra.DEFAULT.mult(freeDOFMatrix, fullDOFVals);
    freeDOFsCheck.assign(freeDOFCenter, Functions.minus);
    freeDOFsCheck.assign(sampFreeDOFs, Functions.minus);//calc deviation in free DOFs
    double freeDOFResid = freeDOFsCheck.zDotProduct(freeDOFsCheck) / numFreeDOFs;
    
    return new double[] {constrResid, freeDOFResid};
}
 

double getScaling(DoubleMatrix1D x){
    //give the scaling factor needed to get x on the ellipse
    //(we are scaling x)
    double co1 = x.zDotProduct( Algebra.DEFAULT.mult(A, x) );
    double co2 = -2*x.zDotProduct( Algebra.DEFAULT.mult(A, c) );
    double co3 = c.zDotProduct( Algebra.DEFAULT.mult(A, c) ) - 1;
    
    if(co1==0 && co2==0 ){
        //this means basically scaling x has no effect on whether we're in the ellipse
        
        
        //because we use getScaling as a measure of how "far out" we are relative to the ellipse,
        //if we're outside (or basically on) the ellipse, we'll consider this to be a scaling of 0.  Otherwise infinity.
        
        if(co3>=-1e-10)//outside
            return 0;
        else
            return Double.POSITIVE_INFINITY;
        
        /*if(Math.abs(co3)<1e-10)//we're already basically on the ellipse: scaling x to 0 is sufficient
            return 0;
        
        System.out.print("Bad quadratic in getScaling! co's: "+co1+" "+co2+" "+co3+" A:");
        System.out.print(A);
        System.out.print("x: ");
        System.out.print(x);
        System.out.print("c: ");
        System.out.println(c);*/
    }
    
    double soln1 = quadraticFormula(co1,co2,co3,true,false);
    double soln2 = quadraticFormula(co1,co2,co3,false,false);
    
    //DEBUG!!
    if(Double.isNaN(co1)||Double.isNaN(co2)||Double.isNaN(co3)){
        System.out.print("NaN in getScaling! co's: "+co1+" "+co2+" "+co3+" A:");
        System.out.print(A);
        System.out.print("x: ");
        System.out.print(x);
        System.out.print("c: ");
        System.out.println(c);
    }
    
    
    if(soln1>0 && soln2>0){
        //System.err.println("ERROR: Ellipse doesn't contain origin!");//This may be fairly routine if origin is at edge...
        return Math.max(soln1,soln2);//when we finally get out of ellipse
    }
        
    if( (!(soln1>0)) && (!(soln2>0)) ){
        //System.err.println("ERROR: Ellipse doesn't contain origin!");//and our ray misses it!//May be a small error for small ellipses...
        return 0;
    }
    
    if(soln1>0)
        return soln1;
    else//(soln2>0)
        return soln2;
}
 

public double ddot(DoubleMatrix1D x, DoubleMatrix1D y) {
	return x.zDotProduct(y);
}
 

/**
 * Inner product of two vectors; <tt>Sum(x[i] * y[i])</tt>.
 * Also known as dot product.
 * <br>
 * Equivalent to <tt>x.zDotProduct(y)</tt>.
 *
 * @param x the first source vector.
 * @param y the second source matrix.
 * @return the inner product.
 *
 * @throws IllegalArgumentException if <tt>x.size() != y.size()</tt>.
 */
public double mult(DoubleMatrix1D x, DoubleMatrix1D y) {
	return x.zDotProduct(y);
}
 

/**
 * Inner product of two vectors; <tt>Sum(x[i] * y[i])</tt>.
 * Also known as dot product.
 * <br>
 * Equivalent to <tt>x.zDotProduct(y)</tt>.
 *
 * @param x the first source vector.
 * @param y the second source matrix.
 * @return the inner product.
 *
 * @throws IllegalArgumentException if <tt>x.size() != y.size()</tt>.
 */
public double mult(DoubleMatrix1D x, DoubleMatrix1D y) {
	return x.zDotProduct(y);
}