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

public double eval( DoubleMatrix1D x, int[] indices ) {
    //The "regular" coordinates to be used have the indicated indices in x

    switch(type){
        case REGULAR:
            return x.get(indices[0]);
        case SUMSQ:
            double sum=0;
            for(int q:indices)
                sum += Math.pow(x.get(q), 2);
            return Math.sqrt(sum);
        case LINCOMB:
            double ans=0;
            for(int a=0; a<coeffs.length; a++)
                ans += coeffs[a]*x.get(indices[a]);
            return ans;
        default:
            throw new Error("Unrecognized generalized coordinate type: "+type);
    }
}
 

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

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 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 double getValForDOF(int d1, double val, DoubleMatrix1D x) {

	// linear term
	double v = coefficients.get(index1(d1));

	// quadratic terms
	for (int d2=0; d2<numDofs; d2++) {
		double x2 = x.get(d2);
		double c = coefficients.get(index2(d1, d2));
		v += c*x2;
	}

	v *= x.get(d1);

	// constant term
	v += coefficients.get(0);

	return v;
}
 

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

DoubleMatrix1D nextSample(){//Sample from the distribution
    DoubleMatrix1D ans = DoubleFactory1D.dense.make(numDOFs);
    for(int dofNum=0; dofNum<numDOFs; dofNum++){
        do {
            double x = center.get(dofNum) + sigmas.get(dofNum) * random.nextGaussian();
            ans.set(dofNum, x);
        } while( ans.get(dofNum) < DOFmin.get(dofNum) 
                || ans.get(dofNum) > DOFmax.get(dofNum) );
    }
    
    return ans;
}
 

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

public boolean valuesInRange(DoubleMatrix1D x){
    //Check if x is in range for this voxel
    for(int dof=0; dof<numDOFs; dof++){
        if(x.get(dof)>DOFmax.get(dof) || x.get(dof)<DOFmin.get(dof))
            return false;
    }
    return true;
}
 

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

double Q(DoubleMatrix1D pt1, DoubleMatrix1D pt2, DoubleMatrix1D scale){
    //(Unnormalized) probability of getting to pt2 from pt1, sampling from a normal distribution
    //with the given standard deviations in each direction (the "sampling scale")
    DoubleMatrix1D diff = pt2.copy().assign(pt1,Functions.minus);

    //ok now we just need a standard normal probability
    //total probability is product of single-variable probabiliies
    double prob = 1;
    for(int dof=0; dof<numDOFs; dof++){
        double v = diff.get(dof)/scale.get(dof);
        prob *= Math.exp(-v*v/2);//standard normal probability, unnormalized
    }
    
    return prob;
}
 

DoubleMatrix1D getAnglesInRange(DoubleMatrix1D u){
    //Given a vector of this minimizer's DOFs, identify those that are angles
    //add multiples of 360 if possible to get them in the [DOFmin,DOFmax] range
    //angles are assumed to be bounded
    //either on both max and min side or not at all
    //this is for preparing initial values
    
    //we will move each out-of-range angle to the first equivalent angle after DOFmin
    DoubleMatrix1D ans = u.copy();
    
    for(int dof=0; dof<numDOFs; dof++){
        
        if( objFcn.isDOFAngle(dof) ){
        
            if(ans.get(dof)<DOFmin.get(dof)-0.001){
                ans.set(dof,
                        DOFmin.get(dof) + (ans.get(dof)-DOFmin.get(dof))%(360) + 360 );
            }
            else if(ans.get(dof)>DOFmax.get(dof)+0.001){
                ans.set(dof,
                        DOFmin.get(dof) + (ans.get(dof)-DOFmin.get(dof))%(360) );
            }
        }
    }

    return ans;
}
 

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

protected void boundEigenvalues(DoubleMatrix1D eigenvalues) {

                for (int i = 0; i < eigenvalues.cardinality(); ++i) {
                    if (eigenvalues.get(i) > MAX_EIGENVALUE) {
                        eigenvalues.set(i, MAX_EIGENVALUE);
                    } else if (eigenvalues.get(i) < MIN_EIGENVALUE) {
                        eigenvalues.set(i, MIN_EIGENVALUE);
                    }
                }
            }
 

static double[] evalSeriesByDegree(double[] coeffs, DoubleMatrix1D x, int nd, boolean includeConst, int order) {
    //like evalSeries, but return terms of each degree in the polynomial separately
    //(so ans[i] is the i-th degree monomials)
    //this only goes up to 4

    if(order>4){
        throw new Error("SeriesFitter.evalSeriesByDegree doesn't support order "+order);
    }
    
    int count = 0;//for counting through coeffs

    int topDegree = order;

    double ans[] = new double[topDegree+1];

    if(includeConst){
        ans[0] = coeffs[0];
        count++;
    }

    for(int dof=0; dof<nd; dof++){
        ans[1] += coeffs[count]*x.get(dof);
        count++;
    }

    for(int dof=0; dof<nd; dof++){
        for(int dof2=0; dof2<dof; dof2++){
            ans[2] += coeffs[count]*x.get(dof)*x.get(dof2);
            count++;
        }
        ans[2] += coeffs[count]*x.get(dof)*x.get(dof);///2;
        count++;
    }


    if(order>=3){

        for(int dof=0; dof<nd; dof++){
            for(int dof2=0; dof2<=dof; dof2++){
                for(int dof3=0; dof3<=dof2; dof3++){

                    ans[3] += coeffs[count]*x.get(dof)*x.get(dof2)*x.get(dof3);
                            //*distingPerm(dof,dof2,dof3)/6;
                    count++;
                }
            }
        }
    }

    if(order>=4){

        for(int dof=0; dof<nd; dof++){
            for(int dof2=0; dof2<=dof; dof2++){
                for(int dof3=0; dof3<=dof2; dof3++){
                    for(int dof4=0; dof4<=dof3; dof4++){
                        ans[4] += coeffs[count]*x.get(dof)*x.get(dof2)*x.get(dof3)*x.get(dof4);
                                //*distingPerm(dof,dof2,dof3,dof4)/24;
                        count++;
                    }
                }
            }
        }
    }

    return ans;
}
 

public static boolean bipartiteMatchingDepth(SparseDoubleMatrix2D candidateList, int[] signatureDepths, int[] functionDepths, int depth) {
	UndirectedGraph<String, DefaultEdge> g = new SimpleGraph<String, DefaultEdge>(DefaultEdge.class);
	List<String> signatureIdcs = new ArrayList<String>();
	int signatureNodesAtDepth = 0;
	for(int i=0; i<signatureDepths.length; ++i) {
		if(signatureDepths[i] == depth) {
			signatureIdcs.add("s"+i);
			g.addVertex("s"+i);
			//System.out.println("bpm:\ts"+i);
			signatureNodesAtDepth++;
		}
	}

	List<String> functionIdcs = new ArrayList<String>();
	for(int j=0; j<functionDepths.length; ++j) {
		if(functionDepths[j] == depth) {
			functionIdcs.add("f"+j);
			g.addVertex("f"+j);
			//System.out.println("bpm:\tf"+j);
		}
	}

	for(int i=0; i<signatureDepths.length; ++i) {
		if(signatureDepths[i] == depth) {
			DoubleMatrix1D row = candidateList.viewRow(i);

			for(int j=0; j<row.size(); ++j) {
				if(row.get(j) == 1.0 && functionDepths[j] == depth) {
					g.addEdge("s"+i, "f"+j);
				}
			}
		}
	}

	Set<String> p1 = new HashSet<String>(signatureIdcs);
       Set<String> p2 = new HashSet<String>(functionIdcs);

	HopcroftKarpBipartiteMatching<String, DefaultEdge> alg = 
		new HopcroftKarpBipartiteMatching<String, DefaultEdge>(g, p1, p2);

	Set<DefaultEdge> match = alg.getMatching();

	// System.out.println(g.toString());
	// System.out.println(match);
	if(match.size() == signatureNodesAtDepth)
		return true;
	else
		return false;

}
 

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

	int n = f.getNumDOFs();
	DoubleMatrix1D herex = startx.copy();
	DoubleMatrix1D nextx = startx.copy();

	// ccd is pretty simple actually
	// just do a line search along each dimension until we stop improving
	// we deal with cycles by just capping the number of iterations
	
	// get the current objective function value
	double herefx = f.getValue(herex);
	
	for (int iter=0; iter<MaxIterations; iter++) {
		
		// update all the dofs using line search
		for (int d=0; d<n; d++) {
			
			LineSearcher lineSearcher = lineSearchers.get(d);
			if (lineSearcher != null) {
				
				// get the next x value for this dof
				double xd = nextx.get(d);
				xd = lineSearcher.search(xd);
				nextx.set(d, xd);
			}
		}
		
		// how much did we improve?
		double nextfx = f.getValue(nextx);
		double improvement = herefx - nextfx;
		
		if (improvement > 0) {
			
			// take the step
			herex.assign(nextx);
			herefx = nextfx;
			
			if (improvement < ConvergenceThreshold) {
				break;
			}
			
		} else {
			break;
		}
	}

	// update the protein conf, one last time
	f.setDOFs(herex);

	return new Minimizer.Result(herex, herefx);
}
 

public static boolean bipartiteMatchingFull(SparseDoubleMatrix2D candidateList, int[] signatureDepths, int[] functionDepths, int depth) {
	UndirectedGraph<String, DefaultEdge> g = new SimpleGraph<String, DefaultEdge>(DefaultEdge.class);
	
	List<String> signatureIdcs = new ArrayList<String>();
	for(int i=0; i<signatureDepths.length; ++i) {
		signatureIdcs.add("s"+i);
		g.addVertex("s"+i);
	}

	List<String> functionIdcs = new ArrayList<String>();
	for(int j=0; j<functionDepths.length; ++j) {
		functionIdcs.add("f"+j);
		g.addVertex("f"+j);
	}

	for(int i=0; i<signatureDepths.length; ++i) {
		DoubleMatrix1D row = candidateList.viewRow(i);

		for(int j=0; j<row.size(); ++j) {
			if(row.get(j) != 0) {
				g.addEdge("s"+i, "f"+j);
			}
		}
	}

	Set<String> p1 = new HashSet<String>(signatureIdcs);
       Set<String> p2 = new HashSet<String>(functionIdcs);

	HopcroftKarpBipartiteMatching<String, DefaultEdge> alg = 
		new HopcroftKarpBipartiteMatching<String, DefaultEdge>(g, p1, p2);

	Set<DefaultEdge> match = alg.getMatching();

	//System.out.println(g.toString());
	//System.out.println(match);
	if(match.size() == signatureDepths.length)
		return true;
	else
		return false;

}
 

private void printFitTests(EPICFitter fitter, RCTuple RCList, double minEnergy,
        MoleculeModifierAndScorer mof, DoubleMatrix1D bestDOFVals, ArrayList<EPoly> series){
    //Do some tests on fit performance, and print the results
    int numDOFs = fitter.numDOFs;

    //TESTING FITS
    System.out.println("RCs: "+RCList.stringListing());
    System.out.println("Minimum energy: "+minEnergy);

    double testScales[] = new double[] { 0.01, 0.5, 5, 100 };//100
    int samplesPerScale = 3;



    double relMax[] = new double[numDOFs];//maximum shifts of degrees of freedom relative to minimum point (startVec)
    double relMin[] = new double[numDOFs];
    DoubleMatrix1D constr[] = mof.getConstraints();
    for(int dof=0; dof<numDOFs; dof++){
        relMax[dof] = constr[1].get(dof) - bestDOFVals.get(dof);
        relMin[dof] = constr[0].get(dof) - bestDOFVals.get(dof);
    }


    for(double scale : testScales){
        for(int s=0; s<samplesPerScale; s++){

            //Generate vector relative to minimum
            double dx[] = new double[numDOFs];
            //and absolute
            DoubleMatrix1D sampAbs = DoubleFactory1D.dense.make(numDOFs);
            for(int dof=0; dof<numDOFs; dof++){
                double top = Math.min(relMax[dof], scale);
                double bottom = Math.max(relMin[dof], -scale);
                dx[dof] = bottom + Math.random()*(top-bottom);
                sampAbs.set(dof, bestDOFVals.get(dof)+dx[dof]);
            }

            double trueVal = mof.getValue(sampAbs) - minEnergy;

            System.out.print("TEST: scale="+scale+" dx=");
            for(int dof=0; dof<numDOFs; dof++)
                System.out.print(dx[dof]+" ");

            System.out.print("TRUE="+trueVal+" FIT=");

            for(EPoly b : series){
                if(b!=null)
                    System.out.print(b.evaluate(sampAbs,false,false)+" ");
            }

            System.out.println();
        }
    }
}
 

private void printFitTests(EPICFitter fitter, RCTuple RCList, double minEnergy,
        MoleculeModifierAndScorer mof, DoubleMatrix1D bestDOFVals, ArrayList<EPoly> series){
    //Do some tests on fit performance, and print the results
    int numDOFs = fitter.numDOFs;

    //TESTING FITS
    System.out.println("RCs: "+RCList.stringListing());
    System.out.println("Minimum energy: "+minEnergy);

    double testScales[] = new double[] { 0.01, 0.5, 5, 100 };//100
    int samplesPerScale = 3;



    double relMax[] = new double[numDOFs];//maximum shifts of degrees of freedom relative to minimum point (startVec)
    double relMin[] = new double[numDOFs];
    DoubleMatrix1D constr[] = mof.getConstraints();
    for(int dof=0; dof<numDOFs; dof++){
        relMax[dof] = constr[1].get(dof) - bestDOFVals.get(dof);
        relMin[dof] = constr[0].get(dof) - bestDOFVals.get(dof);
    }


    for(double scale : testScales){
        for(int s=0; s<samplesPerScale; s++){

            //Generate vector relative to minimum
            double dx[] = new double[numDOFs];
            //and absolute
            DoubleMatrix1D sampAbs = DoubleFactory1D.dense.make(numDOFs);
            for(int dof=0; dof<numDOFs; dof++){
                double top = Math.min(relMax[dof], scale);
                double bottom = Math.max(relMin[dof], -scale);
                dx[dof] = bottom + Math.random()*(top-bottom);
                sampAbs.set(dof, bestDOFVals.get(dof)+dx[dof]);
            }

            double trueVal = mof.getValue(sampAbs) - minEnergy;

            System.out.print("TEST: scale="+scale+" dx=");
            for(int dof=0; dof<numDOFs; dof++)
                System.out.print(dx[dof]+" ");

            System.out.print("TRUE="+trueVal+" FIT=");

            for(EPoly b : series){
                if(b!=null)
                    System.out.print(b.evaluate(sampAbs,false,false)+" ");
            }

            System.out.println();
        }
    }
}