Java Code Examples for cern.colt.matrix.DoubleMatrix2D#set()

/**
 */
public static void doubleTest() {
int rows = 4;
int columns = 5; // make a 4*5 matrix
DoubleMatrix2D master = new DenseDoubleMatrix2D(rows,columns);
System.out.println(master);
master.assign(1); // set all cells to 1
System.out.println("\n"+master);
master.viewPart(2,1,2,3).assign(2); // set [2,1] .. [3,3] to 2
System.out.println("\n"+master);

DoubleMatrix2D copyPart = master.viewPart(2,1,2,3).copy();
copyPart.assign(3); // modify an independent copy
copyPart.set(0,0,4); 
System.out.println("\n"+copyPart); // has changed
System.out.println("\n"+master); // master has not changed

DoubleMatrix2D view1 = master.viewPart(0,3,4,2); // [0,3] .. [3,4]
DoubleMatrix2D view2 = view1.viewPart(0,0,4,1); // a view from a view 
System.out.println("\n"+view1);
System.out.println("\n"+view2);
}
 

private void makeDistConstrJac( DoubleMatrix2D constrJac, int constrNum,
        double coord1[], double coord2[],
        int atomNum1, int atomNum2 ) {
    //Enter the Jacobian of a distance constraint between the atoms as #constrNum
    //in constrJac.  Atom numbers (among free N, CA) and (initial) coordinates given
    //-1 for atomNum means fixed atom, so not part of gradient
    
    for(int dim=0; dim<3; dim++){
        if(atomNum1>=0 && 3*atomNum1<constrJac.columns()){
            constrJac.set(constrNum, 3*atomNum1+dim, 2*coord1[dim]-2*coord2[dim]);
            
            recordJacDeriv(constrNum, 3*atomNum1+dim, 3*atomNum1+dim, 2);
            if(atomNum2>=0 && 3*atomNum2<constrJac.columns())//atomNum2 free to move
                recordJacDeriv(constrNum, 3*atomNum1+dim, 3*atomNum2+dim, -2);
        }
        
        if(atomNum2>=0 && 3*atomNum2<constrJac.columns()){
            constrJac.set(constrNum, 3*atomNum2+dim, 2*coord2[dim]-2*coord1[dim]);
            
            recordJacDeriv(constrNum, 3*atomNum2+dim, 3*atomNum2+dim, 2);
            if(atomNum1>=0 && 3*atomNum1<constrJac.columns())//atomNum1 free to move
                recordJacDeriv(constrNum, 3*atomNum2+dim, 3*atomNum1+dim, -2);
        }
    }
}
 

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

static DoubleMatrix2D getHessian(double coeffs[], int numDOFs, boolean includeConst){
    //Extract the Hessian from the given series coeffs, 
    //which are expected to have order>=2
    
    int count = numDOFs;//start at first index in coeffs for a second-order coefficient
    if(includeConst)
        count++;
    
    DoubleMatrix2D hess = DoubleFactory2D.dense.make(numDOFs,numDOFs);

    for(int a=0; a<numDOFs; a++){
        for(int b=0; b<=a; b++){
            double hessVal = coeffs[count];
            if(b==a)
                hessVal *= 2;
            hess.set(a, b, hessVal);
            hess.set(b, a, hessVal);
            count++;
        }
    }
    
    return hess;
}
 

private void calcTerms(DoubleMatrix2D terms){
	terms.assign(0.0);
	for(int x=0;x<windowLength;++x){
		for(int y=0;y<featureLength;++y){
			for(int i=0;i<k;++i){
				for(int j=0;j<l;++j){
					terms.set(featureLength*x+y,l*i+j,Math.pow(x,i)*Math.pow(y,j));
				}
			}
		}
	}
}
 

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 static void doubleTest(int rows, int columns, int initialCapacity, double minLoadFactor, double maxLoadFactor) {
	DoubleMatrix2D matrix = new SparseDoubleMatrix2D(rows,columns,initialCapacity,minLoadFactor,maxLoadFactor);
	System.out.println(matrix);
	
	System.out.println("adding...");
	int i=0;
	for (int column=0; column < columns; column++) {
		for (int row=0; row < rows; row++) {
			//if (i%1000 == 0) { 
 				matrix.set(row,column,i);
			//}
			i++;
		}
	}
	System.out.println(matrix);

	System.out.println("removing...");
	for (int column=0; column < columns; column++) {
		for (int row=0; row < rows; row++) {
			//if (i%1000 == 0) {
				matrix.set(row,column,0);
			//}
		}
	}

	System.out.println(matrix);
	System.out.println("bye bye.");
}
 

private void calcTerms(DoubleMatrix2D terms){
	terms.assign(0.0);
	for(int x=0;x<windowLength;++x){
		for(int y=0;y<featureLength;++y){
			for(int i=0;i<xDim;++i){
				for(int j=0;j<yDim;++j){
					terms.set(yDim*i+j,featureLength*x+y,Math.pow(x,i)*Math.pow(y,j));
				}
			}
		}
	}
}
 

static DoubleMatrix2D unflattenMatrix(double serCoeffs[], int nd, boolean removeLast, boolean coeffForm){
    //nd = number of dimensions
    
    int count=0;

    DoubleMatrix2D hess = DoubleFactory2D.dense.make(nd,nd);

    for(int a=0; a<nd; a++){
        for(int b=0; b<=a; b++){
            if(b==a){
                if(removeLast && a==nd-1)
                    hess.set(a,a,1);
                else
                    hess.set(a,a,serCoeffs[count]);
            }
            else{
                double co = serCoeffs[count];
                if(!coeffForm)
                    co /= 2;
                hess.set(a,b,co);//Hessian is symmetric
                hess.set(b,a,co);//but computing by Hessian double-counts off-diagonal series coefficients
            }

            count++;
        }
    }

    return hess;
}
 

private void makeLastNCACCConstrJac( DoubleMatrix2D constrJac, int constrNum,
        double CACoord[], double CCoord[], int atomNumN ){ 
    //Last constraint is simpler, because CA and C' are both fixed
    
    for(int dim=0; dim<3; dim++)
        constrJac.set(constrNum, 3*atomNumN+dim, CCoord[dim] - CACoord[dim]);
    
    //no jac derivs because this constraint is linear
}
 

private void calcTerms(DoubleMatrix2D terms){
	terms.assign(0.0);
	for(int x=0;x<windowLength;++x){
		for(int y=0;y<featureLength;++y){
			for(int i=0;i<k;++i){
				for(int j=0;j<l;++j){
					terms.set(l*i+j,featureLength*x+y,Math.pow(x,i)*Math.pow(y,j));
				}
			}
		}
	}
}
 

private void calcTerms(DoubleMatrix2D terms){
	terms.assign(0.0);
	for(int x=0;x<windowLength;++x){
		for(int y=0;y<featureLength;++y){
			for(int i=0;i<k;++i){
				for(int j=0;j<l;++j){
					terms.set(l*i+j,featureLength*x+y,Math.pow(x,i)*Math.pow(y,j));
				}
			}
		}
	}
}
 

static DoubleMatrix2D invertX(DoubleMatrix2D X){
        
        if(X.rows()!=X.columns())
            throw new RuntimeException("ERROR: Can't invert square matrix");
        
        if(X.rows()==1){//easy case of 1-D
            if(Math.abs(X.get(0,0))<1e-8)
                return DoubleFactory2D.dense.make(1,1,0);
            else
                return DoubleFactory2D.dense.make(1,1,1./X.get(0,0));
        }
    
        /*DoubleMatrix2D invX = null;
    
        try{
            invX = Algebra.DEFAULT.inverse(X);
        }
        catch(Exception e){
            System.err.println("ERROR: Singular matrix in MinVolEllipse calculation");
        }*/
        
        //Invert X, but if singular, this probably indicates a lack of points
        //but we can still make a meaningful lower-dimensional ellipsoid for the dimensions we have data
        //so we construct a pseudoinverse by setting those eigencomponents to 0
        //X is assumed to be symmetric (if not would need SVD instead...)
    
        EigenvalueDecomposition edec = new EigenvalueDecomposition(X);
        
        DoubleMatrix2D newD = edec.getD().copy();
        for(int m=0; m<X.rows(); m++){
            if( Math.abs(newD.get(m,m)) < 1e-8 ){
                //System.out.println("Warning: singular X in MinVolEllipse calculation");
                //this may come up somewhat routinely, especially with discrete DOFs in place
                newD.set(m,m,0);
            }
            else
                newD.set(m,m,1./newD.get(m,m));
        }

        DoubleMatrix2D invX = Algebra.DEFAULT.mult(edec.getV(),newD).zMult(edec.getV(), null, 1, 0, false, true);
        
        return invX;
}
 

private void makeNCACCConstrJac( DoubleMatrix2D constrJac, int constrNum,
        double NCoord[], double CACoord[], double NCoordNext[], double CACoordNext[],
        int atomNumN, int pepPlaneNum ){ 
    //constraint the dot product (C'-CA) dot (N-CA)
    //Uses N, CA coord variables for this residue and the next
    //(since C' coords are a linear function of this CA and next N, CA coords)
    //(For this purpose we constrain not the actual C' but its projection into the 
    //CA-N-C plane)
    
    //we'll need expansion coefficients for C'...
    double[] expansionCoeffs = pepPlanes[pepPlaneNum].getProjCAtomCoeffs();
    double a = expansionCoeffs[0];
    double b = expansionCoeffs[1];
    double c = expansionCoeffs[2];
    
    for(int dim=0; dim<3; dim++){
        if(atomNumN>=0){//not the first residue, which has fixed N & CA
            double NGrad = (a-1)*CACoord[dim] + b*NCoordNext[dim] + c*CACoordNext[dim];
            double CAGrad  = (a-1)*NCoord[dim] + 2*(1-a)*CACoord[dim] - b*NCoordNext[dim] - c*CACoordNext[dim];
            constrJac.set(constrNum, 3*atomNumN+dim, NGrad);
            constrJac.set(constrNum, 3*atomNumN+3+dim, CAGrad);
            
            //jac derivs for N
            recordJacDeriv(constrNum, 3*atomNumN+dim, 3*atomNumN+3+dim, a-1);
            recordJacDeriv(constrNum, 3*atomNumN+dim, 3*atomNumN+6+dim, b);
            if(3*atomNumN+9<constrJac.columns())
                recordJacDeriv(constrNum, 3*atomNumN+dim, 3*atomNumN+9+dim, c);
            
            //and for C
            recordJacDeriv(constrNum, 3*atomNumN+3+dim, 3*atomNumN+dim, a-1);
            recordJacDeriv(constrNum, 3*atomNumN+3+dim, 3*atomNumN+3+dim, 2*(1-a));
            recordJacDeriv(constrNum, 3*atomNumN+3+dim, 3*atomNumN+6+dim, -b);
            if(3*atomNumN+9<constrJac.columns())
                recordJacDeriv(constrNum, 3*atomNumN+3+dim, 3*atomNumN+9+dim, -c);
        }
        
        constrJac.set(constrNum, 3*atomNumN+6+dim, b*NCoord[dim]-b*CACoord[dim]);
        
        if(atomNumN>=0){
            recordJacDeriv(constrNum, 3*atomNumN+6+dim, 3*atomNumN+dim, b);
            recordJacDeriv(constrNum, 3*atomNumN+6+dim, 3*atomNumN+3+dim, -b);
        }
        
        
        if(3*atomNumN+9<constrJac.columns()){//i.e., not second-to-last res (which has fixed CANext)
            constrJac.set(constrNum, 3*atomNumN+9+dim, c*NCoord[dim]-c*CACoord[dim]);
            
            if(atomNumN>=0){
                recordJacDeriv(constrNum, 3*atomNumN+9+dim, 3*atomNumN+dim, c);
                recordJacDeriv(constrNum, 3*atomNumN+9+dim, 3*atomNumN+3+dim, -c);
            }
        }
    }
}
 

public static boolean bfsCompare(	SparseDoubleMatrix2D candidateList, 
									SparseDoubleMatrix2D signatureAdjacencyMatrix, 
									SparseDoubleMatrix2D functionAdjacencyMatrix, 
									int[] signatureDepths,
									int[] functionDepths,
									int depth) {
	// System.out.println("DEPTH: " + depth);
	// get vectors for nodes at this depth by retrieving parents (depth - 1)
	DoubleMatrix2D signatureVector = new SparseDoubleMatrix2D(1, signatureDepths.length);
	DoubleMatrix2D functionVector = new SparseDoubleMatrix2D(1, functionDepths.length);
	DoubleMatrix2D parentSigVector = new SparseDoubleMatrix2D(1, signatureDepths.length);
	DoubleMatrix2D parentFunVector = new SparseDoubleMatrix2D(1, functionDepths.length);

	// get parents
	for(int i=0; i<signatureDepths.length; ++i) {
		if(signatureDepths[i] == (depth-1))
			parentSigVector.set(0, i, 1.0);
	}

	for(int i=0; i<functionDepths.length; ++i) {
		if(functionDepths[i] == (depth-1))
			parentFunVector.set(0, i, 1.0);
	}

	// get current nodes
	signatureVector = Algebra.DEFAULT.mult(parentSigVector, signatureAdjacencyMatrix);
	functionVector = Algebra.DEFAULT.mult(parentFunVector, functionAdjacencyMatrix);

	// mark nodes' depths if they are not currently set
	for(int i=0; i<signatureDepths.length; ++i) {
		if(signatureVector.get(0,i) > 0 && signatureDepths[i] < 0)
			signatureDepths[i] = depth;
	}

	for(int i=0; i<functionDepths.length; ++i) {
		if(functionVector.get(0,i) > 0 && functionDepths[i] < 0)
			functionDepths[i] = depth;
	}

	// make signatureVector only have nodes at this depth
	for(int i=0; i<signatureDepths.length; ++i) {
		if(signatureDepths[i] != depth) {
			signatureVector.set(0, i, 0.0);
		}
	}

	// make signatureVector only have nodes at this depth
	for(int i=0; i<functionDepths.length; ++i) {
		if(functionDepths[i] != depth) {
			functionVector.set(0, i, 0.0);
		}
	}

	// check signode count <= funnode count
	if(	Algebra.DEFAULT.mult(signatureVector, signatureAdjacencyMatrix).cardinality() > 
		Algebra.DEFAULT.mult(functionVector, functionAdjacencyMatrix).cardinality()) {
		// System.out.println("Insufficient function basic blocks at depth " + depth + " to satisfy signature!");
		return false;
	}

	// TODO: edge check needed?
	
	// build candidate list 
	candidateList = buildCandidateList(	candidateList, 
										signatureAdjacencyMatrix, 
										functionAdjacencyMatrix, 
										signatureDepths, 
										functionDepths, 
										signatureVector, 
										functionVector, 
										depth);

	if(bipartiteMatchingVector(candidateList, signatureVector, functionVector, signatureDepths, functionDepths)) {
	//if(bipartiteMatchingDepth(candidateList, signatureDepths, functionDepths, depth)) {
		

		// check if we've reached end of graph
		if(Algebra.DEFAULT.mult(signatureVector, signatureAdjacencyMatrix).cardinality() == 0)
			return true;

	} else {
		// System.out.println("Bipartite matching failed!");
		return false;
	}

	if(depth > MAX_DEPTH)
		return false;

	return bfsCompare(candidateList, signatureAdjacencyMatrix, functionAdjacencyMatrix, signatureDepths, functionDepths, depth+1);
}
 

public static SparseDoubleMatrix2D checkIncestConditionParents(	DoubleMatrix2D signatureNodeIncestVector,
																DoubleMatrix2D functionNodeIncestVector,
																SparseDoubleMatrix2D signatureAdjacencyMatrix, 
																SparseDoubleMatrix2D functionAdjacencyMatrix,
																int[] signatureDepths,
																int[] functionDepths,
																int depth,
																SparseDoubleMatrix2D candidateList) {

	SparseDoubleMatrix2D signatureSelector = new SparseDoubleMatrix2D(1, signatureDepths.length);
	SparseDoubleMatrix2D functionSelector = new SparseDoubleMatrix2D(1, functionDepths.length);

	for(int i=0; i<signatureNodeIncestVector.size(); ++i) {
		if(signatureNodeIncestVector.get(0, i) == 0.0)
			continue;

		signatureSelector.assign(0.0);
		signatureSelector.set(0, i, 1.0);

		DoubleMatrix2D signatureNodeParentVector = Algebra.DEFAULT.mult(signatureSelector, Algebra.DEFAULT.transpose(signatureAdjacencyMatrix));

		// only bipartite match on incest parents
		for(int iParent = 0; iParent < signatureNodeParentVector.columns(); ++iParent) {
			if(signatureNodeParentVector.get(0, iParent) > 0 && signatureDepths[iParent] != depth) {
				signatureNodeParentVector.set(0, iParent, 0.0);
			}
		}

		int signatureIncestParentCount = signatureNodeParentVector.cardinality();

		for(int j=0; j<functionNodeIncestVector.size(); ++j) {
			if(candidateList.get(i, j) == 0.0)
				continue;

			functionSelector.assign(0.0);
			functionSelector.set(0, j, 1.0);

			DoubleMatrix2D functionNodeParentVector = Algebra.DEFAULT.mult(functionSelector, Algebra.DEFAULT.transpose(functionAdjacencyMatrix));

			// only bipartite match on incest parents
			for(int iParent = 0; iParent < functionNodeParentVector.columns(); ++iParent) {
				if(functionNodeParentVector.get(0, iParent) > 0 && functionDepths[iParent] != depth) {
					functionNodeParentVector.set(0, iParent, 0.0);
				}
			}

			int functionIncestParentCount = functionNodeParentVector.cardinality();

			// bipartite matching, if !sufficient
			if(	candidateList.get(i, j) > 0 && (signatureIncestParentCount > functionIncestParentCount ||
				!bipartiteMatchingVector(candidateList, signatureNodeParentVector, functionNodeParentVector, signatureDepths, functionDepths))) {
				// remove candidacy
				candidateList.set(i, j, 0.0);
			}
		}
	}

	return candidateList;
}
 

public static SparseDoubleMatrix2D checkPreviouslyVisitedChildren(	IntArrayList signatureNonZeros, 
																	IntArrayList functionNonZeros,
																	SparseDoubleMatrix2D signatureAdjacencyMatrix, 
																	SparseDoubleMatrix2D functionAdjacencyMatrix,
																	int[] signatureDepths,
																	int[] functionDepths,
 																	SparseDoubleMatrix2D candidateList) {
	
	SparseDoubleMatrix2D signatureSelector = new SparseDoubleMatrix2D(1, signatureDepths.length);
	SparseDoubleMatrix2D functionSelector = new SparseDoubleMatrix2D(1, functionDepths.length);

	// reduce candidacy based on previously visited children
	for(int i=0; i<signatureNonZeros.size(); ++i) {
		// for node N, this sets the Nth bit of the vector to 1.0
		signatureSelector.assign(0.0);
		signatureSelector.set(0, signatureNonZeros.get(i), 1.0);

		// get children at depth <= current
		DoubleMatrix2D signatureNodeChildVector = Algebra.DEFAULT.mult(signatureSelector, signatureAdjacencyMatrix);

		// remove signature node's unvisited children
		for(int iChild = 0; iChild < signatureNodeChildVector.columns(); ++iChild) {
			if(signatureNodeChildVector.get(0, iChild) > 0 && signatureDepths[iChild] == -1) { // for the oposite conditional && signatureDepths[iChild] <= depth) {
				signatureNodeChildVector.set(0, iChild, 0.0);
			}
		}

		int signatureVisitedChildCount = signatureNodeChildVector.cardinality();

		for(int j=0; j<functionNonZeros.size(); ++j) {
			functionSelector.assign(0.0);
			functionSelector.set(0, functionNonZeros.get(j), 1.0);
			// get children at depth <= current
			DoubleMatrix2D functionNodeChildVector = Algebra.DEFAULT.mult(functionSelector, functionAdjacencyMatrix);

			// remove function node's unvisited children
			for(int iChild = 0; iChild < functionNodeChildVector.columns(); ++iChild) {
				if(functionNodeChildVector.get(0, iChild) > 0 && functionDepths[iChild] == -1) { // for the oposite conditional && functionDepths[iChild] <= depth) {
					functionNodeChildVector.set(0, iChild, 0.0);
				}
			}

			int functionVisitedChildCount = functionNodeChildVector.cardinality();
			
			// bipartite matching, if !sufficient
			if(	candidateList.get(signatureNonZeros.get(i), functionNonZeros.get(j)) > 0 && 
				(signatureVisitedChildCount > functionVisitedChildCount ||
				!bipartiteMatchingVector(candidateList, signatureNodeChildVector, functionNodeChildVector, signatureDepths, functionDepths))) {
				// remove candidacy
				candidateList.set(signatureNonZeros.get(i), functionNonZeros.get(j), 0.0);
			}
		}
	}

	return candidateList;
}
 

public void normalize() {
	// System.out.println(adjacencyMatrix);
	int adjacencyMatrixSize = adjacencyMatrix.columns();
	
	if(adjacencyMatrixSize < 1)
		return;

	DoubleMatrix2D selector = new SparseDoubleMatrix2D(1, adjacencyMatrixSize);
	selector.set(0, 0, 1.0);
	// System.out.println(selector);
	// get vertices reachable from V0
	int cardinality = selector.cardinality();
	int lastCardinality = 0;

	DoubleDoubleFunction merge = new DoubleDoubleFunction() {
		public double apply(double a, double b) { 
			return (a > 0 || b > 0) ? 1 : 0; 
		}
	};

	while(cardinality != lastCardinality) {
		lastCardinality = cardinality;
		// selector = Algebra.DEFAULT.mult(selector, adjacencyMatrix);
		selector.assign(Algebra.DEFAULT.mult(selector, adjacencyMatrix), merge);
		// System.out.println(selector);
		cardinality = selector.cardinality();
	}

	// System.out.println(selector);

	if(cardinality == adjacencyMatrixSize) {
		return;
	}

	// IntArrayList nonZeros = new IntArrayList();
	// IntArrayList unusedInt = new IntArrayList();
	// DoubleArrayList unusedDouble = new DoubleArrayList();
	// selector.getNonZeros(unusedInt, nonZeros, unusedDouble);
	// SparseDoubleMatrix2D reduced = new SparseDoubleMatrix2D(adjacencyMatrix.viewSelection(nonZeros.elements(), nonZeros.elements()).toArray());

	SparseDoubleMatrix2D reduced = new SparseDoubleMatrix2D(cardinality, cardinality);
	int iCurrentVertex = 0;
	for(int i=0; i<adjacencyMatrixSize; ++i) {
		if(selector.get(0, i) != 0.0) {
			int jCurrentVertex = 0;
			for(int j=0; j<adjacencyMatrixSize; ++j) {
				if(selector.get(0, j) != 0.0){
					reduced.set(iCurrentVertex, jCurrentVertex, adjacencyMatrix.get(i, j));
					++jCurrentVertex;
				}
			}
			++iCurrentVertex;
		}
	}
	// System.out.println(reduced);
	// System.out.println("=======");
	adjacencyMatrix = reduced;
	vertexCount = adjacencyMatrix.columns();
	edgeCount = adjacencyMatrix.cardinality();
}
 

public HashMap<DegreeOfFreedom,Double> calcSpecialCenter(){
    //center for some degrees of freedom may, due to non-box constr, 
    //differ from center of box constr
    //NOTE: the map only contains those degrees of freedom where center is special
    //(not just middle of lb,ub)
    //DEBUG!!! This assumes the particular form of lin constr (pairs of parallel constr,
    //DOFs not in constrained space assumed to center at 0) used in basis big CATS stuff
    
    HashMap<DegreeOfFreedom,Double> ans = new HashMap<>();
    int numSpecialDOFConstr = voxLinConstr.size()/2;
    if(numSpecialDOFConstr==0)
        return ans;
    
    HashSet<Integer> specialDOFSet = new HashSet<>();
    for(int spesh=0; spesh<numSpecialDOFConstr; spesh++){
        RealVector coeff = voxLinConstr.get(2*spesh).getCoefficients();
        if(coeff.getDistance(voxLinConstr.get(2*spesh+1).getCoefficients()) > 0)
            throw new RuntimeException("ERROR: voxLinConstr not paired like expected");
        
        for(int d=0; d<dofIntervals.size(); d++){
            if(coeff.getEntry(d)!=0)
                specialDOFSet.add(d);
        }
    }
    
    int numSpecialDOFs = specialDOFSet.size();
    ArrayList<Integer> specialDOFList = new ArrayList<>();
    specialDOFList.addAll(specialDOFSet);
    
    DoubleMatrix2D specialConstr = DoubleFactory2D.dense.make(numSpecialDOFConstr, numSpecialDOFs);
    //we will constrain values first of the lin combs of DOFs given in voxLinConstr,
    //then of lin combs orth to those
    DoubleMatrix2D target = DoubleFactory2D.dense.make(numSpecialDOFs,1);
    for(int spesh=0; spesh<numSpecialDOFConstr; spesh++){
        RealVector coeffs = voxLinConstr.get(2*spesh).getCoefficients();
        for(int d=0; d<numSpecialDOFs; d++)
            specialConstr.set(spesh, d, coeffs.getEntry(specialDOFList.get(d)));
        target.set( spesh, 0, 0.5*(voxLinConstr.get(2*spesh).getValue()+voxLinConstr.get(2*spesh+1).getValue()) );
    }
    
    //remaining targets (orth components to voxLinConstr) are 0
    DoubleMatrix2D orthComponents = getOrthogVectors(specialConstr);
    DoubleMatrix2D M = DoubleFactory2D.dense.compose(
            new DoubleMatrix2D[][] {
                new DoubleMatrix2D[] {specialConstr},
                new DoubleMatrix2D[] {Algebra.DEFAULT.transpose(orthComponents)}
            }
        );
    
    DoubleMatrix1D specialDOFVals = Algebra.DEFAULT.solve(M, target).viewColumn(0);
    
    for(int d=0; d<numSpecialDOFs; d++){
        ans.put(dofIntervals.get(specialDOFList.get(d)).dof, specialDOFVals.get(d));
    }
    
    return ans;
}
 

private void computeStationaryDistribution(double[] stationaryDistribution) {

        if (allRatesAreZero(switchingRates)) {
            return;
        }

        // Uses an LU decomposition to solve Q^t \pi = 0 and \sum \pi_i = 1
        DoubleMatrix2D mat2 = new DenseDoubleMatrix2D(numBaseModel + 1, numBaseModel);
        int index2 = 0;
        for (int i = 0; i < numBaseModel; ++i) {
            for (int j = i + 1; j < numBaseModel; ++j) {
                mat2.set(j, i, switchingRates.getParameterValue(index2)); // Transposed
                index2++;
            }
        }
        for (int j = 0; j < numBaseModel; ++j) {
            for (int i = j + 1; i < numBaseModel; ++i) {
                mat2.set(j, i, switchingRates.getParameterValue(index2)); // Transposed
                index2++;
            }
        }
        for (int i = 0; i < numBaseModel; ++i) {
            double rowTotal = 0.0;
            for (int j = 0; j < numBaseModel; ++j) {
                if (i != j) {
                    rowTotal += mat2.get(j, i); // Transposed
                }

            }
            mat2.set(i, i, -rowTotal);
        }

        // Add row for sum-to-one constraint
        for (int i = 0; i < numBaseModel; ++i) {
            mat2.set(numBaseModel, i, 1.0);
        }

        LUDecomposition decomposition = new LUDecomposition(mat2);
        DoubleMatrix2D x = new DenseDoubleMatrix2D(numBaseModel + 1, 1);
        x.set(numBaseModel, 0, 1.0);
        DoubleMatrix2D y = decomposition.solve(x);
        for (int i = 0; i < numBaseModel; ++i) {
            stationaryDistribution[i] = y.get(i, 0);
        }
    }