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

static double[] flattenMatrix(DoubleMatrix2D M, boolean removeLast, boolean coeffForm){
    //Given a symmetric matrix M
    //represent it in series form as in SeriesFitter, ConvexityConstraint, etc.
    //(and usable for maximization in maximizeQuadratic)
    //if removeLast, the lower right corner of M is assumed to be 1 and removed from the flattened form
    //if coeffForm, we are representing a linear function of a symmetric matrix, 
    //else just the matrix, so to make the functions line up, !coeffForm means double coefficients
    //as necessary so that flattenMatrix(A,removeLast,true) dot flattenMatrix(B,removeLast,false)
    // + (1 if removeLast) = Frobenius product(A,B)
    int count=0;

    int n = M.rows();
    int paramCount = SeriesFitter.getNumParamsForOrder(n, 2);
    if(removeLast)
        paramCount--;
    double ans[] = new double[paramCount];

    for(int a=0; a<n; a++){
        for(int b=0; b<=a; b++){
            if(b==a){
                ans[count] = M.get(a,a);
            }
            else if(!coeffForm)//M_ab and M_ba collapsed into one series coefficient (using symmetry)
                ans[count] = M.get(a,b)+M.get(b,a);           
            else
                ans[count] = M.get(a,b);

            count++;
            if(count==paramCount)//relevant for removeLast
                break;
        }
    }

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

private static double invMatrixDeriv(DoubleMatrix2D Minv, int i, int j, int u, int v){
    //d(M^-1)_ij / dM_uv
    return - Minv.get(i, u) * Minv.get(v, j);
}
 

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 static boolean bipartiteMatchingVector(DoubleMatrix2D candidateList, DoubleMatrix2D signatureNodeVector, DoubleMatrix2D functionNodeVector, int[] signatureDepths, int[] functionDepths) {
	UndirectedGraph<String, DefaultEdge> g = new SimpleGraph<String, DefaultEdge>(DefaultEdge.class);

	IntArrayList signatureNonZeros = new IntArrayList();
	IntArrayList functionNonZeros = new IntArrayList();
	IntArrayList unusedInt = new IntArrayList();
	DoubleArrayList unusedDouble = new DoubleArrayList();

	// get set column indices for signature vector and function vector
	signatureNodeVector.getNonZeros(unusedInt, signatureNonZeros, unusedDouble);
	functionNodeVector.getNonZeros(unusedInt, functionNonZeros, unusedDouble);
	
	List<String> signatureIdcs = new ArrayList<String>();
	List<String> functionIdcs = new ArrayList<String>();
	int signatureNodeCount = 0;
	// add signature nodes graph
	for(int i=0; i<signatureNonZeros.size(); ++i) {
		int signatureIdx = signatureNonZeros.get(i);
		if(signatureDepths[signatureIdx] != -1) {
			signatureIdcs.add("s"+signatureIdx);
			g.addVertex("s"+signatureIdx);
			signatureNodeCount++;
		}
	}

	// add function nodes graph
	for(int j=0; j<functionNonZeros.size(); ++j) {
		int functionIdx = functionNonZeros.get(j);
		if(functionDepths[functionIdx] != -1) {
			functionIdcs.add("f"+functionNonZeros.get(j));
			g.addVertex("f"+functionNonZeros.get(j));
		}
	}

	// add edges
	for(int i=0; i<signatureNonZeros.size(); ++i) {
		for(int j=0; j<functionNonZeros.size(); ++j) {
			if(candidateList.get(signatureNonZeros.get(i), functionNonZeros.get(j)) != 0) {
				g.addEdge("s"+signatureNonZeros.get(i), "f"+functionNonZeros.get(j));
			}
		}
	}

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

       // bipartite matching!
	HopcroftKarpBipartiteMatching<String, DefaultEdge> alg = 
		new HopcroftKarpBipartiteMatching<String, DefaultEdge>(g, p1, p2);

	Set<DefaultEdge> match = alg.getMatching();
	// sat || unsat
	if(match.size() == signatureNodeCount) {
		return true;
	} else { 
		return false;
	}

}
 

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