size(900,568,P3D);

// The inimmitable graph visualization explorer
public class Adventurer {
	public Graph g;
	public Adventurer( Graph base ) {
		g = base;
		make_all_basic_matrices();
	}

	public void make_all_basic_matrices() {
		int max = g.numbering.size();
		g.matrix = new Matrix[Graph.MatrixTypeIndex.MATRIX_FIELD_COUNT];
		float[] threeArgs = new float[] { (float)Math.random()*max, (float)Math.random()*max, (float)Math.random()*max };
		for(int i = 0; i < Graph.MatrixTypeIndex.MATRIX_FIELD_COUNT;i++) {
			Algebraicist.form(g, i, threeArgs);
		}
	}

	public Matrix make_random_filter(int rand_base_matrix1, int rand_base_matrix2, int rand_power, int rand_binary ) {
		//System.out.println("RF: " + rand_base_matrix1 + "," + rand_base_matrix2 + "," + rand_power + "," + rand_binary);
		Matrix f = new Matrix(g.matrix[rand_base_matrix1]);
		int max = g.numbering.size();
		float[] threeArgs = new float[] { (float)Math.random()*max, (float)Math.random()*max, (float)Math.random()*max };
		f = Arithmetician.fastMultiproduct(f, rand_power);
		f = Arithmetician.ApplyBinary(rand_binary, f, g.matrix[rand_base_matrix2], threeArgs);
		return f;
	}

	public Matrix make_random_Matrix_filter(int mod_or_not, int[] param_Base_Matrix, int[] param_Unary_Ops, int[] param_Binary_Combinators, int[] zeroOne, float[] threeArgs) {
		// We choose a filter length of max 5 matrices
		// We choose five matrices, and five unary filters to apply to each matrix
		// We choose 4 binary combinator ops to combine these 5 matrices
		// We collapse these ops right to left to construct a random filter
		// But it is only random in the sense that we choose our exploration of the filterspace randomly (or using a basic GA)
		// It is deterministic on the parameters passed in here
		// Which in GA are 'genetic code' / feature vector
		Matrix[] theFive = new Matrix[5];
		int i = 5;
		while(--i >= 0)
			theFive[i] = g.matrix[param_Base_Matrix[i]];
		if(mod_or_not == 0) {
			i = 5;
			while(--i >= 0)
				theFive[i] = Arithmetician.ApplyUnary(param_Unary_Ops[i], theFive[i], threeArgs, zeroOne );
			i = 4;
			Matrix result = theFive[4];
			while(--i >= 0)
				result = Arithmetician.ApplyBinary(param_Binary_Combinators[i], theFive[i], result, threeArgs);
			return result;
		}
		else {
			i = 5;
			while(--i >= 0)
				theFive[i] = Arithmetician.ApplyModUnary(param_Unary_Ops[i], theFive[i], threeArgs, zeroOne );
			i = 4;
			Matrix result = theFive[4];
			while(--i >= 0)
				result = Arithmetician.ApplyModBinary(param_Binary_Combinators[i], theFive[i], result, threeArgs);
			return result;
		}
	}
}
import java.util.HashMap;
import java.util.Iterator;
import java.util.Set;

// He makes stuff
public class Algebraicist {
	// 5 primitive matrices

	public static Matrix form( Graph g, int matrixType, float[] threeArgs ) {
		switch(matrixType) {
		case(Graph.MatrixTypeIndex.ADJACENCY):
			return g.matrix[Graph.MatrixTypeIndex.ADJACENCY] = adjacency(g);
		case(Graph.MatrixTypeIndex.DEGREE):
			return g.matrix[Graph.MatrixTypeIndex.DEGREE] = degree(g);
		case(Graph.MatrixTypeIndex.ADJACENCY_2STEP):
			return g.matrix[Graph.MatrixTypeIndex.ADJACENCY_2STEP] = exactlyTwoStepAdjacency(g.matrix[Graph.MatrixTypeIndex.ADJACENCY]);
		case(Graph.MatrixTypeIndex.DEGREE_1STEP):
			return g.matrix[Graph.MatrixTypeIndex.DEGREE_1STEP] = exactlyOneStepDegree(g.matrix[Graph.MatrixTypeIndex.DEGREE], g.matrix[Graph.MatrixTypeIndex.ADJACENCY]);
		case(Graph.MatrixTypeIndex.DEGREE_2STEP):
			return g.matrix[Graph.MatrixTypeIndex.DEGREE_2STEP] = exactlyTwoStepDegree(g.matrix[Graph.MatrixTypeIndex.DEGREE], g.matrix[Graph.MatrixTypeIndex.ADJACENCY_2STEP]);
		case(Graph.MatrixTypeIndex.SEIDELADJACENCY):
			return g.matrix[Graph.MatrixTypeIndex.SEIDELADJACENCY] = seidelAdjacency(g.matrix[Graph.MatrixTypeIndex.ADJACENCY]);
		case(Graph.MatrixTypeIndex.AUGMENTEDADJACENCY):
			return g.matrix[Graph.MatrixTypeIndex.AUGMENTEDADJACENCY] = augmentedAdjacencyMatrix(g.matrix[Graph.MatrixTypeIndex.ADJACENCY] );
		case(Graph.MatrixTypeIndex.SHAREOFDEGREE):
			return g.matrix[Graph.MatrixTypeIndex.SHAREOFDEGREE] =  shareOfDegree( g, g.matrix[Graph.MatrixTypeIndex.DEGREE] );
		case(Graph.MatrixTypeIndex.MINIMUMDISTANCEMATRIX):
			return g.matrix[Graph.MatrixTypeIndex.MINIMUMDISTANCEMATRIX] = minimumDistanceMatrix(g.matrix[Graph.MatrixTypeIndex.ADJACENCY]);
		case(Graph.MatrixTypeIndex.GENERALIZEDADJACENCY):
			return g.matrix[Graph.MatrixTypeIndex.GENERALIZEDADJACENCY] = generalizedAdjacency(  threeArgs[0],  threeArgs[1],  threeArgs[2], g.matrix[Graph.MatrixTypeIndex.ADJACENCY] );
		case(Graph.MatrixTypeIndex.WEIGHTEDDEGREE):
			return g.matrix[Graph.MatrixTypeIndex.WEIGHTEDDEGREE] = weightedDegree( g.matrix[Graph.MatrixTypeIndex.DEGREE], threeArgs[0] );
		case(Graph.MatrixTypeIndex.WEIGHTEDADJACENCY):
			return g.matrix[Graph.MatrixTypeIndex.WEIGHTEDADJACENCY] = weightedAdjacency( g.matrix[Graph.MatrixTypeIndex.ADJACENCY], threeArgs[0] );
		case(Graph.MatrixTypeIndex.SPIELMANLAZYWALK):
			return g.matrix[Graph.MatrixTypeIndex.SPIELMANLAZYWALK] = SpielmanLazywalkMatrix(g.matrix[Graph.MatrixTypeIndex.DEGREE], g.matrix[Graph.MatrixTypeIndex.ADJACENCY]);
		case(Graph.MatrixTypeIndex.LAZYWALKSOFTINVERSE):
			return g.matrix[Graph.MatrixTypeIndex.LAZYWALKSOFTINVERSE] = lazywalkMatrixSoftInverse(g.matrix[Graph.MatrixTypeIndex.DEGREE], g.matrix[Graph.MatrixTypeIndex.ADJACENCY]);
		case(Graph.MatrixTypeIndex.NSTEPWALKS):
			return g.matrix[Graph.MatrixTypeIndex.NSTEPWALKS] = Arithmetician.fastMultiproduct(g.matrix[Graph.MatrixTypeIndex.ADJACENCY], (int)threeArgs[0]);
		case(Graph.MatrixTypeIndex.TRANSITIONMATRIXRANDOMWALK):
			return g.matrix[Graph.MatrixTypeIndex.TRANSITIONMATRIXRANDOMWALK] = transitionMatrixRandomWalk( g.matrix[Graph.MatrixTypeIndex.ADJACENCY], g.matrix[Graph.MatrixTypeIndex.DEGREE] );
		case(Graph.MatrixTypeIndex.LAPLACIAN):
			return g.matrix[Graph.MatrixTypeIndex.LAPLACIAN] = laplacian(g.matrix[Graph.MatrixTypeIndex.DEGREE], g.matrix[Graph.MatrixTypeIndex.ADJACENCY], threeArgs[0]);
		case(Graph.MatrixTypeIndex.SOFTNORMALIZEDLAPLACIAN):
			return g.matrix[Graph.MatrixTypeIndex.SOFTNORMALIZEDLAPLACIAN] = SoftNormalizedLaplacian(g.matrix[Graph.MatrixTypeIndex.DEGREE], g.matrix[Graph.MatrixTypeIndex.LAPLACIAN], threeArgs[0]);
		case(Graph.MatrixTypeIndex.HARDNORMALIZEDLAPLACIAN):
			return g.matrix[Graph.MatrixTypeIndex.HARDNORMALIZEDLAPLACIAN] = HardNormalizedLaplacian(g.matrix[Graph.MatrixTypeIndex.DEGREE], g.matrix[Graph.MatrixTypeIndex.LAPLACIAN], threeArgs[0]);
		case(Graph.MatrixTypeIndex.COMMUTETIME):
			return g.matrix[Graph.MatrixTypeIndex.COMMUTETIME] = commuteTime(g);
		case(Graph.MatrixTypeIndex.DEFORMEDLAPLACIAN):
			return g.matrix[Graph.MatrixTypeIndex.DEFORMEDLAPLACIAN] = deformedLaplacian(g.matrix[Graph.MatrixTypeIndex.DEGREE], g.matrix[Graph.MatrixTypeIndex.ADJACENCY], threeArgs[0]);
		case(Graph.MatrixTypeIndex.PROJECTIONMATRIX):
			return g.matrix[Graph.MatrixTypeIndex.PROJECTIONMATRIX] = projectionMatrix((int)threeArgs[0], (int)threeArgs[1]);
		case(Graph.MatrixTypeIndex.MAXIMUMDISTANCEMATRIX):
			return g.matrix[Graph.MatrixTypeIndex.MAXIMUMDISTANCEMATRIX] = maximumDistanceMatrix(g.matrix[Graph.MatrixTypeIndex.COMMUTETIME]);
		default:
			return new Matrix(Matrix.IDENTITY,g.numbering.size(),g.numbering.size());
		}
	}

	public static Matrix adjacency( Graph g ) {
		int vertices = g.edgesByVertexA.size();
		Set<Graph.Vertex> verticeSet = g.edgesByVertexA.keySet();
		Iterator<Graph.Vertex> vertexIterator = verticeSet.iterator();
		Matrix am = new Matrix(Matrix.NONE, vertices, vertices);
		HashMap<Graph.Vertex,Integer> numbering = g.numbering();
		while(vertexIterator.hasNext()) {
			Graph.Vertex v = vertexIterator.next();
			int vnum = numbering.get(v);
			Set<Graph.Edge> vOutEdges = g.edgesByVertexA.get(v);
			Set<Graph.Edge> vInEdges = g.edgesByVertexB.get(v);
			Iterator<Graph.Edge> vOutIterator = vOutEdges.iterator();
			while(vOutIterator.hasNext()) {
				Graph.Edge e = vOutIterator.next();
				Graph.Vertex b = e.b;
				int bnum = numbering.get(b);
				am.entry[vnum][bnum] = 1.0f;
			}
			Iterator<Graph.Edge> vInIterator = vInEdges.iterator();
			while(vInIterator.hasNext()) {
				Graph.Edge e = vInIterator.next();
				Graph.Vertex a = e.a;
				int anum = numbering.get(a);
				am.entry[anum][vnum] = 1.0f;
			}
		}
		am.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return am;
	}

	public static Matrix incidence( Graph g ) { // VxE
		int vertices = g.numbering.size();
		int edges = g.edgeNumbering.size();
		Matrix im = new Matrix(Matrix.NONE, vertices, edges );
		Set<Graph.Edge> edgeSet = g.edgeNumbering.keySet();
		Iterator<Graph.Edge> edgeIterator = edgeSet.iterator();
		while(edgeIterator.hasNext()) {
			Graph.Edge e = edgeIterator.next();
			int edgnum = g.edgeNumbering.get(e);
			Graph.Vertex a = e.a;
			Graph.Vertex b = e.b;
			int anum = g.numbering.get(a);
			int bnum = g.numbering.get(b);
			im.entry[anum][edgnum] = 1.0f;
			im.entry[bnum][edgnum] = 1.0f;
		}
		return im;
	}

	public static Matrix SignedIncidence( Graph g ) { // VxE
		int vertices = g.numbering.size();
		int edges = g.edgeNumbering.size();
		Matrix im = new Matrix(Matrix.NONE, vertices, edges );
		Set<Graph.Edge> edgeSet = g.edgeNumbering.keySet();
		Iterator<Graph.Edge> edgeIterator = edgeSet.iterator();
		while(edgeIterator.hasNext()) {
			Graph.Edge e = edgeIterator.next();
			int edgnum = g.edgeNumbering.get(e);
			Graph.Vertex a = e.a;
			Graph.Vertex b = e.b;
			int anum = g.numbering.get(a);
			int bnum = g.numbering.get(b);
			im.entry[anum][edgnum] = 1.0f;
			im.entry[bnum][edgnum] = -1.0f;
		}
		return im;
	}

	public static Matrix degree( Graph g ) {
		int vertices = g.edgesByVertexA.size();
		Matrix dg = new Matrix(Matrix.NONE, vertices, vertices);
		HashMap<Graph.Vertex,Integer> numbering = g.numbering();
		Set<Graph.Vertex> verticeSet = g.edgesByVertexA.keySet();
		Iterator<Graph.Vertex> vertexIterator = verticeSet.iterator();
		g.volume = 0.0f;
		while(vertexIterator.hasNext()) {
			Graph.Vertex v = vertexIterator.next();
			int vnum = numbering.get(v);
			Set<Graph.Edge> vOutEdges = g.edgesByVertexA.get(v);
			Set<Graph.Edge> vInEdges = g.edgesByVertexB.get(v);
			g.volume += dg.entry[vnum][vnum] = vOutEdges.size() + vInEdges.size();
		}
		dg.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return dg;
	}

	public static Matrix edgeAdjacency( Graph g ) { // ExE
		Set<Graph.Vertex> verticeSet = g.edgesByVertexA.keySet();
		Iterator<Graph.Vertex> vertexIterator = verticeSet.iterator();
		int edges = g.edgeNumbering.size();
		Matrix eda = new Matrix(Matrix.NONE, edges, edges);
		while(vertexIterator.hasNext()) {
			Graph.Vertex stupidFuckFace = vertexIterator.next();
			Set<Graph.Edge> fuckingMoron = g.edgesByVertexA.get(stupidFuckFace);
			Iterator<Graph.Edge> edgeIterator = fuckingMoron.iterator();
			while(edgeIterator.hasNext()) {
				Graph.Edge e = edgeIterator.next();
				int edgnum = g.edgeNumbering.get(e);
				Graph.Vertex a = e.a;
				Set<Graph.Edge> heIsStillAFuckFace = g.edgesByVertexA.get(a);
				Iterator<Graph.Edge> fuckIterator = heIsStillAFuckFace.iterator();
				while(fuckIterator.hasNext()) {
					Graph.Edge eneighbour = fuckIterator.next();
					int eneighbournum = g.edgeNumbering.get(eneighbour);
					if(edgnum != eneighbournum)
						eda.entry[edgnum][eneighbournum] = 1.0f;
				}
				if(a!=e.b) {
					Graph.Vertex b = e.b;
					heIsStillAFuckFace = g.edgesByVertexA.get(b);
					fuckIterator = heIsStillAFuckFace.iterator();
					while(fuckIterator.hasNext()) {
						Graph.Edge eneighbour = fuckIterator.next();
						int eneighbournum = g.edgeNumbering.get(eneighbour);
						if(edgnum != eneighbournum)
							eda.entry[edgnum][eneighbournum] = 1.0f;
					}	
				}
			}
		}
		eda.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return eda;
	}

	public static Matrix edgeDegree( Graph g ) { // ExE
		Set<Graph.Vertex> verticeSet = g.edgesByVertexA.keySet();
		Iterator<Graph.Vertex> vertexIterator = verticeSet.iterator();
		int edges = g.edgeNumbering.size();
		Matrix edd = new Matrix(Matrix.NONE, edges, edges);
		while(vertexIterator.hasNext()) {
			Graph.Vertex stupidFuckFace = vertexIterator.next();
			Set<Graph.Edge> fuckingMoron = g.edgesByVertexA.get(stupidFuckFace);
			Iterator<Graph.Edge> edgeIterator = fuckingMoron.iterator();
			while(edgeIterator.hasNext()) {
				Graph.Edge e = edgeIterator.next();
				int edgnum = g.edgeNumbering.get(e);
				Graph.Vertex a = e.a;
				Set<Graph.Edge> heIsStillAFuckFace = g.edgesByVertexA.get(a);
				Iterator<Graph.Edge> fuckIterator = heIsStillAFuckFace.iterator();
				while(fuckIterator.hasNext()) {
					Graph.Edge eneighbour = fuckIterator.next();
					int eneighbournum = g.edgeNumbering.get(eneighbour);
					if(edgnum != eneighbournum)
						edd.entry[edgnum][edgnum] += 1.0f;
				}
				if(a!=e.b) {
					Graph.Vertex b = e.b;
					heIsStillAFuckFace = g.edgesByVertexA.get(b);
					fuckIterator = heIsStillAFuckFace.iterator();
					while(fuckIterator.hasNext()) {
						Graph.Edge eneighbour = fuckIterator.next();
						int eneighbournum = g.edgeNumbering.get(eneighbour);
						if(edgnum != eneighbournum)
							edd.entry[edgnum][edgnum] += 1.0f;
					}	
				}
			}
		}
		edd.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return edd;
	}


	// 5 secondary matrix operators

	public static Matrix complement( Matrix a ) {
		int i = a.rows;
		int j = a.columns;
		Matrix c = new Matrix(Matrix.NONE, i, j);
		while( --i >= 0 ) {
			while( --j >= 0 ) {
				if( a.entry[i][j] != 0.0f) 
					c.entry[i][j] = 0.0f;
				else
					c.entry[i][j] = 1.0f;
			}
			j = a.columns;
		}
		return c;
	}

	public static Matrix exactlyTwoStepAdjacency( Matrix a ) {
		int i = a.rows;
		int j = a.columns;
		int shortest = a.shortest;
		Matrix a2 = new Matrix(Matrix.NONE, i, j);
		while( --i >= 0 ) {
			while( --j >= 0 ) {
				int k = shortest;
				while ( --k >= 0 ) {
					if( a.entry[i][k] != 0.0f && a.entry[k][j] != 0.0f )
						a2.entry[i][j] = a.entry[i][k] + a.entry[k][j];
				}
			}
			j = a.columns;
		}
		return a2;
	}

	public static Matrix exactlyOneStepDegree( Matrix d, Matrix a ) {
		int i = a.rows;
		int j = a.columns;
		Matrix d1 = new Matrix(Matrix.NONE, i, j);
		while( --i >= 0 ) {
			while( --j >= 0 ) {
				if( a.entry[i][j] != 0.0f) 
					d1.entry[i][i] += d.entry[j][j];
			}
			j = a.columns;
		}
		return d1;
	}

	public static Matrix exactlyTwoStepDegree( Matrix d, Matrix a2 ) {
		int i = a2.rows;
		int j = a2.columns;
		Matrix d2 = new Matrix(Matrix.NONE, i, j);
		while( --i >= 0 ) {
			while( --j >= 0 ) {
				if( a2.entry[i][j] != 0.0f) 
					d2.entry[i][i] += d.entry[j][j];
			}
			j = a2.columns;
		}
		return d2;
	}

	public static Matrix seidelAdjacency( Matrix a ) {
		int i = a.rows;
		int j = a.columns;
		Matrix sa = new Matrix(Matrix.NONE, i, j);
		while( --i >= 0 ) {
			while( --j >= 0 ) {
				if(i==j)
					sa.entry[i][i] = 0.0f;
				else if(a.entry[i][j] != 0.0f)
					sa.entry[i][j] = -1.0f;
				else
					sa.entry[i][j] = 1.0f;
			}
			j = a.columns;
		}
		return sa;
	}

	// 22 Matrix combinators 

	public static Matrix augmentedAdjacencyMatrix( Matrix incidence ) {
		return Arithmetician.productResult(incidence, incidence.transpose());
	}

	public static Matrix FanChungBoundaryOperator( Matrix SignedIncidence, Matrix degree ) {
		int i = SignedIncidence.rows;
		int j = SignedIncidence.columns;
		Matrix sd = Arithmetician.entrywiseSquareRootResult(degree);
		Matrix b = new Matrix(Matrix.NONE, j, i);
		while( --i >= 0 ) {
			while( --j >= 0 ) {
				b.entry[j][i] = SignedIncidence.entry[i][j]/sd.entry[i][i];
			}
			j = SignedIncidence.columns;
		}
		return b;
	}

	public static Matrix shareOfDegree( Graph g, Matrix degree ) {
		int i = degree.rows;
		int j = degree.columns;
		Matrix sd = new Matrix(Matrix.NONE, i, j);
		while( --i >= 0 ) {
			while( --j >= 0 ) {
				sd.entry[i][j] = degree.entry[i][j]/g.volume;
			}
			j = degree.columns;
		}
		return sd;
	}

	public static Matrix minimumDistanceMatrix(Matrix connectivityMetric) {
		return Arithmetician.FloydWarshall(connectivityMetric, 1.0f, 50.0f, connectivityMetric.shortest*4.0f*50.0f, 0, Arithmetician.MIN );
	}

	public static Matrix generalizedAdjacency( float x, float y, float z, Matrix a ) {
		//(xI +yA + z(J-I-A))
		return Arithmetician.pairwiseSumResult(
				Arithmetician.scalarMultiplyResult( 
						new Matrix(Matrix.IDENTITY, a.rows, a.columns), x), 
						Arithmetician.pairwiseSumResult(
								Arithmetician.scalarMultiplyResult(a, y),
								Arithmetician.scalarMultiplyResult(
										Arithmetician.pairwiseSubtractResult(
												new Matrix(Matrix.COUNTERIDENTITY, a.rows, a.columns),
												Arithmetician.pairwiseSumResult(new Matrix(Matrix.IDENTITY, a.rows, a.columns), a)),
												z)
								));

	}

	public static Matrix weightedDegree( Matrix degree, float weight ) {
		return Arithmetician.scalarMultiplyResult(degree, weight);
	}

	public static Matrix weightedAdjacency( Matrix adjacency, float weight ) {
		return Arithmetician.scalarMultiplyResult(adjacency, weight);
	}

	public static Matrix SpielmanLazywalkMatrix(Matrix degree, Matrix adjacency) {
		return Arithmetician.scalarMultiplyResult(
				Arithmetician.pairwiseSumResult(
						new Matrix(Matrix.IDENTITY, adjacency.rows, adjacency.columns),
						Arithmetician.productResult(adjacency, 
								Arithmetician.Inverse(degree))), 
								0.5f);
	}

	public static Matrix lazywalkMatrixSoftInverse(Matrix degree, Matrix adjacency) {
		return Arithmetician.scalarMultiplyResult(
				Arithmetician.pairwiseSumResult(
						new Matrix(Matrix.IDENTITY, adjacency.rows, adjacency.columns),
						Arithmetician.productResult(degree, 
								Arithmetician.Inverse(adjacency))), 
								2.0f);
	}

	public static Matrix nStepWalks( Matrix adjacency, int nsteps ) {
		return Arithmetician.fastMultiproduct(adjacency, nsteps);
	}

	public static Matrix transitionMatrixRandomWalk( Matrix a, Matrix degree ) {
		int i = a.rows;
		int j = a.columns;
		Matrix t = new Matrix(Matrix.NONE, i, j);
		while( --i >= 0 ) {
			while( --j >= 0 ) {
				if( a.entry[i][j] != 0.0f )
					t.entry[i][j] = 1.0f/degree.entry[i][i];
			}
			j = a.columns;
		}
		return t;
	}

	public static Matrix laplacian(Matrix degree, Matrix adjacency, float weight) {
		if(weight != 1.0f)
			return Arithmetician.scalarMultiplyResult(
					Arithmetician.pairwiseSubtractResult(degree, adjacency),
					weight);
		else
			return 	Arithmetician.pairwiseSubtractResult(degree, adjacency);
	}

	public static Matrix laplaceOperator( Matrix d ) {
		int i = d.rows;
		int j = d.columns;
		Matrix lo = new Matrix(Matrix.NONE, i, j);
		while( --i >= 0 ) {
			while( --j >= 0 ) {
				if(i==j)
					lo.entry[i][i] = 1.0f;
				else
					lo.entry[i][j] = -1.0f/d.entry[i][j]; 
			}
			j = d.columns;
		}
		return lo;
	}

	public static Matrix SoftNormalizedLaplacian(Matrix degree, Matrix laplacian, float weight ) {
		Matrix entryWiseRoot = Arithmetician.entrywiseSquareRootResult(degree);
		Matrix entryWiseRootInverse = Arithmetician.entrywiseInversionResult(entryWiseRoot);
		if(weight != 1.0f)
			return Arithmetician.productResult(
					Arithmetician.productResult(entryWiseRootInverse, Arithmetician.scalarMultiplyResult(laplacian, weight)),
					entryWiseRoot);
		else
			return Arithmetician.productResult(
					Arithmetician.productResult(entryWiseRootInverse, laplacian),
					entryWiseRoot);
	}

	public static Matrix HardNormalizedLaplacian(Matrix degree, Matrix laplacian, float weight ) {
		Matrix[] rootAndInverse = Arithmetician.squareRootInverseSquareRootResult(degree);
		if(weight != 1.0f)
			return Arithmetician.productResult(
					Arithmetician.productResult(rootAndInverse[1], Arithmetician.scalarMultiplyResult(laplacian, weight)),
					rootAndInverse[0]);

		else
			return Arithmetician.productResult(
					Arithmetician.productResult(rootAndInverse[1], laplacian),
					rootAndInverse[0]);
	}

	public static Matrix eigenvectors( Matrix m ) { 
		Arithmetician.PowerIterationDeflate(m, m.shortest, 0.000000001, 0 );
		Matrix e = new Matrix(m.eigenvector);
		System.arraycopy(m.eigenvalue, 0, e.eigenvalue, 0, m.eigenvalue.length);
		return e;
	}

	public static void show_v (float[] v) {
		int i = v.length;
		while(--i >= 0)
			System.out.print(v[i] + " ");
		System.out.println();
	}


	public static Matrix commuteTime( Graph g ) {  
		Matrix nl = g.matrix[Graph.MatrixTypeIndex.SOFTNORMALIZEDLAPLACIAN];
		Matrix d = Arithmetician.entrywiseSquareRootResult(g.matrix[Graph.MatrixTypeIndex.DEGREE]);
		Matrix eigns = eigenvectors(nl);
		Matrix ct = new Matrix(Matrix.NONE, d.rows, d.columns);
		float vol = g.volume;
		int k = nl.shortest;
		while(--k >= 0 ) {
			int i = d.rows;
			while(--i >= 0) {
				int j = d.columns;
				while(--j >= 0) {
					float eki = eigns.entry[k][i]/d.entry[i][i];
					float ekj = eigns.entry[k][j]/d.entry[j][j];
					ct.entry[i][j] += vol*(eki-ekj)*(eki-ekj)/eigns.eigenvalue[k];
				}
			}
		}
		return ct;
	}

	public static Matrix deformedLaplacian(Matrix d, Matrix a, float deformation ) {
		//I - sA + s2(D - I)
		return Arithmetician.pairwiseSumResult(
				Arithmetician.pairwiseSubtractResult(
						new Matrix(Matrix.IDENTITY, a.rows, a.columns ),
						Arithmetician.scalarMultiplyResult(a, deformation)),
						Arithmetician.scalarMultiplyResult(
								Arithmetician.pairwiseSubtractResult(d, new Matrix(Matrix.IDENTITY, d.rows, d.columns)),
								deformation*deformation));

	}

	public static Matrix projectionMatrix( int rows, int columns ) {
		float ns = (float)Math.sqrt(rows*columns);
		float diag = (ns-1.0f)/ns;
		float nondiag = -1.0f/ns;
		int i = rows;
		int j = columns;
		Matrix pm = new Matrix(Matrix.NONE, i, j);
		while( --i >= 0 ) {
			while( --j >= 0 ) {
				if(i==j)
					pm.entry[i][i] = diag;
				else
					pm.entry[i][j] = nondiag; 
			}
			j = columns;
		}
		return pm;
	}

	public static Matrix projectedMatrix(Matrix embeddingMetric, float projectionCoefficient ) { // 0.5f
		// p = 0.5PEP
		return Arithmetician.productResult(
				Arithmetician.scalarMultiplyResult(projectionMatrix(embeddingMetric.rows, embeddingMetric.columns), projectionCoefficient),
				Arithmetician.productResult( embeddingMetric, projectionMatrix(embeddingMetric.rows, embeddingMetric.columns)));
	}

	public static Matrix RPaugmentedDegreeMatrix(Matrix d, Matrix mindistance) {
		Matrix mindistanceRoot = Arithmetician.entrywiseSquareRootResult(mindistance);
		int i = d.rows;
		int j = d.columns;
		Matrix augd = new Matrix(Matrix.NONE, i, j);
		while( --i >= 0 ) {
			while( --j >= 0 ) {
				if(i==j)
					augd.entry[i][i] = d.entry[i][i];
				else if ((1 << (int)(mindistanceRoot.entry[i][j])) != 0.0f )
					augd.entry[i][j] = d.entry[j][j] / (1 << (int)(mindistanceRoot.entry[i][j])) ;
			}
			j = d.columns;
		}
		return augd;
	}

	public static Matrix maximumDistanceMatrix(Matrix connectivityMetric) {
		return Arithmetician.FloydWarshall(connectivityMetric, 1.0f, 50.0f, connectivityMetric.shortest*4.0f*50.0f, 1, Arithmetician.MAX );
	}












}
import java.util.Arrays;

public final class Arithmetician {
	public static final int MIN = 0;
	public static final int MAX = 1;
	public int parametrization = (2*3*10*18*9*32);

	public class MatrixOperation {
		public static final int UNARY = 0;
		public static final int SCALAR = 1;
		public static final int BINARY = 2;
		public static final int MOD = 3;

		public class MatrixBinary {
			public static final int PAIRWISESUM = 0;
			public static final int PAIRWISESUBTRACT = 1;
			public static final int HADARMARDPRODUCT = 2;
			public static final int PRODUCTPROPAGATED = 3;
			public static final int PRODUCT = 4;
			public static final int HADAMARDDIVISION = 5;
			public static final int AND = 6;
			public static final int OR = 7;
			public static final int XOR = 8;
			public static final int KOR = 9;
			public static final int FROBENIUSINNERPRODUCT = 10;
			public static final int BINARY_FIELD_COUNT = 11;

		}
		public class MatrixUnary {
			public static final int ENTRYWISEINVERSION = 1;
			public static final int EXPONENTIALTRACE = 2;
			public static final int SQUAREROOTINVERSESQUAREROOT = 3;
			public static final int SOFTSQUAREROOT = 4;
			public static final int ENTRYWISENATURALLOG = 5;
			public static final int RANDOMFILL = 6;
			public static final int L0NORMALIZE = 7;
			public static final int L1NORMALIZE = 8;
			public static final int L2NORMALIZE = 9;
			public static final int LINFNORMALIZE = 10;
			public static final int FLOYDWARSHALL = 11;
			public static final int POWERITERATIONDEFLATE = 12;
			public static final int LUDOLITTLEDECOMPOSITION = 13;
			public static final int DETERMINANT = 14;
			public static final int INVERSE = 15;
			public static final int NOT = 16;
			public static final int COMPLEMENT = 17;
			public static final int UNARY_FIELD_COUNT = 18;
		}
		public class MatrixScalar {
			public static final int SCALARADD = 0;
			public static final int SCALARSUBTRACT = 1;
			public static final int SCALARMULTIPLY = 2;
			public static final int SCALARDIVIDE = 3;
			public static final int HADAMARDPOWER = 4;		
			public static final int FASTPROPAGATEDMULTIPRODUCT = 5;
			public static final int FASTMULTIPRODUCT = 6;
			public static final int ENTRYWISEMAPTORANGE = 7;
			public static final int THRESHOLD = 8;

		}

		public class Mod {
			public static final int UNARY = 0;
			public static final int SCALAR = 1;
			public static final int BINARY = 2;
			public class MatrixBinary {
				public static final int PAIRWISESUM = 0;
				public static final int PAIRWISESUBTRACT = 1;
				public static final int HADARMARDPRODUCT = 2;
				public static final int PRODUCTPROPAGATED = 3;
				public static final int PRODUCT = 4;
				public static final int HADAMARDDIVISION = 5;
				public static final int AND = 6;
				public static final int OR = 7;
				public static final int XOR = 8;
				public static final int KOR = 9;

			}
			public class MatrixUnary {
				public static final int FROBENIUSINNERPRODUCT = 0;
				public static final int ENTRYWISEINVERSION = 1;
				public static final int EXPONENTIALTRACE = 2;
				public static final int SQUAREROOTINVERSESQUAREROOT = 3;
				public static final int SOFTSQUAREROOT = 4;
				public static final int ENTRYWISENATURALLOG = 5;
				public static final int RANDOMFILL = 6;
				public static final int L0NORMALIZE = 7;
				public static final int L1NORMALIZE = 8;
				public static final int L2NORMALIZE = 9;
				public static final int LINFNORMALIZE = 10;
				public static final int FLOYDWARSHALL = 11;
				public static final int POWERITERATIONDEFLATE = 12;
				public static final int LUDOLITTLEDECOMPOSITION = 13;
				public static final int DETERMINANT = 14;
				public static final int INVERSE = 15;
				public static final int NOT = 16;
				public static final int COMPLEMENT = 17;
			}
			public class MatrixScalar {
				public static final int SCALARADD = 0;
				public static final int SCALARSUBTRACT = 1;
				public static final int SCALARMULTIPLY = 2;
				public static final int SCALARDIVIDE = 3;
				public static final int HADAMARDPOWER = 4;		
				public static final int FASTPROPAGATEDMULTIPRODUCT = 5;
				public static final int FASTMULTIPRODUCT = 6;
				public static final int ENTRYWISEMAPTORANGE = 7;
				public static final int THRESHOLD = 8;

			}
		}

	}

	public static void vectorScalarAdd( float[] v, float add ) {
		int i = v.length;
		while(--i >= 0)
			v[i]+=add;
	}

	public static Matrix ApplyScalar( int opNum, Matrix m, float[] threeArgs) {
		switch(opNum) {
		case(MatrixOperation.MatrixScalar.SCALARADD):
			return scalarAddResult(m, threeArgs[0]);
		case(MatrixOperation.MatrixScalar.SCALARSUBTRACT):
			return scalarSubtractResult(m, threeArgs[0]);
		case(MatrixOperation.MatrixScalar.SCALARMULTIPLY):
			return scalarMultiplyResult(m, threeArgs[0]);
		case(MatrixOperation.MatrixScalar.SCALARDIVIDE):
			return scalarDivideResult(m, threeArgs[0]);
		case(MatrixOperation.MatrixScalar.HADAMARDPOWER):		
			return hadamardPowerResult(m, (int)threeArgs[0]);
		case(MatrixOperation.MatrixScalar.FASTPROPAGATEDMULTIPRODUCT):
			fastPropagatedMultiproduct(m, (int)threeArgs[0]);
		return m;
		case(MatrixOperation.MatrixScalar.FASTMULTIPRODUCT):
			return fastMultiproduct(m, (int)threeArgs[0]);
		case(MatrixOperation.MatrixScalar.ENTRYWISEMAPTORANGE):
			return scaleToRangeResult(m, threeArgs[0], threeArgs[1]);
		case(MatrixOperation.MatrixScalar.THRESHOLD):
			return thresholdResult(m, threeArgs[0], (int)threeArgs[1]);
		default:
			return new Matrix(m);
		}
	}

	public static Matrix ApplyModScalar( int opNum, Matrix m, float[] threeArgs) {
		switch(opNum) {
		case(MatrixOperation.MatrixScalar.SCALARADD):
			return ModscalarAddResult(m, threeArgs[0], threeArgs[2]);
		case(MatrixOperation.MatrixScalar.SCALARSUBTRACT):
			return ModscalarSubtractResult(m, threeArgs[0], threeArgs[2]);
		case(MatrixOperation.MatrixScalar.SCALARMULTIPLY):
			return ModscalarMultiplyResult(m, threeArgs[0], threeArgs[2]);
		case(MatrixOperation.MatrixScalar.SCALARDIVIDE):
			return ModscalarDivideResult(m, threeArgs[0], threeArgs[2]);
		case(MatrixOperation.MatrixScalar.HADAMARDPOWER):		
			return ModhadamardPowerResult(m, (int)threeArgs[0], threeArgs[2]);
		case(MatrixOperation.MatrixScalar.FASTPROPAGATEDMULTIPRODUCT):
			ModfastPropagatedMultiproduct(m, (int)threeArgs[0], threeArgs[2]);
		return m;
		case(MatrixOperation.MatrixScalar.FASTMULTIPRODUCT):
			return ModfastMultiproduct(m, (int)threeArgs[0], threeArgs[2]);
		case(MatrixOperation.MatrixScalar.ENTRYWISEMAPTORANGE):
			return ModscaleToRangeResult(m, threeArgs[0], threeArgs[1], threeArgs[2]);
		case(MatrixOperation.MatrixScalar.THRESHOLD):
			return ModthresholdResult(m, threeArgs[0], (int)threeArgs[1], threeArgs[2]);
		default:
			return new Matrix(m);
		}
	}


	public  static Matrix ApplyUnary( int opNum, Matrix m, float[] threeArgs, int[] zeroOne) {
		switch(opNum) {
		case(MatrixOperation.MatrixUnary.ENTRYWISEINVERSION):
			return entrywiseInversionResult(m);
		case(MatrixOperation.MatrixUnary.EXPONENTIALTRACE):
			float result = exponentialTrace(m);
		return scalarMultiplyResult(m,result);
		case(MatrixOperation.MatrixUnary.SQUAREROOTINVERSESQUAREROOT):
			return squareRootInverseSquareRootResult(m)[0];
		case(MatrixOperation.MatrixUnary.SOFTSQUAREROOT):
			return entrywiseSquareRootResult(m);
		case(MatrixOperation.MatrixUnary.ENTRYWISENATURALLOG):
			return entrywiseLogarithmResult(m);
		case(MatrixOperation.MatrixUnary.RANDOMFILL):
			fill_random(m);
		return m;
		case(MatrixOperation.MatrixUnary.L0NORMALIZE):
			L0normalize(m);
		return m;
		case(MatrixOperation.MatrixUnary.L1NORMALIZE):
			L1normalize(m);
		return m;
		case(MatrixOperation.MatrixUnary.L2NORMALIZE):
			L2normalize(m);
		return m;
		case(MatrixOperation.MatrixUnary.LINFNORMALIZE):
			Linfnormalize(m);
		return m;
		case(MatrixOperation.MatrixUnary.FLOYDWARSHALL):
			return FloydWarshall(m, 0.00000001f, threeArgs[0], threeArgs[1], zeroOne[0], zeroOne[1]);
		case(MatrixOperation.MatrixUnary.POWERITERATIONDEFLATE):
			m.eigenvalue = new float[m.rows];
		m.eigenvector = new float[m.rows][m.columns];
		PowerIterationDeflate(m, m.shortest, (double)0.00000001f, 1); // 1 populates krylov subspace
		return new Matrix(m.eigenvector);
		case(MatrixOperation.MatrixUnary.LUDOLITTLEDECOMPOSITION):
			return LU_Dolittle_DecompositionResult(m);
		case(MatrixOperation.MatrixUnary.DETERMINANT):
			result = Determinant(m);
		return scalarMultiplyResult(m,result);
		case(MatrixOperation.MatrixUnary.INVERSE):
			return Inverse(m);
		case(MatrixOperation.MatrixUnary.NOT):
			return NOTResult(m);
		case(MatrixOperation.MatrixUnary.COMPLEMENT):
			return complement(m);
		default:
			return new Matrix(m);
		}
	}

	public static Matrix ApplyModUnary( int opNum, Matrix m, float[] threeArgs, int[] zeroOne) {
		switch(opNum) {
		case(MatrixOperation.MatrixUnary.ENTRYWISEINVERSION):
			return ModentrywiseInversionResult(m, threeArgs[2]);
		case(MatrixOperation.MatrixUnary.EXPONENTIALTRACE):
			float result = ModexponentialTrace(m, threeArgs[2]);
		return ModscalarMultiplyResult(m, result, threeArgs[2]);
		case(MatrixOperation.MatrixUnary.SQUAREROOTINVERSESQUAREROOT):
			return ModsquareRootInverseSquareRootResult(m, threeArgs[2])[0];
		case(MatrixOperation.MatrixUnary.SOFTSQUAREROOT):
			return ModentrywiseSquareRootResult(m, threeArgs[2]);
		case(MatrixOperation.MatrixUnary.ENTRYWISENATURALLOG):
			return ModentrywiseLogarithmResult(m, threeArgs[2]);
		case(MatrixOperation.MatrixUnary.RANDOMFILL):
			Modfill_random(m, threeArgs[2]);
		return m;
		case(MatrixOperation.MatrixUnary.L0NORMALIZE):
			ModL0normalize(m, threeArgs[2]);
		return m;
		case(MatrixOperation.MatrixUnary.L1NORMALIZE):
			ModL1normalize(m, threeArgs[2]);
		return m;
		case(MatrixOperation.MatrixUnary.L2NORMALIZE):
			ModL2normalize(m, threeArgs[2]);
		return m;
		case(MatrixOperation.MatrixUnary.LINFNORMALIZE):
			ModLinfnormalize(m, threeArgs[2]);
		return m;
		case(MatrixOperation.MatrixUnary.FLOYDWARSHALL):
			return ModFloydWarshall(m, 0.00000001f, threeArgs[0], threeArgs[1], zeroOne[0], zeroOne[1], threeArgs[2]);
		case(MatrixOperation.MatrixUnary.POWERITERATIONDEFLATE):
			ModPowerIterationDeflate(m, m.shortest, (double)0.00000001f, 1, threeArgs[2]); // 1 populates krylov subspace
		return new Matrix(m.eigenvector);
		case(MatrixOperation.MatrixUnary.LUDOLITTLEDECOMPOSITION):
			return ModLU_Dolittle_DecompositionResult(m, threeArgs[2]);
		case(MatrixOperation.MatrixUnary.DETERMINANT):
			result = ModDeterminant(m, threeArgs[2]);
		return ModscalarMultiplyResult(m,result, threeArgs[2]);
		case(MatrixOperation.MatrixUnary.INVERSE):
			return ModInverse(m, threeArgs[2]);
		case(MatrixOperation.MatrixUnary.NOT):
			return ModNOTResult(m, threeArgs[2]);
		case(MatrixOperation.MatrixUnary.COMPLEMENT):
			return complement(m);
		default:
			return new Matrix(m);
		}
	}


	public  static Matrix ApplyBinary( int opNum, Matrix m, Matrix m2, float[] threeArgs) {
		switch(opNum) {
		case(MatrixOperation.MatrixBinary.FROBENIUSINNERPRODUCT):
			float result = FrobeniusInnerProduct(m, m2);
		return scalarMultiplyResult(m,result);
		case(MatrixOperation.MatrixBinary.PAIRWISESUM):
			return pairwiseSumResult(m, m2);
		case(MatrixOperation.MatrixBinary.PAIRWISESUBTRACT):
			return pairwiseSubtractResult(m, m2);
		case(MatrixOperation.MatrixBinary.HADARMARDPRODUCT):
			return hadamardProductResult(m, m2);
		case(MatrixOperation.MatrixBinary.PRODUCTPROPAGATED):
			productPropagated(m, m2);
		return m;
		case(MatrixOperation.MatrixBinary.PRODUCT):
			return productResult(m, m2);
		case(MatrixOperation.MatrixBinary.HADAMARDDIVISION):
			return hadamardDivisionResult(m, m2);
		case(MatrixOperation.MatrixBinary.AND):
			return ANDResult(m, m2);
		case(MatrixOperation.MatrixBinary.OR):
			return ORResult(m, m2);
		case(MatrixOperation.MatrixBinary.XOR):
			return XORResult(m, m2);
		case(MatrixOperation.MatrixBinary.KOR):
			return KORResult(m, m2);
		default:
			return new Matrix(m);
		}
	}
	//Addition


	public  static Matrix ApplyModBinary( int opNum, Matrix m, Matrix m2, float[] threeArgs) {
		switch(opNum) {
		case(MatrixOperation.MatrixBinary.FROBENIUSINNERPRODUCT):
			float result = ModFrobeniusInnerProduct(m, m2, threeArgs[2]);
		return ModscalarMultiplyResult(m,result, threeArgs[2]);
		case(MatrixOperation.MatrixBinary.PAIRWISESUM):
			return ModpairwiseSumResult(m, m2, threeArgs[2]);
		case(MatrixOperation.MatrixBinary.PAIRWISESUBTRACT):
			return ModpairwiseSubtractResult(m, m2, threeArgs[2]);
		case(MatrixOperation.MatrixBinary.HADARMARDPRODUCT):
			return ModhadamardProductResult(m, m2, threeArgs[2]);
		case(MatrixOperation.MatrixBinary.PRODUCTPROPAGATED):
			ModproductPropagated(m, m2, threeArgs[2]);
		return m;
		case(MatrixOperation.MatrixBinary.PRODUCT):
			return ModproductResult(m, m2, threeArgs[2]);
		case(MatrixOperation.MatrixBinary.HADAMARDDIVISION):
			return ModhadamardDivisionResult(m, m2, threeArgs[2]);
		case(MatrixOperation.MatrixBinary.AND):
			return ModANDResult(m, m2, threeArgs[2]);
		case(MatrixOperation.MatrixBinary.OR):
			return ModORResult(m, m2, threeArgs[2]);
		case(MatrixOperation.MatrixBinary.XOR):
			return ModXORResult(m, m2, threeArgs[2]);
		case(MatrixOperation.MatrixBinary.KOR):
			return ModKORResult(m, m2, threeArgs[2]);
		default:
			return new Matrix(m);
		}
	}

	public static void scalarAdd(Matrix accumulator, float addend) {
		if(addend == Matrix.ZEROF) {
			accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int shortestColumn = accumulator.columns;
		int i = accumulator.rows;
		int j = shortestColumn;
		while(--i >= 0) {
			while(--j >= 0)
				accumulator.entry[i][j] += addend;
		}
		accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static void scalarAdd(Matrix accumulator, float addend, int row) {
		if(addend == Matrix.ZEROF) {
			accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int shortestColumn = accumulator.columns;
		int i = row;
		int j = shortestColumn;
		while(--j >= 0)
			accumulator.entry[i][j] += addend;
		accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix scalarAddResult(Matrix adder, float addend) {
		Matrix result;
		if(addend == Matrix.ZEROF) {
			result = new Matrix(adder);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		int shortestRow = adder.rows;
		int shortestColumn = adder.columns;
		result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				result.entry[i][j] = adder.entry[i][j] + addend;
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static void pairwiseSum(Matrix accumulator, Matrix summand) {
		int shortestRow = accumulator.rows < summand.rows ? accumulator.rows : summand.rows;
		int shortestColumn = accumulator.columns < summand.columns ? accumulator.columns : summand.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				accumulator.entry[i][j] += summand.entry[i][j];
			}  
		}
		accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix pairwiseSumResult(Matrix adder, Matrix summand) {
		int shortestRow = adder.rows < summand.rows ? adder.rows : summand.rows;
		int shortestColumn = adder.columns < summand.columns ? adder.columns : summand.columns;
		Matrix result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				result.entry[i][j] = adder.entry[i][j] + summand.entry[i][j];
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	// Subtraction

	public static void scalarSubtract(Matrix accumulator, float subtractand) {
		if(subtractand == Matrix.ZEROF) {
			accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int shortestRow = accumulator.rows;
		int shortestColumn = accumulator.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				accumulator.entry[i][j] -= subtractand;
			}  
		}
		accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static void scalarSubtract(Matrix accumulator, float subtractand, int row) {
		if(subtractand == Matrix.ZEROF) {
			accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int shortestColumn = accumulator.columns;
		int i = row;
		int j = shortestColumn;
		while(--j >= 0)
			accumulator.entry[i][j] -= subtractand;
		accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix scalarSubtractResult(Matrix subtractor, float subtractand) {
		Matrix result;
		if(subtractand == Matrix.ZEROF) {
			result = new Matrix(subtractor);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		int shortestRow = subtractor.rows;
		int shortestColumn = subtractor.columns;
		result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				result.entry[i][j] = subtractor.entry[i][j] - subtractand;
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static void pairwiseSubtract(Matrix accumulator, Matrix subtractand) {
		int shortestRow = accumulator.rows < subtractand.rows ? accumulator.rows : subtractand.rows;
		int shortestColumn = accumulator.columns < subtractand.columns ? accumulator.columns : subtractand.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				accumulator.entry[i][j] -= subtractand.entry[i][j];
			}  
		}
		accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix pairwiseSubtractResult(Matrix subtractor, Matrix subtractand) {
		int shortestRow = subtractor.rows < subtractand.rows ? subtractor.rows : subtractand.rows;
		int shortestColumn = subtractor.columns < subtractand.columns ? subtractor.columns : subtractand.columns;
		Matrix result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				result.entry[i][j] = subtractor.entry[i][j] - subtractand.entry[i][j];
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	// Multiplication

	public static void scalarMultiply(Matrix accumulator, float multiplier) {
		int shortestRow = accumulator.rows;
		int shortestColumn = accumulator.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				accumulator.entry[i][j] *= multiplier;
			}  
		}
		accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static void scalarMultiply(Matrix accumulator, float multiplier, int row) {
		int shortestColumn = accumulator.columns;
		int i = row;
		int j = shortestColumn;
		while(--j >= 0)
			accumulator.entry[i][j] *= multiplier;
		accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix scalarMultiplyResult(Matrix multiplicand, float multiplier) {
		int shortestRow = multiplicand.rows;
		int shortestColumn = multiplicand.columns;
		Matrix result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				result.entry[i][j] = multiplicand.entry[i][j] * multiplier;
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static float FrobeniusInnerProduct( Matrix multiplicand, Matrix multiplier ) {
		int shortestRow = multiplicand.rows < multiplier.rows ? multiplicand.rows : multiplier.rows;
		int shortestColumn = multiplicand.columns < multiplier.columns ? multiplicand.columns : multiplier.columns;
		int i = shortestRow;
		float accumulator = 0.0f;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				accumulator += multiplicand.entry[i][j]*multiplier.entry[i][j];
			}  
		}
		return accumulator;
	}

	public static void hadamardProduct(Matrix accumulator, Matrix multiplicand) {
		int shortestRow = accumulator.rows < multiplicand.rows ? accumulator.rows : multiplicand.rows;
		int shortestColumn = accumulator.columns < multiplicand.columns ? accumulator.columns : multiplicand.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				accumulator.entry[i][j] *= multiplicand.entry[i][j];
			}  
		}
		accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix hadamardProductResult(Matrix multiplier, Matrix multiplicand) {
		int shortestRow = multiplier.rows < multiplicand.rows ? multiplier.rows : multiplicand.rows;
		int shortestColumn = multiplier.columns < multiplicand.columns ? multiplier.columns : multiplicand.columns;
		Matrix result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				result.entry[i][j] = multiplier.entry[i][j] * multiplicand.entry[i][j];
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static void productPropagated(Matrix accumulator, Matrix multiplier) {
		int shortestRow = accumulator.rows < multiplier.rows ? accumulator.rows : multiplier.rows;
		int shortestColumn = accumulator.columns < multiplier.columns ? accumulator.columns : multiplier.columns;
		int shortest = shortestRow < shortestColumn ? shortestRow : shortestColumn;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				int k = shortest;
				float entrySum = 0.0f;
				while(--k >= 0) {
					entrySum += accumulator.entry[i][k]*multiplier.entry[k][j];
				}
				accumulator.entry[i][j] = entrySum;
			}  
		}
		accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix productResult(Matrix multiplicand, Matrix multiplier) { 
		int shortest = multiplicand.columns < multiplier.rows ? multiplicand.columns : multiplier.rows;
		int i = multiplicand.rows;
		int j = multiplier.columns;
		Matrix result = new Matrix(Matrix.NONE, i, j);
		while(--i >= 0) {
			j = multiplier.columns;
			while(--j >= 0) {
				int k = shortest;
				float entrySum = 0.0f;
				while(--k >= 0) {
					entrySum += multiplicand.entry[i][k]*multiplier.entry[k][j];
				}
				result.entry[i][j] = entrySum;
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	// Division

	public static void entrywiseInversion( Matrix divisor ) {
		int shortestRow = divisor.rows;
		int shortestColumn = divisor.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float divisorij = divisor.entry[i][j];
				if(divisorij != Matrix.ZEROF)
					divisor.entry[i][j] = Matrix.IDENTITYF / divisorij;
				else
					divisor.entry[i][j] = Matrix.IDENTITYF;
			}  
		}
	}

	public static Matrix entrywiseInversionResult ( Matrix divisor ) {
		int shortestRow = divisor.rows;
		int shortestColumn = divisor.columns;
		Matrix result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float divisorij = divisor.entry[i][j];
				if(divisorij != Matrix.ZEROF)
					result.entry[i][j] = Matrix.IDENTITYF / divisorij;
				else
					result.entry[i][j] = Matrix.IDENTITYF;
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static void scalarDivide(Matrix accumulator, float divisor) {
		int shortestRow = accumulator.rows;
		int shortestColumn = accumulator.columns;
		if(divisor == Matrix.ZEROF) {
			Arrays.fill(accumulator.entry, Matrix.IDENTITYF);
			accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) 
				accumulator.entry[i][j] /= divisor;
		}
	}

	public static Matrix scalarDivideResult(Matrix dividend, float divisor) {
		int shortestRow = dividend.rows;
		int shortestColumn = dividend.columns;
		Matrix result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				if(divisor != Matrix.ZEROF)
					result.entry[i][j] = dividend.entry[i][j] / divisor;
				else
					result.entry[i][j] = Matrix.IDENTITYF;
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static void hadamardDivision(Matrix accumulator, Matrix divisor) {
		int shortestRow = accumulator.rows < divisor.rows ? accumulator.rows : divisor.rows;
		int shortestColumn = accumulator.columns < divisor.columns ? accumulator.columns : divisor.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				if(divisor.entry[i][j] != Matrix.ZEROF)
					accumulator.entry[i][j] /= divisor.entry[i][j];
				else
					accumulator.entry[i][j] = Matrix.ZEROF;
			}  
		}
		accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix hadamardDivisionResult(Matrix dividend, Matrix divisor) {
		int shortestRow = dividend.rows < divisor.rows ? dividend.rows : divisor.rows;
		int shortestColumn = dividend.columns < divisor.columns ? dividend.columns : divisor.columns;
		Matrix result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				if(divisor.entry[i][j] != Matrix.ZEROF)
					result.entry[i][j] = dividend.entry[i][j] / divisor.entry[i][j];
				else
					result.entry[i][j] = Matrix.IDENTITYF;
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	// Powers

	public static void hadamardPower(Matrix accumulator, int power) {
		int shortestRow = accumulator.rows;
		int shortestColumn = accumulator.columns;
		if(power == 0) {
			accumulator.make_ones();
			return;  
		}
		if(power == 1)
			return;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				int p = power;
				while(--p > 0)
					accumulator.entry[i][j] *= accumulator.entry[i][j];
			}  
		}
		accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix hadamardPowerResult(Matrix multiplier, int power) {

		int shortestRow = multiplier.rows;
		int shortestColumn = multiplier.columns;
		if(power == 1)
			return new Matrix(multiplier);
		Matrix result = new Matrix(Matrix.ONES, shortestRow, shortestColumn);
		if(power == 0)
			return result;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				int p = power;
				while(--p >= 0)
					result.entry[i][j] *= multiplier.entry[i][j];
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static int[] radix(int num, int base) {
		int b = 1; int b1 = base;
		int rem = 0; int rem1 = num % b1;
		int bit = 0;
		int[] bits = new int[10];
		while(b <= num) {
			bits[bit] = (rem1 - rem)/b;
			b=b1;b1*=base;
			rem=rem1; rem1 = num % b1;
			bit++;
			if(bit >= bits.length) {
				int[] temp = new int[bit<<1];
				System.arraycopy(temp, 0, bits, 0, bit);
				bits = temp;
			}
		}
		return bits;
	}

	public static float exponentialTrace( Matrix m ) {
		float result = 1.0f;
		float trace = m.trace();
		int itrace = (int)trace;
		float diff = trace - itrace;
		int[] bits = radix(itrace, 2);
		int iteration = 0;
		float lastPow = 2.71828183f;
		int lastPowIndex = 0;
		while(++iteration < bits.length) {
			int bit = bits[iteration];
			if(bit == 1) {
				int squaresRequired = iteration - lastPowIndex;
				while(--squaresRequired >= 0)
					lastPow *= lastPow;
				result *= lastPow;
			}
		}
		float extra = -(float)Math.log(diff);
		if(extra != 0.0f)
			result *= extra;
		return result;
	}

	public static void polysquarePropagated(Matrix base, int times) {
		while(--times >= 0)
			productPropagated(base, base);
	}

	public static void fastPropagatedMultiproduct(Matrix base, int exponent) { 
		int[] operations = radix(exponent, 2);
		int lastProductIndex = 0, jmax = operations.length;
		Matrix partialPowers;
		if(operations[0] == 0)
			partialPowers = new Matrix(Matrix.IDENTITY, base.rows, base.columns);
		else
			partialPowers = new Matrix(base);
		for(int j = 1; j < jmax;j++) {
			int bit = operations[j];
			if(bit == 1) {
				polysquarePropagated(partialPowers, j-lastProductIndex);
				lastProductIndex = j;
				productPropagated(base, partialPowers);
			}
		}
	}

	public static Matrix polysquare(Matrix base, int times) {
		while(--times >= 0)
			base = productResult(base, base);
		return base;
	}

	public static Matrix fastMultiproduct(Matrix base, int exponent) { 
		int[] operations = radix(exponent, 2);
		int lastProductIndex = 0, jmax = operations.length;
		Matrix partialPowers;
		if(operations[0] == 0)
			partialPowers = new Matrix(Matrix.IDENTITY, base.rows, base.columns);
		else
			partialPowers = new Matrix(base);
		for(int j = 1; j < jmax;j++) {
			int bit = operations[j];
			if(bit == 1) {
				polysquare(partialPowers, j-lastProductIndex);
				lastProductIndex = j;
				base = productResult(base, partialPowers);
			}
		}
		return base;
	}

	public static Matrix squareRootReturnInverseSquareRoot(Matrix y) {
		Matrix z = new Matrix(Matrix.IDENTITY, y.rows, y.columns);
		int j = (int)Math.sqrt(y.rows) + 1;
		while(--j >= 0) {
			scalarMultiply(pairwiseSumResult(y,Inverse(z)), 0.5f);
			scalarMultiply(pairwiseSumResult(z,Inverse(y)), 0.5f);
		}
		y.updates &= (Matrix.all_changed^Matrix.entry_changed);
		z.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return z;
	}

	public static Matrix[] squareRootInverseSquareRootResult(Matrix x) {
		Matrix y = new Matrix(x);
		Matrix z = new Matrix(Matrix.IDENTITY, y.rows, y.columns);
		int j = (int)Math.sqrt(y.rows) + 1;
		while(--j >= 0) {
			scalarMultiply(pairwiseSumResult(y,Inverse(z)), 0.5f);
			scalarMultiply(pairwiseSumResult(z,Inverse(y)), 0.5f);
		}
		y.updates &= (Matrix.all_changed^Matrix.entry_changed);
		z.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return new Matrix[] { y, z };
	}

	public static void entrywiseSquareRoot(Matrix x) {
		int shortestRow = x.rows;
		int shortestColumn = x.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float xij = x.entry[i][j];
				float sign;
				if(xij != 0.0f)
					sign = xij/Math.abs(xij);
				else 
					sign = 1.0f;
				x.entry[i][j] = (float)( sign*Math.sqrt(Math.abs(xij)) );
			}  
		}
		x.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix entrywiseSquareRootResult(Matrix x) {
		int shortestRow = x.rows;
		int shortestColumn = x.columns;
		int i = shortestRow;
		Matrix result = new Matrix(Matrix.NONE, x.rows, x.columns);
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float xij = x.entry[i][j];
				float sign;
				if(xij != 0.0f)
					sign = xij/Math.abs(xij);
				else 
					sign = 1.0f;
				result.entry[i][j] = (float)( sign*Math.sqrt(Math.abs(xij)) );
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}	

	// Logarithm

	public static void entrywiseLogarithm( Matrix m ) {
		int shortestRow = m.rows;
		int shortestColumn = m.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float eij = m.entry[i][j];
				float sign;
				if(eij != 0.0f)
					sign = eij/Math.abs(eij);
				else 
					sign = 1.0f;
				m.entry[i][j] = (float)(sign*Math.log(eij));
			}  
		}
		m.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix entrywiseLogarithmResult( Matrix m ) {
		int shortestRow = m.rows;
		int shortestColumn = m.columns;
		int i = shortestRow;
		Matrix result = new Matrix(m);
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float xij = m.entry[i][j];
				float sign;
				if(xij != 0.0f)
					sign = xij/Math.abs(xij);
				else 
					sign = 1.0f;
				result.entry[i][j] = (float)(sign*Math.log(xij));
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	// Mapping

	public static void scaleToRange(Matrix m, float low, float high) {
		float[] minmax = m.minmax();
		float newrange = high-low;
		float oldrange = minmax[1]-minmax[0];
		float range = oldrange == 0.0f ? 1.0f : newrange/oldrange;
		scalarSubtract(m, minmax[0]);
		scalarMultiply(m, range );
		scalarAdd(m, low );
	}

	public static void scaleToRangeRows(Matrix m, float low, float high) {
		int i = m.rows;
		while(--i >= 0) {
			float[] minmax = m.minmaxRow(i);
			float newrange = high-low;
			float oldrange = minmax[1]-minmax[0];
			float range = oldrange == 0.0f ? 1.0f : newrange/oldrange;
			scalarSubtract(m, minmax[0], i);
			scalarMultiply(m, range, i );
			scalarAdd(m, low, i );
		}
	}

	public static void scaleToRangeRow(Matrix m, int row, float low, float high) {
		int i = row;
		float[] minmax = m.minmaxRow(i);
		float newrange = high-low;
		float oldrange = minmax[1]-minmax[0];
		float range = oldrange == 0.0f ? 1.0f : newrange/oldrange;
		scalarSubtract(m, minmax[0], i);
		scalarMultiply(m, range, i );
		scalarAdd(m, low, i );

	}

	public static Matrix scaleToRangeResult(Matrix m, float low, float high) {
		float[] minmax = m.minmax();
		float newrange = high-low;
		float oldrange = minmax[1]-minmax[0];
		float range = oldrange == 0.0f ? 1.0f : newrange/oldrange;
		Matrix result = new Matrix(Matrix.NONE, m.rows, m.columns);
		result = scalarSubtractResult(m, minmax[0]);
		scalarMultiply(result, range );
		scalarAdd(result, low );
		return result;
	}


	// Random

	public static void fill_random(Matrix m) {
		int i = m.rows;
		while(--i >= 0 ) {
			int j = m.columns;
			while(--j >= 0 )
				m.entry[i][j] = (float)Math.random();
		}
	}

	public static void fill_random( float[] v) {
		int j = v.length;
		while(--j >= 0 ) v[j] = (float)Math.random();
	}


	// Norms & normalization vectors and matrices : 0, 1, 2, inf

	public static void Linfnormalize( Matrix m ) {
		int i = m.rows;
		int j = m.columns;
		if( i == 0 || j == 0)
			return;
		float max = Math.abs(m.entry[0][0]);
		while( --i >= 0 ) {
			while( --j >= 0 ) {
				float eij = Math.abs(m.entry[i][j]);
				if(eij > max)
					max = eij;
			}
			j = m.columns;
		}
		i = m.rows;
		if(max != 0.0f) {
			while( --i >= 0 ) {
				j = m.columns;
				while( --j >= 0 ) {
					m.entry[i][j] /= max;
				}
			}
		}
	}

	public static void Linfnormalize( float[] v ) {
		int j = v.length;
		if(j == 0)
			return;
		float max = Math.abs(v[0]);
		while(--j > 0) {
			float vj = Math.abs(v[j]);
			if(vj > max)
				max = vj;
		}
		j = v.length;
		if(max == 0.0f)
			return;
		while(--j >= 0) v[j] /= max;
	}

	public static float Linfnorm( Matrix m ) {
		int i = m.rows;
		int j = m.columns;
		if( i == 0 || j == 0)
			return 0.0f;
		float max = Math.abs(m.entry[0][0]);
		while( --i >= 0 ) {
			while( --j >= 0 ) {
				float eij = Math.abs(m.entry[i][j]);
				if(eij > max)
					max = eij;
			}
			j = m.columns;
		}
		if(m.norms == null)
			m.norms = new float[4];
		m.norms[3] = max;
		return max;
	}

	public static float Linfnorm( float[] v ) {
		int j = v.length;
		if(j == 0)
			return 0.0f;
		float max = Math.abs(v[0]);
		while(--j > 0) {
			float vj = Math.abs(v[j]);
			if(vj > max)
				max = vj;
		}
		return max;
	}

	public static void L2normalize( Matrix m ) {
		int i = m.rows;
		float sumSquares = 0.0f;
		while( --i >= 0 ) {
			int j = m.columns;
			while( --j >= 0 ) {
				sumSquares += m.entry[i][j]*m.entry[i][j];
			}
		}
		sumSquares = (float)Math.sqrt(sumSquares);
		i = m.rows;
		if(sumSquares == 0.0f)
			return;
		while( --i >= 0 ) {
			int j = m.columns;
			while( --j >= 0 ) {
				m.entry[i][j] /= sumSquares;
			}
		}
	}

	public static void L2normalize( float [] v ) {
		int j = v.length;
		float sumSquares = 0.0f;
		while(--j >= 0) sumSquares += v[j]*v[j];
		j = v.length;
		sumSquares = (float)Math.sqrt(sumSquares);
		if(sumSquares == 0.0f)
			return;
		while(--j >= 0) v[j] /= sumSquares;
	}

	public static float L2norm( Matrix m ) {
		int i = m.rows;
		float sumSquares = 0.0f;
		while( --i >= 0 ) {
			int j = m.columns;
			while( --j >= 0 ) {
				sumSquares += m.entry[i][j];
			}
		}
		sumSquares = (float)Math.sqrt(sumSquares);
		if(m.norms == null)
			m.norms = new float[4];
		m.norms[2] = sumSquares;
		return sumSquares;
	}

	public static float L2norm( float [] v ) {
		int j = v.length;
		float sumSquares = 0.0f;
		while(--j >= 0) sumSquares += v[j]*v[j];
		return (float)Math.sqrt(sumSquares);
	}

	public static void L1normalize( Matrix m ) {
		int i = m.rows;
		float sum = 0.0f;
		while( --i >= 0 ) {
			int j = m.columns;
			while( --j >= 0 ) {
				sum += Math.abs(m.entry[i][j]);
			}
		}
		i = m.rows;
		if(sum == 0.0f)
			return;
		while( --i >= 0 ) {
			int j = m.columns;
			while( --j >= 0 ) {
				m.entry[i][j] /= sum;
			}
		}
	}

	public static void L1normalize( float[] v ) {
		int j = v.length;
		float sum = 0.0f;
		while(--j >= 0) sum += Math.abs(v[j]);
		j = v.length;
		if(sum == 0.0f)
			return;
		while(--j >= 0) v[j] /= sum;
	}

	public static float L1norm( Matrix m ) {
		int i = m.rows;
		float sum = 0.0f;
		while( --i >= 0 ) {
			int j = m.columns;
			while( --j >= 0 ) {
				sum += Math.abs(m.entry[i][j]);
			}
		}
		if(m.norms == null)
			m.norms = new float[4];
		m.norms[1] = sum;
		return sum;
	}

	public static float L1norm( float[] v ) {
		int j = v.length;
		float sum = 0.0f;
		while(--j >= 0) sum += Math.abs(v[j]);
		return sum;
	}

	public static void L0normalize( Matrix m ) {
		int i = m.rows;
		float sum = 0.0f;
		while( --i >= 0 ) {
			int j = m.columns;
			while( --j >= 0 ) {
				if(m.entry[i][j] != Matrix.ZEROF)
					sum += 1.0f;
			}
		}
		i = m.rows;
		if(sum == 0.0f)
			return;
		while( --i >= 0 ) {
			int j = m.columns;
			while( --j >= 0 ) {
				m.entry[i][j] /= sum;
			}
		}
	}

	public static void L0normalize( float[] v ) {
		int j = v.length;
		float sum = 0.0f;
		while(--j >= 0) 
			if(v[j] != 0.0f)
				sum += 1.0f;
		j = v.length;
		if(sum == 0.0f)
			return;
		while(--j >= 0) v[j] /= sum;
	}

	public static float L0norm( Matrix m ) {
		int i = m.rows;
		float sum = 0.0f;
		while( --i >= 0 ) {
			int j = m.columns;
			while( --j >= 0 ) {
				if(m.entry[i][j] != Matrix.ZEROF)
					sum += 1.0f;
			}
		}
		if(m.norms == null)
			m.norms = new float[4];
		m.norms[0] = sum;
		return sum;
	}

	public static float L0norm( float[] v ) {
		int j = v.length;
		float sum = 0.0f;
		while(--j >= 0) 
			if(v[j] != 0.0f)
				sum += 1.0f;
		return sum;
	}

	// Algorithms

	public static Matrix FloydWarshall( Matrix connectivityMetric, float epsilon, float connected_bound, float disconnected_bound, int strictness, int minmax ) { 
		epsilon = Math.abs(epsilon);
		float connected = connected_bound;
		float disconnected = disconnected_bound;
		int i = connectivityMetric.rows;
		int maxspan = connectivityMetric.shortest;
		Matrix fw = new Matrix(Matrix.NONE, i, connectivityMetric.columns);
		while(--i >= 0) {
			int j = connectivityMetric.columns;
			int max = (minmax == MIN) ? 0 : i;
			while(--j >= max ) {
				float eij = connectivityMetric.entry[i][j];
				if( strictness == 0 && j < maxspan && i < maxspan ) {
					eij += connectivityMetric.entry[j][i];
				}
				if( eij != Matrix.ZEROF && eij >= -epsilon && eij <= epsilon )
					fw.entry[i][j] = connected;
				else
					fw.entry[i][j] = disconnected;
			}
		}
		int k = maxspan;
		while( --k >= 0 ) {      
			i = connectivityMetric.rows;
			while( --i >= 0 ) {
				float current_distance_part = fw.entry[i][k];
				int j = connectivityMetric.columns;
				while( -- j >= 0 ) {
					if(i==j)
						fw.entry[i][j] = 0.0f;
					else {
						float current_distance = fw.entry[i][j];
						float current_distance_whole;
						switch(minmax) {
						case(MIN):
							current_distance_whole = current_distance_part + fw.entry[k][j];
						fw.entry[i][j] = current_distance < current_distance_whole ? current_distance : current_distance_whole;
						break;
						case(MAX):
							current_distance_whole = current_distance_part < fw.entry[k][j] ? current_distance_part : fw.entry[k][j];
						fw.entry[i][j] = current_distance > current_distance_whole ? current_distance : current_distance_whole;
						break;
						}
					}
				}  
			}
		}
		i = connectivityMetric.rows;
		while( --i >= 0 ) {
			int j = connectivityMetric.columns;
			while( --j >= 0 ) 
				fw.entry[i][j] *= fw.entry[i][j];
		}
		return fw;
	}

	public static void PowerIterationDeflate(Matrix m, int eigens_to_compute, double epsilons, int populateKrylovSubspace ) { 
		int iteration = 0, i = -1;
		int maxIteration = m.shortest;
		m.krylovs = new float[populateKrylovSubspace][m.columns];
		while(++i < eigens_to_compute ) {
			Matrix next_eigen_approximation, current_eigen_approximation = new Matrix(Matrix.NONE, 1, m.columns );
			double prev, current = epsilons; float eigenvalue_approximation;
			fill_random( current_eigen_approximation );
			L2normalize( current_eigen_approximation );
			int iIteration = 0;
			do {
				prev = current;
				next_eigen_approximation = new Matrix(current_eigen_approximation);
				current_eigen_approximation = productResult(m, next_eigen_approximation.transpose()).transpose();
				if(iteration<populateKrylovSubspace)
					System.arraycopy(current_eigen_approximation.entry[0],0,m.krylovs[iteration++],0,m.columns);
				eigenvalue_approximation = FrobeniusInnerProduct(next_eigen_approximation, current_eigen_approximation);
				L2normalize(current_eigen_approximation);
				current = FrobeniusInnerProduct(next_eigen_approximation,current_eigen_approximation);
			} while( (current < (1-epsilons) && Math.abs( current/prev ) > (1 + epsilons)) && iIteration++ < maxIteration );
			m.eigenvalue[i] = eigenvalue_approximation;
			Deflate(m, current_eigen_approximation, eigenvalue_approximation );
			scalarMultiply(current_eigen_approximation, (float)Math.sqrt(Math.abs(eigenvalue_approximation)));
			int shortest = current_eigen_approximation.entry[0].length < m.eigenvector[i].length ? current_eigen_approximation.entry[0].length :  m.eigenvector[i].length;
			System.arraycopy( current_eigen_approximation.entry[0], 0, m.eigenvector[i], 0, shortest );
		}
	}

	public static void show_v (float[] v) {
		int i = v.length;
		while(--i >= 0)
			System.out.print(v[i] + " ");
		System.out.println();
	}

	public static void Deflate( Matrix base, Matrix eigen_vector, float eigen_value ) {
		Matrix r1 = scalarMultiplyResult(eigen_vector, eigen_value);
		Matrix result = productResult(  eigen_vector.transpose(), r1 );
		pairwiseSubtract( base, result );
	}

	public static Matrix DeflateResult( Matrix base, Matrix eigen_vector, float eigen_value ) {
		Matrix result = productResult(  eigen_vector, scalarMultiplyResult(eigen_vector, eigen_value).transpose() );
		result = pairwiseSubtractResult( base, result );
		return result;
	}

	public static void LU_Dolittle_Decomposition( Matrix m ) { 
		int rows = m.rows; int columns = m.columns;
		for(int i = 0; i < rows;i++ ) {
			for(int j = i; j < columns; j++ ) {
				for(int k = 0; k < i-1; k++ )
					m.entry[i][j] -= m.entry[i][k]*m.entry[k][j];
			}
			for(int j = i+1; j < rows; j++ ) {
				for(int k = 0; k < i-1; k++ ) 
					m.entry[j][i] -= m.entry[j][k]*m.entry[k][i];
				if(m.entry[i][i] != 0.0f)
					m.entry[j][i] /= m.entry[i][i];
			}
		}
	}

	public static Matrix LU_Dolittle_DecompositionResult( Matrix m ) { 
		int rows = m.rows; int columns = m.columns;
		Matrix lowerupper = new Matrix(Matrix.IDENTITY, rows, columns );
		for(int i = 0; i < rows;i++ ) {
			for(int j = i; j < columns; j++ ) {
				lowerupper.entry[i][j] = m.entry[i][j];
				for(int k = 0; k < i-1; k++ )
					lowerupper.entry[i][j] -= lowerupper.entry[i][k]*lowerupper.entry[k][j];
			}
			for(int j = i+1; j < rows; j++ ) {
				lowerupper.entry[j][i] = m.entry[j][i];
				for(int k = 0; k < i-1; k++ ) 
					lowerupper.entry[j][i] -= lowerupper.entry[j][k]*lowerupper.entry[j][i];
				if(lowerupper.entry[i][i] != 0.0f)
					lowerupper.entry[j][i] /= lowerupper.entry[i][i];
			}
		}
		return lowerupper;
	}

	public static Matrix[] LU_Dolittle_DecompositionResultPair( Matrix m ) { 
		int rows = m.rows; int columns = m.columns;
		Matrix lower = new Matrix(Matrix.IDENTITY, rows, columns );
		Matrix upper = new Matrix(Matrix.NONE, rows, columns );
		for(int i = 0; i < rows;i++ ) {
			for(int j = i; j < columns; j++ ) {
				upper.entry[i][j] = m.entry[i][j];
				for(int k = 0; k < i-1; k++ )
					upper.entry[i][j] -= lower.entry[i][k]*upper.entry[k][j];
			}
			for(int j = i+1; j < rows; j++ ) {
				lower.entry[j][i] = m.entry[j][i];
				for(int k = 0; k < i-1; k++ ) 
					lower.entry[j][i] -= lower.entry[j][k]*upper.entry[j][i];
				if(lower.entry[i][i] != 0.0f)
					lower.entry[j][i] /= upper.entry[i][i];
			}
		}
		return new Matrix[] { lower, upper };
	}

	public static float diagonalDet( Matrix triangular ) {
		float det = 1.0f;
		int i = triangular.shortest;
		while(--i >= 0 )
			det *= triangular.entry[i][i];
		triangular.det = det;
		return det;
	}

	public static float Determinant( Matrix m ) { 
		Matrix[] lowerupper = LU_Dolittle_DecompositionResultPair(m);
		float detLower = diagonalDet(lowerupper[0]);
		float detUpper = diagonalDet(lowerupper[1]);
		m.det = detLower*detUpper;
		return m.det;
	}

	public static Matrix Inverse( Matrix a ) {
		Matrix m = new Matrix(a);
		Matrix w = new Matrix( Matrix.NONE, m.rows, m.columns );
		Matrix l = new Matrix( Matrix.NONE, m.rows, m.columns );
		Matrix u = new Matrix( Matrix.NONE, m.rows, m.columns );
		for(int k = 0; k < m.shortest; k++ ) {
			l.entry[k][k] = Matrix.IDENTITY;
			for(int i = k+1; i < m.shortest; i++ ) {
				float coeff = Matrix.IDENTITYF; 
				if(m.entry[k][k] != 0.0f)
					coeff = m.entry[i][k]/m.entry[k][k];
				l.entry[i][k] = coeff;
				for(int j = k+1; j < m.shortest; j++ )
					m.entry[i][j] -= coeff*m.entry[k][j];
			}
		}
		for(int j = 0; j < m.shortest; j++ ) {
			for(int i = 0; i < j; i++ ) {
				u.entry[i][j] = m.entry[i][j];
			}
		}
		float b[] = new float[m.shortest];
		float d[] = new float[m.shortest];
		float x[] = new float[m.shortest];
		for(int k = 0; k < m.columns; k++ ) {

			for(int i = 1; i < m.shortest; i++ ) {
				b[k] = 1.0f;
				d[1] = b[1];
				for(int j = 1; j < i-1; j++ ) {
					d[i] = b[i];
					d[i] -= l.entry[i][j]*d[j];
				}
			}
			if(u.entry[m.shortest-1][m.shortest-1] != 0.0f)
				x[m.shortest-1] = d[m.shortest-1]/u.entry[m.shortest-1][m.shortest-1];
			int i = m.shortest;
			while(--i >= 0 ) {
				x[i] = d[i];
				int j = m.shortest;
				while(--j > i ) {
					x[i] = x[i] - u.entry[i][j]*x[j];
				}
				if(u.entry[i][i] != 0.0f) 
					x[i] /= u.entry[i][i]; 
			}
			for(i = 0; i < m.shortest; i++ ) {
				w.entry[i][k] = x[i];
			}
			b[k] = 0.0f;
		}
		return w;
	}

	// Thresholding

	public static void threshold ( Matrix a, float threshold, int directed ) {
		int shortestRow = a.rows;
		int shortestColumn = a.columns ;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float aij = a.entry[i][j];
				if(directed==0)
					aij = Math.abs(a.entry[i][j]);
				if(aij < threshold)
					a.entry[i][j] = 0.0f;
				else
					a.entry[i][j] = 1.0f;
			}  
		}
		a.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix thresholdResult ( Matrix a, float threshold, int directed ) {
		int shortestRow = a.rows;
		int shortestColumn = a.columns;
		Matrix result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float aij = a.entry[i][j];
				if(directed==0)
					aij = Math.abs(a.entry[i][j]);
				if(aij < threshold)
					result.entry[i][j] = 0.0f;
				else
					result.entry[i][j] = 1.0f;
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	// Pair/Entry-Wise Fuzzy Binary arithmetic
	// a AND 0 = 0
	// a or 0 = a
	// a XOR 0 = a
	// a XOR a = 0
	// a AND b = a*b
	// a OR b = a+b
	// a XOR b = (a/b + b/a)
	// a KOR b = -ab + b/a 
	// NOT a = 0
	// a NAND B = a/b
	// a NOR b = a + 1/b

	public static void AND( Matrix a, Matrix b ) {
		int shortestRow = a.rows < b.rows ? a.rows : b.rows;
		int shortestColumn = a.columns < b.columns ? a.columns : b.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				a.entry[i][j] = a.entry[i][j]*b.entry[i][j];
			}  
		}
		a.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix ANDResult ( Matrix a, Matrix b ) {
		int shortestRow = a.rows < b.rows ? a.rows : b.rows;
		int shortestColumn = a.columns < b.columns ? a.columns : b.columns;
		Matrix result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				result.entry[i][j] = a.entry[i][j]*b.entry[i][j];
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static void OR( Matrix a, Matrix b ) {
		int shortestRow = a.rows < b.rows ? a.rows : b.rows;
		int shortestColumn = a.columns < b.columns ? a.columns : b.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				a.entry[i][j] = a.entry[i][j]+b.entry[i][j];
			}  
		}
		a.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix ORResult ( Matrix a, Matrix b ) {
		int shortestRow = a.rows < b.rows ? a.rows : b.rows;
		int shortestColumn = a.columns < b.columns ? a.columns : b.columns;
		Matrix result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				result.entry[i][j] = a.entry[i][j]+b.entry[i][j];
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static void XOR( Matrix a, Matrix b ) {
		int shortestRow = a.rows < b.rows ? a.rows : b.rows;
		int shortestColumn = a.columns < b.columns ? a.columns : b.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float aij = a.entry[i][j];
				float bij = b.entry[i][j];
				if(aij == 0) {
					if(bij == 0 )
						a.entry[i][j] = 0.0f;
					else
						a.entry[i][j] = bij;
				}
				else {
					if(bij != 0 )
						a.entry[i][j] = aij/bij + bij/aij;
				}
			}  
		}
		a.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix XORResult ( Matrix a, Matrix b ) {
		int shortestRow = a.rows < b.rows ? a.rows : b.rows;
		int shortestColumn = a.columns < b.columns ? a.columns : b.columns;
		Matrix result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float aij = a.entry[i][j];
				float bij = b.entry[i][j];
				if(aij == 0) {
					if(bij == 0 )
						result.entry[i][j] = 0.0f;
					else
						result.entry[i][j] = bij;
				}
				else {
					if(bij != 0 )
						result.entry[i][j] = aij/bij + bij/aij;
				}
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static void KOR( Matrix a, Matrix b ) {
		int shortestRow = a.rows < b.rows ? a.rows : b.rows;
		int shortestColumn = a.columns < b.columns ? a.columns : b.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float aij = a.entry[i][j];
				float bij = b.entry[i][j];
				if(aij != 0)
					a.entry[i][j] = -aij*bij + bij/aij;
			}
		}
		a.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix KORResult ( Matrix a, Matrix b ) {
		int shortestRow = a.rows < b.rows ? a.rows : b.rows;
		int shortestColumn = a.columns < b.columns ? a.columns : b.columns;
		Matrix result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float aij = a.entry[i][j];
				float bij = b.entry[i][j];
				if(aij != 0)
					result.entry[i][j] = -aij*bij + bij/aij;
			}
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static void NOT ( Matrix a ) {
		int shortestRow = a.rows;
		int shortestColumn = a.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float aij = a.entry[i][j];
				if(aij != 0)
					a.entry[i][j] = 1/-aij;
			}
		}
		a.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix NOTResult ( Matrix a ) {
		int shortestRow = a.rows;
		int shortestColumn = a.columns;
		Matrix result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float aij = a.entry[i][j];
				if(aij != 0)
					result.entry[i][j] = 1.0f/-aij;
			}
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	// General

	public static Matrix complement ( Matrix m ) {
		float[] mm = m.minmax();
		float min = mm[0];
		float max = mm[1];
		float avg = (min+max)/2.0f;
		Matrix result = new Matrix(Matrix.NONE, m.rows, m.columns);
		int shortestRow = m.rows;
		int shortestColumn = m.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float mij = m.entry[i][j];
				if(mij < avg)
					result.entry[i][j] = max;
				else
					result.entry[i][j] = min;
			}
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	// Modular arithmetic 

	//Addition

	public static void ModscalarAdd(Matrix accumulator, float addend, float modulus ) {
		if(addend == Matrix.ZEROF || modulus == Matrix.ZEROF) {
			accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int shortestRow = accumulator.rows;
		int shortestColumn = accumulator.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				accumulator.entry[i][j] += addend;
				accumulator.entry[i][j] %= modulus;
			}  
		}
		accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix ModscalarAddResult(Matrix adder, float addend, float modulus ) {
		Matrix result;
		if(addend == Matrix.ZEROF || modulus == Matrix.ZEROF) {
			result = new Matrix(adder);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		int shortestRow = adder.rows;
		int shortestColumn = adder.columns;
		result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				result.entry[i][j] = adder.entry[i][j] + addend;
				result.entry[i][j] %= modulus;
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static void ModpairwiseSum(Matrix accumulator, Matrix summand, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int shortestRow = accumulator.rows < summand.rows ? accumulator.rows : summand.rows;
		int shortestColumn = accumulator.columns < summand.columns ? accumulator.columns : summand.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				accumulator.entry[i][j] += summand.entry[i][j];
				accumulator.entry[i][j] %= modulus;
			}  
		}
		accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix ModpairwiseSumResult(Matrix adder, Matrix summand, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			adder.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return adder;
		}int shortestRow = adder.rows < summand.rows ? adder.rows : summand.rows;
		int shortestColumn = adder.columns < summand.columns ? adder.columns : summand.columns;
		Matrix result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				result.entry[i][j] = adder.entry[i][j] + summand.entry[i][j];
				result.entry[i][j] %= modulus;
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	// Subtraction

	public static void ModscalarSubtract(Matrix accumulator, float subtractand, float modulus ) {
		if(modulus == Matrix.ZEROF || subtractand == Matrix.ZEROF) {
			accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int shortestRow = accumulator.rows;
		int shortestColumn = accumulator.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				accumulator.entry[i][j] -= subtractand;
				accumulator.entry[i][j] %= modulus;
			}  
		}
		accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix ModscalarSubtractResult(Matrix subtractor, float subtractand, float modulus ) {
		Matrix result;
		if(modulus == Matrix.ZEROF || subtractand == Matrix.ZEROF) {
			result = new Matrix(subtractor);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		int shortestRow = subtractor.rows;
		int shortestColumn = subtractor.columns;
		result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				result.entry[i][j] = subtractor.entry[i][j] - subtractand;
				result.entry[i][j] %= modulus;
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static void ModpairwiseSubtract(Matrix accumulator, Matrix subtractand, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int shortestRow = accumulator.rows < subtractand.rows ? accumulator.rows : subtractand.rows;
		int shortestColumn = accumulator.columns < subtractand.columns ? accumulator.columns : subtractand.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				accumulator.entry[i][j] -= subtractand.entry[i][j];
				accumulator.entry[i][j] %= modulus;
			}  
		}
		accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix ModpairwiseSubtractResult(Matrix subtractor, Matrix subtractand, float modulus ) {
		Matrix result;
		if(modulus == Matrix.ZEROF) {
			result = new Matrix(subtractor);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		int shortestRow = subtractor.rows < subtractand.rows ? subtractor.rows : subtractand.rows;
		int shortestColumn = subtractor.columns < subtractand.columns ? subtractor.columns : subtractand.columns;
		result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				result.entry[i][j] = subtractor.entry[i][j] - subtractand.entry[i][j];
				result.entry[i][j] %= modulus;
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	// Multiplication

	public static void ModscalarMultiply(Matrix accumulator, float multiplier, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int shortestRow = accumulator.rows;
		int shortestColumn = accumulator.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				accumulator.entry[i][j] *= multiplier;
				accumulator.entry[i][j] %= modulus;
			}  
		}
		accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix ModscalarMultiplyResult(Matrix multiplicand, float multiplier, float modulus ) {
		Matrix result;
		if(modulus == Matrix.ZEROF) {
			result = new Matrix(multiplicand);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		int shortestRow = multiplicand.rows;
		int shortestColumn = multiplicand.columns;
		result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				result.entry[i][j] = multiplicand.entry[i][j] * multiplier;
				result.entry[i][j] %= modulus;
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static float ModFrobeniusInnerProduct( Matrix multiplicand, Matrix multiplier, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			return modulus;
		}
		int shortestRow = multiplicand.rows < multiplier.rows ? multiplicand.rows : multiplier.rows;
		int shortestColumn = multiplicand.columns < multiplier.columns ? multiplicand.columns : multiplier.columns;
		int i = shortestRow;
		float accumulator = 0.0f;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				accumulator += multiplicand.entry[i][j]*multiplier.entry[i][j];
			}  
		}
		return accumulator % modulus;
	}

	public static void ModhadamardProduct(Matrix accumulator, Matrix multiplicand, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int shortestRow = accumulator.rows < multiplicand.rows ? accumulator.rows : multiplicand.rows;
		int shortestColumn = accumulator.columns < multiplicand.columns ? accumulator.columns : multiplicand.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				accumulator.entry[i][j] *= multiplicand.entry[i][j];
				accumulator.entry[i][j] %= modulus;
			}  
		}
		accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix ModhadamardProductResult(Matrix multiplier, Matrix multiplicand, float modulus ) {
		Matrix result;
		if(modulus == Matrix.ZEROF) {
			result = new Matrix(multiplicand);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		int shortestRow = multiplier.rows < multiplicand.rows ? multiplier.rows : multiplicand.rows;
		int shortestColumn = multiplier.columns < multiplicand.columns ? multiplier.columns : multiplicand.columns;
		result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				result.entry[i][j] = multiplier.entry[i][j] * multiplicand.entry[i][j];
				result.entry[i][j] %= modulus;
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static void ModproductPropagated(Matrix accumulator, Matrix multiplier, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int shortestRow = accumulator.rows < multiplier.rows ? accumulator.rows : multiplier.rows;
		int shortestColumn = accumulator.columns < multiplier.columns ? accumulator.columns : multiplier.columns;
		int shortest = shortestRow < shortestColumn ? shortestRow : shortestColumn;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				int k = shortest;
				float entrySum = 0.0f;
				while(--k >= 0) {
					entrySum += accumulator.entry[i][k]*multiplier.entry[k][j];
				}
				accumulator.entry[i][j] = entrySum;
				accumulator.entry[i][j] %= modulus;
			}  
		}
		accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix ModproductResult(Matrix multiplicand, Matrix multiplier, float modulus ) { 
		Matrix result;
		if(modulus == Matrix.ZEROF) {
			result = new Matrix(multiplicand);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		int shortest = multiplicand.columns < multiplier.rows ? multiplicand.columns : multiplier.rows;
		int i = multiplicand.rows;
		int j = multiplier.columns;
		result = new Matrix(Matrix.NONE, i, j);
		while(--i >= 0) {
			j = multiplier.columns;
			while(--j >= 0) {
				int k = shortest;
				float entrySum = 0.0f;
				while(--k >= 0) {
					entrySum += multiplicand.entry[i][k]*multiplier.entry[k][j];
				}
				result.entry[i][j] = entrySum;
				result.entry[i][j] %= modulus;
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	// Division

	public static void ModentrywiseInversion( Matrix divisor, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			divisor.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int shortestRow = divisor.rows;
		int shortestColumn = divisor.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float divisorij = divisor.entry[i][j];
				if(divisorij != Matrix.ZEROF)
					divisor.entry[i][j] = Matrix.IDENTITYF / divisorij;
				else
					divisor.entry[i][j] = Matrix.IDENTITYF;
				divisor.entry[i][j] %= modulus;
			}  
		}
	}

	public static Matrix ModentrywiseInversionResult ( Matrix divisor, float modulus ) {
		int shortestRow = divisor.rows;
		int shortestColumn = divisor.columns;
		Matrix result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float divisorij = divisor.entry[i][j];
				if(divisorij != Matrix.ZEROF)
					result.entry[i][j] = Matrix.IDENTITYF / divisorij;
				else
					result.entry[i][j] = Matrix.IDENTITYF;
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static void ModscalarDivide(Matrix accumulator, float divisor, float modulus) {
		if(modulus == Matrix.ZEROF) {
			accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int shortestRow = accumulator.rows;
		int shortestColumn = accumulator.columns;
		if(divisor == Matrix.ZEROF) {
			Arrays.fill(accumulator.entry, Matrix.IDENTITYF);
			accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) 
				accumulator.entry[i][j] /= divisor;
		}
	}

	public static Matrix ModscalarDivideResult(Matrix dividend, float divisor, float modulus) {
		Matrix result;
		if(modulus == Matrix.ZEROF) {
			result = new Matrix(dividend);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		int shortestRow = dividend.rows;
		int shortestColumn = dividend.columns;
		result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				if(divisor != Matrix.ZEROF)
					result.entry[i][j] = dividend.entry[i][j] / divisor;
				else
					result.entry[i][j] = Matrix.IDENTITYF;
				result.entry[i][j] %= modulus;
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static void ModhadamardDivision(Matrix accumulator, Matrix divisor, float modulus) {
		if(modulus == Matrix.ZEROF) {
			accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int shortestRow = accumulator.rows < divisor.rows ? accumulator.rows : divisor.rows;
		int shortestColumn = accumulator.columns < divisor.columns ? accumulator.columns : divisor.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				if(divisor.entry[i][j] != Matrix.ZEROF)
					accumulator.entry[i][j] /= divisor.entry[i][j];
				else
					accumulator.entry[i][j] = Matrix.ZEROF;
				accumulator.entry[i][j] %= modulus;
			}  
		}
		accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix ModhadamardDivisionResult(Matrix dividend, Matrix divisor, float modulus) {
		Matrix result;
		if(modulus == Matrix.ZEROF) {
			result = new Matrix(dividend);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		int shortestRow = dividend.rows < divisor.rows ? dividend.rows : divisor.rows;
		int shortestColumn = dividend.columns < divisor.columns ? dividend.columns : divisor.columns;
		result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				if(divisor.entry[i][j] != Matrix.ZEROF)
					result.entry[i][j] = dividend.entry[i][j] / divisor.entry[i][j];
				else
					result.entry[i][j] = Matrix.ZEROF;
				result.entry[i][j] %= modulus;
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	// Powers

	public static void ModhadamardPower(Matrix accumulator, int power, float modulus) {
		if(modulus == Matrix.ZEROF) {
			accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int shortestRow = accumulator.rows;
		int shortestColumn = accumulator.columns;
		if(power == 0) {
			accumulator.make_ones();
			return;  
		}
		if(power == 1)
			return;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				int p = power;
				while(--p > 0) {
					accumulator.entry[i][j] *= accumulator.entry[i][j];
					accumulator.entry[i][j] %= modulus;
				}
			}  
		}
		accumulator.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix ModhadamardPowerResult(Matrix multiplier, int power, float modulus) {
		Matrix result;
		if(modulus == Matrix.ZEROF) {
			result = new Matrix(multiplier);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		int shortestRow = multiplier.rows;
		int shortestColumn = multiplier.columns;
		if(power == 1)
			return new Matrix(multiplier);
		result = new Matrix(Matrix.ONES, shortestRow, shortestColumn);
		if(power == 0)
			return result;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				int p = power;
				while(--p >= 0) {
					result.entry[i][j] *= multiplier.entry[i][j];
					result.entry[i][j] %= modulus;
				}
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static int[] Modradix(int num, int base, float modulus) {
		if(modulus == Matrix.ZEROF) {
			return new int[10];
		}
		int b = 1; int b1 = base;
		int rem = 0; int rem1 = num % b1;
		int bit = 0;
		int[] bits = new int[10];
		while(b <= num) {
			bits[bit] = ((rem1 - rem)/b) % (int)modulus;
			b=b1;b1*=base;
			rem=rem1; rem1 = num % b1;
			bit++;
			if(bit >= bits.length) {
				int[] temp = new int[bit<<1];
				System.arraycopy(bits, 0, temp, 0, bit);
				bits = temp;
			}
		}
		return bits;
	}

	public static float ModexponentialTrace( Matrix m, float modulus ) {
		if(modulus == Matrix.ZEROF)
			return modulus;
		float result = 1.0f;
		float trace = m.trace();
		int itrace = (int)trace;
		float diff = trace - itrace;
		int[] bits = radix(itrace, 2);
		int iteration = 0;
		float lastPow = 2.71828183f;
		int lastPowIndex = 0;
		while(iteration++ < bits.length) {
			int bit = bits[iteration];
			if(bit == 1) {
				int squaresRequired = iteration - lastPowIndex;
				while(--squaresRequired >= 0)
					lastPow = (lastPow * lastPow) % modulus;
				result = (result * lastPow) % modulus;
			}
		}
		float extra = -(float)Math.log(diff);
		if(extra != 0.0f)
			result *= extra;
		return result % modulus;
	}

	public static void mod(Matrix base, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			base.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int shortestRow = base.rows;
		int shortestColumn = base.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				base.entry[i][j] %= modulus;
			}  
		}
		base.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix modResult( Matrix base, float modulus ) {
		Matrix result;
		if(modulus == Matrix.ZEROF) {
			result = new Matrix(base);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		int shortestRow = base.rows;
		int shortestColumn = base.columns;
		result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0)
				result.entry[i][j] %= modulus;
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static void ModpolysquarePropagated(Matrix base, int times, float modulus) {
		if(modulus == Matrix.ZEROF) {
			base.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		while(--times >= 0) {
			productPropagated(base, base);
			mod(base, modulus);
		}
	}

	public static void ModfastPropagatedMultiproduct(Matrix base, int exponent, float modulus) { 
		if(modulus == Matrix.ZEROF) {
			base.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int[] operations = radix(exponent, 2);
		int lastProductIndex = 0, jmax = operations.length;
		Matrix partialPowers;
		if(operations[0] == 0)
			partialPowers = new Matrix(Matrix.IDENTITY, base.rows, base.columns);
		else
			partialPowers = new Matrix(base);
		for(int j = 1; j < jmax;j++) {
			int bit = operations[j];
			if(bit == 1) {
				ModpolysquarePropagated(partialPowers, j-lastProductIndex, modulus);
				lastProductIndex = j;
				ModproductPropagated(base, partialPowers, modulus);
				mod(base,modulus);
			}
		}
	}

	public static Matrix Modpolysquare(Matrix base, int times, float modulus) {
		if(modulus == Matrix.ZEROF) {
			Matrix result = new Matrix(base);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		while(--times >= 0) {
			base = productResult(base, base);
			base = modResult(base, modulus);
		}
		base.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return base;
	}

	public static Matrix ModfastMultiproduct(Matrix base, int exponent, float modulus) { 
		if(modulus == Matrix.ZEROF) {
			Matrix result = new Matrix(base);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		int[] operations = radix(exponent, 2);
		int lastProductIndex = 0, jmax = operations.length;
		Matrix partialPowers;
		if(operations[0] == 0)
			partialPowers = new Matrix(Matrix.IDENTITY, base.rows, base.columns);
		else
			partialPowers = new Matrix(base);
		for(int j = 1; j < jmax;j++) {
			int bit = operations[j];
			if(bit == 1) {
				Modpolysquare(partialPowers, j-lastProductIndex, modulus);
				lastProductIndex = j;
				base = ModproductResult(base, partialPowers, modulus);
				base = modResult(base,modulus);
			}
		}
		return base;
	}

	public static Matrix ModsquareRootReturnInverseSquareRoot(Matrix y, float modulus) {
		if(modulus == Matrix.ZEROF) {
			y.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return ModInverse(y, modulus);
		}
		Matrix z = new Matrix(Matrix.IDENTITY, y.rows, y.columns);
		int j = (int)Math.sqrt(y.rows) + 1;
		while(--j >= 0) {
			scalarMultiply(pairwiseSumResult(y,Inverse(z)), 0.5f);
			scalarMultiply(pairwiseSumResult(z,Inverse(y)), 0.5f);
		}
		y.updates &= (Matrix.all_changed^Matrix.entry_changed);
		z.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return z;
	}

	public static Matrix[] ModsquareRootInverseSquareRootResult(Matrix x, float modulus) {
		if(modulus == Matrix.ZEROF) {
			Matrix result = new Matrix(x);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return new Matrix[] {result, ModInverse(result, modulus)};
		}
		Matrix y = new Matrix(x);
		Matrix z = new Matrix(Matrix.IDENTITY, y.rows, y.columns);
		int j = (int)Math.sqrt(y.rows) + 1;
		while(--j >= 0) {
			scalarMultiply(pairwiseSumResult(y,Inverse(z)), 0.5f);
			scalarMultiply(pairwiseSumResult(z,Inverse(y)), 0.5f);
		}
		y.updates &= (Matrix.all_changed^Matrix.entry_changed);
		z.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return new Matrix[] { y, z };
	}

	public static void ModentrywiseSquareRoot(Matrix x, float modulus) {
		if(modulus == Matrix.ZEROF)
			return;
		int shortestRow = x.rows;
		int shortestColumn = x.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float xij = x.entry[i][j] % modulus;
				float sign = xij/Math.abs(xij);
				x.entry[i][j] = (float)( sign*Math.sqrt(Math.abs(xij)) );
			}  
		}
		x.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix ModentrywiseSquareRootResult(Matrix x, float modulus) {
		if(modulus == Matrix.ZEROF)
			return new Matrix(x);
		int shortestRow = x.rows;
		int shortestColumn = x.columns;
		int i = shortestRow;
		Matrix result = new Matrix(Matrix.NONE, x.rows, x.columns);
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float xij = x.entry[i][j] % modulus;
				float sign = xij/Math.abs(xij);
				result.entry[i][j] = (float)( sign*Math.sqrt(Math.abs(xij)) );
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	// Logarithm

	public static void ModentrywiseLogarithm( Matrix m, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			return;
		}
		int shortestRow = m.rows;
		int shortestColumn = m.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float eij = m.entry[i][j];
				float sign = eij/Math.abs(eij);
				m.entry[i][j] = (float)(sign*Math.log(eij));
			}  
		}
		m.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix ModentrywiseLogarithmResult( Matrix m, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			m.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return m;
		}
		int shortestRow = m.rows;
		int shortestColumn = m.columns;
		int i = shortestRow;
		Matrix result = new Matrix(m);
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float eij = m.entry[i][j];
				float sign = eij/Math.abs(eij);
				result.entry[i][j] = (float)(sign*Math.log(eij));
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}


	// Mapping 

	public static void ModscaleToRange(Matrix m, float low, float high, float modulus) {
		if(modulus == Matrix.ZEROF) {
			m.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		float[] minmax = m.minmax();
		float newrange = high-low;
		float oldrange = minmax[1]-minmax[0];
		float range = oldrange == 0.0f ? 1.0f : newrange/oldrange;
		mod(m,modulus);
		scalarSubtract(m, minmax[0]);
		scalarMultiply(m, range );
		scalarAdd(m, low );

	}

	public static Matrix ModscaleToRangeResult(Matrix m, float low, float high, float modulus) {
		if(modulus == Matrix.ZEROF) {
			Matrix result = new Matrix(m);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		float[] minmax = m.minmax();
		float newrange = high-low;
		float oldrange = minmax[1]-minmax[0];
		float range = oldrange == 0.0f ? 1.0f : newrange/oldrange;
		Matrix result = new Matrix(Matrix.NONE, m.rows, m.columns);
		result = ModscalarSubtractResult(m, minmax[0], modulus);
		ModscalarMultiply(result, range, modulus );
		ModscalarAdd(result, low, modulus );
		return result;
	}


	// Random

	public static void Modfill_random(Matrix m, float modulus) {
		if(modulus == Matrix.ZEROF) {
			m.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int i = m.rows;
		while(--i >= 0 ) {
			int j = m.columns;
			while(--j >= 0 )
				m.entry[i][j] = (float)Math.random() % modulus;
		}
	}

	public static void Modfill_random( float[] v, float modulus) {
		if(modulus == Matrix.ZEROF) {
			return;
		}
		int j = v.length;
		while(--j >= 0 ) v[j] = (float)Math.random() % modulus;
	}


	// Norms & normalization vectors and matrices : 0, 1, 2, inf

	public static void ModLinfnormalize( Matrix m, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			m.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int i = m.rows;
		int j = m.columns;
		if( i == 0 || j == 0)
			return;
		float max = Math.abs(m.entry[0][0]);
		while( --i >= 0 ) {
			while( --j >= 0 ) {
				float eij = Math.abs(m.entry[i][j]);
				if(eij > max)
					max = eij;
			}
			j = m.columns;
		}
		i = m.rows;
		while( --i >= 0 ) {
			j = m.columns;
			while( --j >= 0 ) {
				m.entry[i][j] /= max;
			}
		}
	}

	public static void ModLinfnormalize( float[] v, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			return;
		}
		int j = v.length;
		if(j == 0)
			return;
		float max = Math.abs(v[0]);
		while(--j > 0) {
			float vj = Math.abs(v[j]);
			if(vj > max)
				max = vj;
		}
		j = v.length;
		while(--j >= 0) v[j] /= max;
	}

	public static float ModLinfnorm( Matrix m, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			m.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return modulus;
		}
		int i = m.rows;
		int j = m.columns;
		if( i == 0 || j == 0)
			return 0.0f;
		float max = Math.abs(m.entry[0][0]);
		while( --i >= 0 ) {
			while( --j >= 0 ) {
				float eij = Math.abs(m.entry[i][j]);
				if(eij > max)
					max = eij;
			}
			j = m.columns;
		}
		if(m.norms == null)
			m.norms = new float[4];
		m.norms[3] = max;
		return max;
	}

	public static float ModLinfnorm( float[] v, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			return modulus;
		}
		int j = v.length;
		if(j == 0)
			return 0.0f;
		float max = Math.abs(v[0]);
		while(--j > 0) {
			float vj = Math.abs(v[j]);
			if(vj > max)
				max = vj;
		}
		return max;
	}

	public static void ModL2normalize( Matrix m, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			m.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int i = m.rows;
		float sumSquares = 0.0f;
		while( --i >= 0 ) {
			int j = m.columns;
			while( --j >= 0 ) {
				sumSquares += m.entry[i][j]*m.entry[i][j];
			}
		}
		sumSquares = (float)Math.sqrt(sumSquares);
		i = m.rows;
		while( --i >= 0 ) {
			int j = m.columns;
			while( --j >= 0 ) {
				m.entry[i][j] /= sumSquares;
			}
		}
	}

	public static void ModL2normalize( float [] v, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			return;
		}
		int j = v.length;
		float sumSquares = 0.0f;
		while(--j >= 0) sumSquares += v[j]*v[j];
		j = v.length;
		sumSquares = (float)Math.sqrt(sumSquares);
		while(--j >= 0) v[j] /= sumSquares;
	}

	public static float ModL2norm( Matrix m, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			m.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return modulus;
		}
		int i = m.rows;
		float sumSquares = 0.0f;
		while( --i >= 0 ) {
			int j = m.columns;
			while( --j >= 0 ) {
				sumSquares += (m.entry[i][j]*m.entry[i][j])%modulus;
			}
		}
		sumSquares = (float)Math.sqrt(sumSquares);
		if(m.norms == null)
			m.norms = new float[4];
		m.norms[2] = sumSquares;
		return sumSquares;
	}

	public static float ModL2norm( float [] v, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			return modulus;
		}
		int j = v.length;
		float sumSquares = 0.0f;
		while(--j >= 0) sumSquares += (v[j]*v[j])%modulus;
		return (float)Math.sqrt(sumSquares);
	}

	public static void ModL1normalize( Matrix m, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			m.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int i = m.rows;
		float sum = 0.0f;
		while( --i >= 0 ) {
			int j = m.columns;
			while( --j >= 0 ) {
				sum += Math.abs(m.entry[i][j]);
			}
		}
		i = m.rows;
		while( --i >= 0 ) {
			int j = m.columns;
			while( --j >= 0 ) {
				m.entry[i][j] /= sum;
			}
		}
	}

	public static void ModL1normalize( float[] v, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			return;
		}
		int j = v.length;
		float sum = 0.0f;
		while(--j >= 0) sum += Math.abs(v[j]);
		j = v.length;
		while(--j >= 0) v[j] /= sum;
	}

	public static float ModL1norm( Matrix m, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			m.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return modulus;
		}
		int i = m.rows;
		float sum = 0.0f;
		while( --i >= 0 ) {
			int j = m.columns;
			while( --j >= 0 ) {
				sum += Math.abs(m.entry[i][j]);
			}
		}
		if(m.norms == null)
			m.norms = new float[4];
		m.norms[1] = sum;
		return sum;
	}

	public static float ModL1norm( float[] v, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			return modulus;
		}
		int j = v.length;
		float sum = 0.0f;
		while(--j >= 0) sum += Math.abs(v[j]);
		return sum;
	}

	public static void ModL0normalize( Matrix m, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			m.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int i = m.rows;
		float sum = 0.0f;
		while( --i >= 0 ) {
			int j = m.columns;
			while( --j >= 0 ) {
				if(m.entry[i][j] != Matrix.ZEROF)
					sum += 1.0f;
			}
		}
		i = m.rows;
		while( --i >= 0 ) {
			int j = m.columns;
			while( --j >= 0 ) {
				m.entry[i][j] /= sum;
			}
		}
	}

	public static void ModL0normalize( float[] v, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			return;
		}
		int j = v.length;
		float sum = 0.0f;
		while(--j >= 0) 
			if(v[j] != 0.0f)
				sum += 1.0f;
		j = v.length;
		while(--j >= 0) v[j] /= sum;
	}

	public static float ModL0norm( Matrix m, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			m.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return modulus;
		}
		int i = m.rows;
		float sum = 0.0f;
		while( --i >= 0 ) {
			int j = m.columns;
			while( --j >= 0 ) {
				if(m.entry[i][j] != Matrix.ZEROF)
					sum += 1.0f;
			}
		}
		if(m.norms == null)
			m.norms = new float[4];
		m.norms[0] = sum;
		return sum;
	}

	public static float ModL0norm( float[] v, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			return modulus;
		}
		int j = v.length;
		float sum = 0.0f;
		while(--j >= 0) 
			if(v[j] != 0.0f)
				sum += 1.0f;
		return sum;
	}

	// Algorithms

	public static Matrix ModFloydWarshall( Matrix base, float epsilon, float connected_bound, float disconnected_bound, int strictness, int minmax, float modulus ) { 
		if(modulus == Matrix.ZEROF) {
			Matrix result = new Matrix(base);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		epsilon = Math.abs(epsilon);
		float connected = connected_bound;
		float disconnected = disconnected_bound;
		int i = base.rows;
		int maxspan = base.shortest;
		Matrix fw = new Matrix(Matrix.NONE, i, base.columns);
		while(--i >= 0) {
			int j = base.columns;
			while(--j >= 0 ) {
				float eij = base.entry[i][j];
				if( strictness == 0 && j < maxspan && i < maxspan ) {
					eij += base.entry[j][i];
				}
				if( eij > -epsilon && eij < epsilon )
					fw.entry[i][j] = connected;
				else
					fw.entry[i][j] = disconnected;
			}
		}
		int k = maxspan;
		while( --k >= 0 ) {      
			i = base.rows;
			while( --i >= 0 ) {
				float current_distance_part = fw.entry[i][k];
				int j = base.columns;
				while( -- j >= 0 ) {
					if(i==j)
						fw.entry[i][j] = 0.0f;
					else {
						float current_distance = fw.entry[i][j];
						float current_distance_whole = current_distance_part + fw.entry[k][j];
						switch(minmax) {
						case(MIN):
							fw.entry[i][j] = current_distance < current_distance_whole ? current_distance : current_distance_whole;
						break;
						case(MAX):
							fw.entry[i][j] = current_distance > current_distance_whole ? current_distance : current_distance_whole;
							break;
						}
						fw.entry[i][j] %= modulus;
					}
				}  
			}
		}
		i = base.rows;
		while( --i >= 0 ) {
			int j = base.columns;
			while( --j >= 0 ) 
				fw.entry[i][j] *= fw.entry[i][j];
		}
		return fw;
	}

	public static void ModPowerIterationDeflate(Matrix m, int eigens_to_compute, double epsilons, int populateKrylovSubspace, float modulus ) { 
		if(modulus == Matrix.ZEROF) {
			m.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int iteration = 0;
		m.krylovs = new float[populateKrylovSubspace][m.columns];
		while(--eigens_to_compute >= 0) {
			Matrix next_eigen_approximation, current_eigen_approximation = new Matrix(Matrix.NONE, 1, m.columns );
			double prev, current = epsilons; float eigenvalue_approximation;
			fill_random( current_eigen_approximation );
			L2normalize( current_eigen_approximation );
			do {
				prev = current;
				next_eigen_approximation = new Matrix(current_eigen_approximation);
				current_eigen_approximation = ModproductResult(m, next_eigen_approximation, modulus);
				if(iteration<populateKrylovSubspace)
					System.arraycopy(current_eigen_approximation.entry[0],0,m.krylovs[iteration++],0,m.columns);
				eigenvalue_approximation = ModFrobeniusInnerProduct(next_eigen_approximation, current_eigen_approximation, modulus);
				ModL2normalize(current_eigen_approximation, modulus);
				current = ModFrobeniusInnerProduct(next_eigen_approximation,current_eigen_approximation, modulus);
			} while( (current < (1-epsilons) && Math.abs( current/prev ) > (1 + epsilons)) );
			ModscalarMultiply(current_eigen_approximation, (float)Math.sqrt(Math.abs(eigenvalue_approximation)), modulus);
			m.eigenvalue[eigens_to_compute] = eigenvalue_approximation;
			System.arraycopy( current_eigen_approximation.entry[0], 0,  m.eigenvector[eigens_to_compute],0, m.columns );
			if(eigens_to_compute > 0)
				m = ModDeflateResult(m, current_eigen_approximation, eigenvalue_approximation, modulus );
		}
	}

	public static Matrix ModDeflateResult( Matrix base, Matrix eigen_vector, float eigen_value, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			Matrix result = new Matrix(base);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		Matrix result = ModproductResult(  eigen_vector, scalarMultiplyResult(eigen_vector, eigen_value).transpose(), modulus );
		result = ModpairwiseSubtractResult( base, result, modulus );
		return result;
	}

	public static void ModLU_Dolittle_Decomposition( Matrix m, float modulus ) { 
		if(modulus == Matrix.ZEROF) {
			m.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int rows = m.rows; int columns = m.columns;
		for(int i = 0; i < rows;i++ ) {
			for(int j = i; j < columns; j++ ) {
				for(int k = 0; k < i-1; k++ ) {
					m.entry[i][j] -= m.entry[i][k]*m.entry[k][j];
					m.entry[i][j] %= modulus;
				}
			}
			for(int j = i+1; j < rows; j++ ) {
				for(int k = 0; k < i-1; k++ ) 
					m.entry[j][i] -= m.entry[j][k]*m.entry[k][i];
				m.entry[j][i] /= m.entry[i][i];
				m.entry[j][i] %= modulus;
			}
		}
	}

	public static Matrix ModLU_Dolittle_DecompositionResult( Matrix m, float modulus ) { 
		if(modulus == Matrix.ZEROF) {
			Matrix result = new Matrix(m);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		int rows = m.rows; int columns = m.columns;
		Matrix lowerupper = new Matrix(Matrix.IDENTITY, rows, columns );
		for(int i = 0; i < rows;i++ ) {
			for(int j = i; j < columns; j++ ) {
				lowerupper.entry[i][j] = m.entry[i][j];
				for(int k = 0; k < i-1; k++ ) {
					lowerupper.entry[i][j] -= lowerupper.entry[i][k]*lowerupper.entry[k][j];
					lowerupper.entry[i][j] %= modulus;
				}
			}
			for(int j = i+1; j < rows; j++ ) {
				lowerupper.entry[j][i] = m.entry[j][i];
				for(int k = 0; k < i-1; k++ ) 
					lowerupper.entry[j][i] -= lowerupper.entry[j][k]*lowerupper.entry[j][i];
				lowerupper.entry[j][i] /= lowerupper.entry[i][i];
				lowerupper.entry[j][i] %= modulus;
			}
		}
		return lowerupper;
	}

	public static Matrix[] ModLU_Dolittle_DecompositionResultPair( Matrix m, float modulus ) { 
		if(modulus == Matrix.ZEROF) {
			Matrix result = new Matrix(m);
			Matrix id = new Matrix(Matrix.IDENTITY, m.rows, m.columns );
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			id.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return new Matrix[] { result, id };
		}
		int rows = m.rows; int columns = m.columns;
		Matrix lower = new Matrix(Matrix.IDENTITY, rows, columns );
		Matrix upper = new Matrix(Matrix.NONE, rows, columns );
		for(int i = 0; i < rows;i++ ) {
			for(int j = i; j < columns; j++ ) {
				upper.entry[i][j] = m.entry[i][j];
				for(int k = 0; k < i-1; k++ ) {
					upper.entry[i][j] -= lower.entry[i][k]*upper.entry[k][j];
					upper.entry[i][j] %= modulus;
				}
			}
			for(int j = i+1; j < rows; j++ ) {
				lower.entry[j][i] = m.entry[j][i];
				for(int k = 0; k < i-1; k++ ) 
					lower.entry[j][i] -= lower.entry[j][k]*upper.entry[j][i];
				lower.entry[j][i] /= upper.entry[i][i];
				lower.entry[j][i] %= modulus;
			}
		}
		return new Matrix[] { lower, upper };
	}

	public static float ModdiagonalDet( Matrix triangular, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			triangular.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return modulus;
		}
		float det = 1.0f;
		int i = triangular.shortest;
		while(--i >= 0 )
			det = (det * triangular.entry[i][i]) % modulus;
		triangular.det = det;
		return det;
	}

	public static float ModDeterminant( Matrix m, float modulus ) { 
		if(modulus == Matrix.ZEROF) {
			m.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return modulus;
		}
		Matrix[] lowerupper = ModLU_Dolittle_DecompositionResultPair(m, modulus);
		float detLower = ModdiagonalDet(lowerupper[0], modulus);
		float detUpper = ModdiagonalDet(lowerupper[1], modulus);
		m.det = (detLower*detUpper) % modulus;
		return m.det;
	}

	public static Matrix ModInverse( Matrix a, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			Matrix result = new Matrix(a);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		Matrix m = new Matrix(a);
		Matrix w = new Matrix( Matrix.NONE, m.rows, m.columns );
		Matrix l = new Matrix( Matrix.NONE, m.rows, m.columns );
		Matrix u = new Matrix( Matrix.NONE, m.rows, m.columns );
		for(int k = 0; k < m.shortest; k++ ) {
			l.entry[k][k] = Matrix.IDENTITY;
			for(int i = k+1; i < m.shortest; i++ ) {
				float coeff = m.entry[i][k]/m.entry[k][k];
				l.entry[i][k] = coeff;
				for(int j = k+1; j < m.shortest; j++ ) {
					m.entry[i][j] -= coeff*m.entry[k][j];
					m.entry[i][j] %= modulus;
				}
			}
		}
		for(int j = 0; j < m.shortest; j++ ) {
			for(int i = 0; i < j; i++ ) {
				u.entry[i][j] = m.entry[i][j] % modulus;
			}
		}
		float b[] = new float[m.shortest];
		float d[] = new float[m.shortest];
		float x[] = new float[m.shortest];
		for(int k = 0; k < m.columns; k++ ) {
			d[1] = b[1];
			b[k] = 1.0f;
			for(int i = 1; i < m.shortest; i++ ) {
				for(int j = 1; j < i-1; j++ ) {
					d[i] = b[i];
					d[i] -= l.entry[i][j]*d[j];
					d[i] %= modulus;
				}
			}
			x[m.shortest-1] = d[m.shortest-1]/u.entry[m.shortest-1][m.shortest-1];
			for(int i = m.shortest; i >= 0; i-- ) {
				x[i] = d[i];
				for(int j = m.shortest; j > i; j-- ) {
					x[i] = x[i] - u.entry[i][j]*x[j];
					x[i] %= modulus;
				}
				x[i] /= u.entry[i][i]; 
			}
			for(int i = 0; i < m.shortest; i++ ) {
				w.entry[i][k] = x[i];
				w.entry[i][k] %= modulus;
			}
			b[k] = 0.0f;
		}
		return w;
	}

	// Thresholding

	public static void Modthreshold ( Matrix a, float threshold, int directed, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			a.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int shortestRow = a.rows;
		int shortestColumn = a.columns ;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float aij = a.entry[i][j];
				if(directed==0)
					aij = Math.abs(a.entry[i][j]);
				if((aij % modulus ) < threshold)
					a.entry[i][j] = 0.0f;
				else
					a.entry[i][j] = 1.0f;
			}  
		}
		a.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix ModthresholdResult ( Matrix a, float threshold, int directed, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			Matrix result = new Matrix(a);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		int shortestRow = a.rows;
		int shortestColumn = a.columns;
		Matrix result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float aij = a.entry[i][j];
				if(directed==0)
					aij = Math.abs(a.entry[i][j]);
				if((aij % modulus) < threshold)
					result.entry[i][j] = 0.0f;
				else
					result.entry[i][j] = 1.0f;
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	// Pair/Entry-Wise Fuzzy Binary arithmetic
	// a AND 0 = 0
	// a or 0 = a
	// a XOR 0 = a
	// a XOR a = 0
	// a AND b = a*b
	// a OR b = a+b
	// a XOR b = (a/b + b/a)
	// a KOR b = -ab + b/a 
	// NOT a = 0
	// a NAND B = a/b
	// a NOR b = a + 1/b

	public static void ModAND( Matrix a, Matrix b, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			a.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int shortestRow = a.rows < b.rows ? a.rows : b.rows;
		int shortestColumn = a.columns < b.columns ? a.columns : b.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				a.entry[i][j] = a.entry[i][j]*b.entry[i][j];
				a.entry[i][j] %= modulus;
			}  
		}
		a.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix ModANDResult ( Matrix a, Matrix b, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			Matrix result = new Matrix(a);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		int shortestRow = a.rows < b.rows ? a.rows : b.rows;
		int shortestColumn = a.columns < b.columns ? a.columns : b.columns;
		Matrix result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				result.entry[i][j] = a.entry[i][j]*b.entry[i][j];
				result.entry[i][j] %= modulus;
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static void ModOR( Matrix a, Matrix b, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			a.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int shortestRow = a.rows < b.rows ? a.rows : b.rows;
		int shortestColumn = a.columns < b.columns ? a.columns : b.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				a.entry[i][j] = (a.entry[i][j]+b.entry[i][j]) % modulus;
			}  
		}
		a.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix ModORResult ( Matrix a, Matrix b, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			Matrix result = new Matrix(a);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		int shortestRow = a.rows < b.rows ? a.rows : b.rows;
		int shortestColumn = a.columns < b.columns ? a.columns : b.columns;
		Matrix result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				result.entry[i][j] = a.entry[i][j]+b.entry[i][j];
				result.entry[i][j] %= modulus;
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static void ModXOR( Matrix a, Matrix b, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			a.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int shortestRow = a.rows < b.rows ? a.rows : b.rows;
		int shortestColumn = a.columns < b.columns ? a.columns : b.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float aij = a.entry[i][j];
				float bij = b.entry[i][j];
				if(aij == 0) {
					if(bij == 0 )
						a.entry[i][j] = 0.0f;
					else
						a.entry[i][j] = bij;
				}
				else {
					if(bij != 0 )
						a.entry[i][j] = aij/bij + bij/aij;
				}
				a.entry[i][j] %= modulus;
			}  
		}
		a.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix ModXORResult ( Matrix a, Matrix b, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			Matrix result = new Matrix(a);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		int shortestRow = a.rows < b.rows ? a.rows : b.rows;
		int shortestColumn = a.columns < b.columns ? a.columns : b.columns;
		Matrix result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float aij = a.entry[i][j];
				float bij = b.entry[i][j];
				if(aij == 0) {
					if(bij == 0 )
						result.entry[i][j] = 0.0f;
					else
						result.entry[i][j] = bij;
				}
				else {
					if(bij != 0 )
						result.entry[i][j] = aij/bij + bij/aij;
				}
				result.entry[i][j] %= modulus;
			}  
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static void ModKOR( Matrix a, Matrix b, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			a.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int shortestRow = a.rows < b.rows ? a.rows : b.rows;
		int shortestColumn = a.columns < b.columns ? a.columns : b.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float aij = a.entry[i][j];
				float bij = b.entry[i][j];
				if(aij != 0) {
					a.entry[i][j] = -aij*bij + bij/aij;
					a.entry[i][j] %= modulus;
				}
			}
		}
		a.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix ModKORResult ( Matrix a, Matrix b, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			Matrix result = new Matrix(a);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		int shortestRow = a.rows < b.rows ? a.rows : b.rows;
		int shortestColumn = a.columns < b.columns ? a.columns : b.columns;
		Matrix result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float aij = a.entry[i][j];
				float bij = b.entry[i][j];
				if(aij != 0) {
					result.entry[i][j] = -aij*bij + bij/aij;
					result.entry[i][j] %= modulus;
				}
			}
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

	public static void ModNOT ( Matrix a, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			a.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return;
		}
		int shortestRow = a.rows;
		int shortestColumn = a.columns;
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float aij = a.entry[i][j];
				if(aij != 0) {
					a.entry[i][j] = 1/-aij;
					a.entry[i][j] %= modulus;
				}
			}
		}
		a.updates &= (Matrix.all_changed^Matrix.entry_changed);
	}

	public static Matrix ModNOTResult ( Matrix a, float modulus ) {
		if(modulus == Matrix.ZEROF) {
			Matrix result = new Matrix(a);
			result.updates &= (Matrix.all_changed^Matrix.entry_changed);
			return result;
		}
		int shortestRow = a.rows;
		int shortestColumn = a.columns;
		Matrix result = new Matrix(Matrix.NONE, shortestRow, shortestColumn);
		int i = shortestRow;
		while(--i >= 0) {
			int j = shortestColumn;
			while(--j >= 0) {
				float aij = a.entry[i][j];
				if(aij != 0) {
					result.entry[i][j] = 1.0f/-aij;
					result.entry[i][j] %= modulus;
				}
			}
		}
		result.updates &= (Matrix.all_changed^Matrix.entry_changed);
		return result;
	}

}
import processing.core.*; 

import java.util.ArrayList; 

public class AspektWorld extends PApplet {
	public static final int BIG_STROKE = 10;
	public static final int MIDDLE_STROKE = 3;
	public static final int SMALL_STROKE = 3;
	// private static final long serialVersionUID = 61803L;
	public final static int FORWARD = 1;
	public final static int BACKWARD = 2;
	public final static int RIGHTWARD = 4;
	public final static int LEFTWARD = 8;
	public final static int UPWARD = 16;
	public final static int DOWNWARD = 32;
	static public void main(String args[]) {
		PApplet.main(new String[] { "--bgcolor=#D4D0C8", "AspektWorld" });
	}
	public int direction = 0;
	public int last_direction = 0;
	public float lmx,lmy,dx,dy;
	public final static float MAX_SPEED = 20.0f;
	public final static float MIN_SPEED = 1.0f;
	public float speed = 20.0f;
	public float wd = 1600;
	public float hh = 900;
	public float dp = 1000;
	public PVector forward_heading;
	public PVector back_heading,left_heading,right_heading,up_heading,down_heading;
	public PVector straight_up, straight_left, straight_down, straight_right, straight_forward, straight_back;
	public int appletMode=1;
	public int diameter = 100;
	public float scalar;
	public int longest;
	public double ttheta;
	// public Rotation r;
	public int double_click;
	public float sbxs,sbys;
	public ArrayList<Line> lines;
	public int edgeColor = 0;
	public int mouseLookOn = 0;
	public int verticesOn = 1;
	public int edgesOn = 1;


	public Graph.Vertex hoverNode;
	public PMatrix3D currCameraMatrix;

	public PGraphics3D g3;
	public int in_a_control = 0;
	public int basics_panel_active = 0;
	public int matrix_panel_active = 0;
	public int ev_panel_active = 0;
	public int instructions_panel_active = 0;
	public int ev1_sz, ev2_sz, ev3_sz, ev4_sz, ev5_sz;
	public int ev_sz = 138, ev_sz_actv = 150;
	public int mpa_sz, mpl_sz, mpd_sz;
	public int m_sz = 100, m_sz_actv = 168;
	public int mpa_sd, mpl_sd, mpd_sd;
	public float m_actv_delta = (m_sz_actv-m_sz)/2.0f;
	public int counter = 0;
	public PFont face; 
	public float th_constant;
	public Visionary cs;

	PVector pos;


	public Graph g1;

	public void attach_graph_to_local_coordinate_master() {
		int i = g1.v.length;
		while(--i>=0)
			g1.v[i].posa = pos;
		i = g1.e.length;
	}

	/**
  public void check_nodes_press() {
    boolean result = false;
    float cz = (float)cam.getDistance();
    float[] la = cam.getLookAt();
    float[] cam_rots = cam.getRotations();
    float mx = mouseX-width/2.0f;
    float my = mouseY-height/2.0f;
    peasy.org.apache.commons.math.geometry.PVector cl = new peasy.org.apache.commons.math.geometry.PVector(la[0],la[1],la[2]);
    Rotation inv_cam = new Rotation(RotationOrder.ZYX, -cam_rots[2], -cam_rots[1], -cam_rots[0]);
    int ns = g1.nodes.size();
    float ndx,ndy,ndz,dx,dy,shifts;
    Node lh = null;
    float lhndz = -height;

    for(int i = 0; i < ns;i++) {
      Node nd = (Node)g1.nodes.get(i);
      if(nd == null)
        continue;
      peasy.org.apache.commons.math.geometry.PVector ndp = new peasy.org.apache.commons.math.geometry.PVector(nd.position.x,nd.position.y,nd.position.z);
      ndp = inv_cam.applyTo(ndp.subtract(cl));
      ndx = (float)(((ndp.x))); ndy = (float)(((ndp.y))); ndz = (float)(((ndp.z)));
      //applet
      //float shifts = (float)((height-64)/(cz-ndp.z));
      shifts = (float)((height-100)/(cz-ndz));
      dx = (float)mx - ndx*shifts;
      dy = (float)my - ndy*shifts;
      if(g1.hoverOnlyOn == 0)
        result = shifts*shifts*nd.radiusq > dx*dx + dy*dy;
      else
        result = 0.5*shifts*nd.radiusq > dx*dx + dy*dy;
      if(result) {
        if(ndz > lhndz && (cz-ndz)>40) {
          lh = nd;  
          lhndz = ndz;
        }
      }
    }
  if(lh != null) {
      Node nd = lh;
      counter++;
      if(nd.partial_display >= 0 && g1.hideOn == 0) {
        //g1.load_more(nd);
        nd.show_hidden();
      }
      else if(g1.hideOn == 1) {
        nd.hide();  
      }
      else if(nd.selected == 0) {
        if(g1.allPaths == 1 || (g1.pointToPointOn == 1 && g1.directedOn == 0))
          g1.select_all_paths(nd,2);
        else if(g1.selectNode != null && g1.pointToPointOn == 1 && g1.directedOn == 1){
          g1.silence_all_paths(2);
          g1.select_all_paths(g1.selectNode, nd,2);
        }
        g1.selectNode = nd;
        nd.selected = 1;
        nd.counter = counter;
        spread_activity(nd);
      }
      else {
        if(g1.allPaths == 1 || g1.pointToPointOn == 1)
          g1.silence_all_paths(2);
        g1.selectNode = null;
        nd.hovered = 0;
        nd.selected = 0;
        nd.counter = 0;
        retract_activity(nd);
      }
    }
    else if(mouseEvent.getClickCount()==2) {
      silence_grid();
      g1.silence_all_paths(2);
    }
    else if(mouseEvent.getClickCount()==3) {
      cam.reset();
    } 
  }
	 **/
	public void check_matrices_pass() {
		float mx = mouseX;
		float my = mouseY;
		float m1x =  (longest == 1)? 53 : 53;
		float m1y =  (longest==1)? height-253 : 53;
		float m2x = (longest == 1)? width/2 : 53;
		float m2y = (longest==1)? height-253 : height/2-100;
		float m3x = (longest == 1)? width-253 : 53;
		float m3y = (longest==1)? height-253 : height-253;
		if( mx < m1x + mpa_sz && mx > m1x && my < m1y + mpa_sz && my > m1y ) {
			if(mpa_sz == m_sz) 
				mpa_sz = m_sz_actv;
		}
		else if(mpa_sd == 0)
			mpa_sz = m_sz;
		if( mx < m2x + mpl_sz && mx > m2x && my < m2y + mpl_sz && my > m2y ) {
			if(mpl_sz == m_sz)
				mpl_sz = m_sz_actv;
		}
		else if(mpl_sd == 0)
			mpl_sz = m_sz;
		if( mx < m3x + mpd_sz && mx > m3x &&  my < m3y + mpd_sz && my > m3y ) {
			if(mpd_sz == m_sz)
				mpd_sz = m_sz_actv;
		}
		else if(mpd_sd == 0)
			mpd_sz = m_sz;
	}

	public void check_matrices_press() {
		float mx = mouseX;
		float my = mouseY;
		float m1x =  (longest == 1)? 53 : 53;
		float m1y =  (longest==1)? height-253 : 53;
		float m2x = (longest == 1)? width/2 : 53;
		float m2y = (longest==1)? height-253 : height/2-100;
		float m3x = (longest == 1)? width-253 : 53;
		float m3y = (longest==1)? height-253 : height-253;
		if( mx < m1x + mpa_sz && mx > m1x ) {
			if( my < m1y + mpa_sz && my > m1y ) {
				if(mpa_sd == 0) {
					mpa_sz = m_sz_actv;
					mpa_sd = 1;
				} 
				else {
					mpa_sz = m_sz;
					mpa_sd = 0;
				}
				return; 
			}  
		}
		if( mx < m2x + mpl_sz && mx > m2x ) {
			if( my < m2y + mpl_sz && my > m2y ) {
				if(mpl_sd == 0) {
					mpl_sz = m_sz_actv;
					mpl_sd = 1;
				} 
				else {
					mpl_sz = m_sz;
					mpl_sd = 0;
				}
				return;    
			}  
		}
		if( mx < m3x + mpd_sz && mx > m3x ) {
			if( my < m3y + mpd_sz && my > m3y ) {
				if(mpd_sd == 0) {
					mpd_sz = m_sz_actv;
					mpd_sd = 1;
				} 
				else {
					mpd_sz = m_sz;
					mpd_sd = 0;
				}
				return; 
			}  
		}
	}


	public void check_nodes_pass() {
		boolean result = false;
		float cz = pos.mag();
		float mx = mouseX-width/2.0f;
		float my = mouseY-height/2.0f;
		PVector cl = new PVector(0,0,0);
		cl.set(forward_heading);
		//inv_cam.invert();
		int ns = g1.v.length;
		float ndx,ndy,ndz,dx,dy,shifts;
		PMatrix3D inv_cam = new PMatrix3D(g3.camera);
		if(hoverNode != null) {
			hoverNode.hover = 0;
			hoverNode.radius/=1.618;
			hoverNode = null;
		}
		PVector ndt = new PVector(), ndp = new PVector();
		for(int i = 0; i < ns;i++) {
			Graph.Vertex nd = g1.v[i];
			ndp.set(nd.pos);
			ndp.add(pos);
			ndp.sub(cl);
			cz = ndp.mag();
			//ndt.set(ndp);
			inv_cam.mult(ndp, ndt);
			ndx = (((ndt.x))); ndy = (((ndt.y))); ndz = (((ndt.z)));
			//applet
			//float shifts = (float)((height-64)/(cz-ndp.getZ()));
			shifts = ((2*height-112)/(cz-ndz));
			dx = mx - ndx*shifts;
			dy = my - ndy*shifts;
			/**
                        currCameraMatrix = new PMatrix3D(g3.camera);
        		camera();
        	        fill(nd.color);
                        rect(ndx*shifts+width/2,ndy*shifts+height/2,nd.radius*shifts,nd.radius*shifts);
                  	g3.camera = currCameraMatrix;
  		        //popMatrix();
			 **/
			result = shifts*nd.radius*nd.radius > dx*dx + dy*dy;
			if(result) {
				if(hoverNode != null) {
					hoverNode.radius /= 1.618;
					hoverNode.hover = 0;
					hoverNode = null;
				}
				hoverNode = nd;
				nd.hover = 1;
				nd.radius *= 1.618;
			}
		}
	}

	@Override
	public void draw() {
		in_a_control = 0;
		background(0);
		draw_grid();
		g1.draw();
		//currCameraMatrix = new PMatrix3D(g3.camera);
		//camera();
		//if(matrix_panel_active != 0)
		//	draw_matrix_panel();
		//draw_controls();
		//g3.camera = currCameraMatrix;
		update_heading();
	}

	public void draw_ev(int index, ArrayList y, float xc, float yc, int sz ) {
		rectMode(CORNER);
		strokeWeight(1);
		stroke(360,255,250,255);
		noFill();
		rect(xc-sz-3,yc-105,sz,105);
		plot_v(y, xc+3, yc, sz, -1, true);
	}

	public void draw_grid() {

		int linum = lines.size();
		for(int i = 0; i < linum;i++) {
			Line lin = lines.get(i);
			lin.draw();
		}  
	}

	public void draw_instructions() {
		float th = 1.618f*th_constant;
		float xo = width/8;
		float yo = height/4;
		textSize(20);
		fill(0,0,50,150);
		rectMode(CORNER);
		rect(xo,yo,screenWidth*2/3,22*th);
		rectMode(CENTER);
		fill(120,250,250,230);
		text("k - toggle directed",xo+33,yo+th);  
		text(", - layout",xo+33,yo+2*th);  
		text("b - re-box",xo+33,yo+3*th);  
		text("f - display an xml or xml map file",xo+33,yo+4*th);  
		text("j - highlight connected component (use with d)",xo+33,yo+5*th);  
		text("p - highlight all paths. To use: select a vertex, toggle p, hover or select another vertex",xo+33,yo+6*th);  
		text("esc - quit",xo+33,yo+7*th);  
		text("r - new random graph",xo+33,yo+8*th);
		text("mouse : left-drag to rotate, centre-button-drag to translate, right-drag to zoom",xo+33,yo+9*th);
		text("mouse : click to select/unselect, double click to unselect all",xo+33,yo+10*th);
		text("mouse : triple click to reset model.",xo+33,yo+11*th);
		text("h - toggle Focus mode",xo+33,yo+12*th);  
		text("e - toggle edges display",xo+33,yo+13*th);  
		text("v - toggle vertices display",xo+33,yo+14*th);  
		text("i - toggle hide mode (click on a vertext to hide it)",xo+33,yo+15*th);
		text("    when off, clicking on partially displayed vertices unhides their neighbours",xo+33,yo+16*th);
		text("left menu panel - control of matrices chosen for layout : adjacency, laplacian, shortest distance",xo+33,yo+17*th);
		text("bottom menu panel - control of edge length, vertext size",xo+33,yo+18*th);
		text("                    toggle size norming (normed by default), color range",xo+33,yo+19*th);
		text("                    toggle color norming (normed by default)",xo+33,yo+20*th);
		text("right menu panel - view of eigen values (left), and first 5 eigenvectors",xo+33,yo+21*th);
	}


	public void draw_matrix_panel() {
		float m1x =  (longest == 1)? 53 : 53;
		float m1y =  (longest==1)? height-253 : 53;
		show_matrix_f( 0, m1x, m1y, mpa_sz ); 
		float m2x = (longest == 1)? width/2 : 53;
		float m2y = (longest==1)? height-253 : height/2-100;
		show_matrix_f( 1, m2x, m2y, mpl_sz );
		float m3x = (longest == 1)? width-253 : 53;
		float m3y = (longest==1)? height-253 : height-253;
		show_matrix_f( 2, m3x, m3y, mpd_sz );
	}





	public void dvals() {
		pos = new PVector(0,0,0);
		lmx = 0.0f; lmy = 0.0f;
		//int boxnum = boxes.size();
		forward_heading = new PVector(250,10,250);
		straight_up = new PVector(0,1,0);
		straight_down = new PVector(0,-1,0);
		straight_left = new PVector(-1,0,0);
		straight_right = new PVector(1,0,0);
		straight_forward = new PVector(0,0,1);
		straight_back = new PVector(0,0,-1);
		make_normals();
		//cam.lookAt(forward_heading.x, forward_heading.y, forward_heading.z, 1.0f);

		camera(0,0,0, forward_heading.x, forward_heading.y, forward_heading.z, 0.0f, 1.0f, 0.0f);
	}


	public void input() {
		float mx = mouseX;
		float my = mouseY;
		//if(result < height*height/4 ) {
		if(mouseLookOn==1)
			cursor(CROSS);  
		else
			cursor(HAND);
		mx = map(mx, 0, width, -1, 1);
		my = map(my, 0, height, -1, 1); 
		//}
		//else {
		// cursor(HAND);
		// return; 
		//}
		float xang = mx * Math.abs(mx) / (10*PI);
		float yang = my * Math.abs(my) / (12*PI);
		if(mouseLookOn==1) {
			right_heading.mult(speed*sin(xang));
			forward_heading.mult(cos(xang));
			forward_heading.add(right_heading);
			straight_up.mult(speed*sin(yang));
			forward_heading.mult(cos(yang));
			forward_heading.add(straight_up);
		}
		//cam.lookAt(forward_heading.x, forward_heading.y, forward_heading.z, 1.0f);
		camera(0,0,0, forward_heading.x, forward_heading.y, forward_heading.z, 0.0f, 1.0f, 0.0f);
		forward_heading.normalize();
		update_normals();
		forward_heading.mult(speed);
	}

	/**
  public void draw_touchpads ( ) {
    float cz = (float)cam.getDistance();
    float[] pos = cam.getPosition();
    float[] la = cam.getLookAt();
    float[] cam_rots = cam.getRotations();
    peasy.org.apache.commons.math.geometry.PVector cl = new peasy.org.apache.commons.math.geometry.PVector(la[0],la[1],la[2]);
    peasy.org.apache.commons.math.geometry.PVector cp = new peasy.org.apache.commons.math.geometry.PVector(-pos[0],-pos[1],-pos[2]);
    Rotation inv_cam = new Rotation(RotationOrder.ZYX, -cam_rots[2], -cam_rots[1], -cam_rots[0]);
    int ns = g1.nodes.size();
    for(int i = 0; i < ns;i++) {
      Node nd = (Node)g1.nodes.get(i);
      if(nd == null)
        continue;
      peasy.org.apache.commons.math.geometry.PVector ndp = new peasy.org.apache.commons.math.geometry.PVector(nd.position.x,nd.position.y,nd.position.z);
      ndp = inv_cam.applyTo(ndp.subtract(cl));
      fill(nd.selectcolor);
      float ndx = (float)(((ndp.x))); 
      float ndy = (float)(((ndp.y)));
      float ndz = (float)(((ndp.z)));
      // for applet
      //float shifts = 1.0f;
      //if(online)
      //  float shifts = (float)((height-64)/(cz-ndp.z));
      //else
      float shifts = (float)((height-100)/(cz-ndp.z));
      rect((float)(ndx*shifts+width/2.0),(float)(ndy*shifts+height/2.0),(float)nd.radius*shifts, (float)nd.radius*shifts);
      text(nd.num+"",(float)(ndx*shifts+width/2.0),(float)ndy*shifts+height/2.0-nd.radius*shifts/1.618);
    }  
  }
	 **/

	/**
  public void draw_texts ( ) {
    float cz = (float)cam.getDistance();
    float[] pos = cam.getPosition();
    float[] la = cam.getLookAt();
    float[] cam_rots = cam.getRotations();
    peasy.org.apache.commons.math.geometry.PVector cl = new peasy.org.apache.commons.math.geometry.PVector(la[0],la[1],la[2]);
    peasy.org.apache.commons.math.geometry.PVector cp = new peasy.org.apache.commons.math.geometry.PVector(-pos[0],-pos[1],-pos[2]);
    Rotation inv_cam = new Rotation(RotationOrder.ZYX, -cam_rots[2], -cam_rots[1], -cam_rots[0]);
    int ns = g1.nodes.size();
    for(int i = 0; i < ns;i++) {
      Node nd = (Node)g1.nodes.get(i);
      if(nd == null)
        continue;
      if(nd.selected + nd.hovered == 0)
        continue;
      if(nd.label.length() + nd.special.length() == 0)
        continue;
      peasy.org.apache.commons.math.geometry.PVector ndp = new peasy.org.apache.commons.math.geometry.PVector(nd.position.x,nd.position.y,nd.position.z);
      ndp = inv_cam.applyTo(ndp.subtract(cl));
      fill(nd.selectcolor);
      float ndx = (float)(((ndp.x))); 
      float ndy = (float)(((ndp.y)));
      float ndz = (float)(((ndp.z)));
      // for applet
      //float shifts = 1.0f;
      //if(online)
      //  float shifts = (float)((height-64)/(cz-ndp.z));
      //else
      float shifts = (float)((height-100)/(cz-ndp.z));
      strokeWeight(1);
      stroke(color(0,0,255,255));
      textSize(20);
      float tw = textWidth(nd.label) <textWidth(nd.special) ? textWidth(nd.special)*1.2f : textWidth(nd.label)*1.2f;
      float th = textAscent();
      //textSize(10);
      fill(0,0,50,255);
      rectMode(CORNER);
      float tsx = tw*0.309;
      float tsh = th*0.618;
      rect((float)(ndx*shifts+width/2.0f)+sbxs,(float)ndy*shifts+height/2.0f-th+sbys, tw*1.618, 3.236*th);
      rectMode(CENTER);
      fill(color(0,0,255,200));
      text(nd.label,(float)(ndx*shifts+width/2.0f)+tsx+sbxs,(float)ndy*shifts+height/2.0f+sbys+tsh);
      textSize(14);
      text(nd.special,(float)(ndx*shifts+width/2.0f)+tsx+sbxs,(float)ndy*shifts+height/2.0f+sbys+tsh+th);
      if(nd.partial_display >= 0) {
        if(nd.storage_type == nd.MEM_STORAGE)
          text(" +" + nd.undirectededges.size() );
        else if(nd.storage_type == nd.LEGACY_STORAGE)
          text(" +" + nd.serial.getChildCount() );
      }
    }  
  }
	 **/
	@Override
	public void keyPressed(){
		if(key != CODED)
			switch(key) {
			case('w'):case('W'):direction |=FORWARD; last_direction = FORWARD; break;
			case('d'):case('D'):direction |=RIGHTWARD; last_direction = RIGHTWARD; break;
			case('s'):case('S'):direction |=BACKWARD; last_direction = BACKWARD; break;
			case('a'):case('A'):direction |=LEFTWARD; last_direction = LEFTWARD; break;
			case(' '):direction |=UPWARD; last_direction = UPWARD; break;
			case('c'):case('C'):direction |=DOWNWARD; last_direction = DOWNWARD; break;
			case('r'):case('R'):
				cs.reEmbedWithSameFilter();
			break;
			case('1'):
				cs.basicMatrixEmbed();
			break;
			case('2'):
				cs.basicMatrixProjectedEmbed();
			break;
			case('3'):
				cs.eigeneeringEmbed();
			break;
			case('4'):
				cs.randomFilterEmbed();
			break;
			case('['):
				cs.offset = cs.offset - 1;
			if(cs.offset < 0)
				cs.offset = 0;
			//println(cs.offset);
			cs.reEmbedWithSameFilter();
			break;
			case(']'):
				cs.offset = cs.offset + 1;
			if(cs.offset > cs.g1.v.length-5)
				cs.offset = 0;
			//println(cs.offset);
			cs.reEmbedWithSameFilter();
			break;
			case('`'):case('~'):
				mouseLookOn ^= 1;
			break;
			case('v'):case('V'):
				verticesOn ^= 1;
			break;
			case('b'):case('B'):
				cs.boxit();
			break;
			case('e'):case('E'):
				edgesOn ^= 1;
			break;
			}
		else {
			switch(keyCode) {
			case(SHIFT):
				if(speed > MIN_SPEED) {
					speed = MIN_SPEED; 
					break;
				}
			}
		}
	}

	@Override
	public void keyReleased(){  
		if(key != CODED)
			switch(key) {
			case('w'):case('W'):direction ^=FORWARD; last_direction = FORWARD; break;
			case('d'):case('D'):direction ^=RIGHTWARD; last_direction = RIGHTWARD; break;
			case('s'):case('S'):direction ^=BACKWARD; last_direction = BACKWARD; break;
			case('a'):case('A'):direction ^=LEFTWARD; last_direction = LEFTWARD; break;
			case(' '):direction ^=UPWARD; last_direction = UPWARD; break;
			case('c'):case('C'):direction ^=DOWNWARD; last_direction = DOWNWARD; break;

			}  
		else
			switch(keyCode) {
			case(SHIFT):
				if(speed < MAX_SPEED) {
					speed = MAX_SPEED; 
					break;
				}
			}
	}

	public void make_grid() {
		float spacing = 100;
		lines = new ArrayList((int)(2*(wd+dp)/spacing + 1));
		stroke(200,255,255,255);
		for(float i = -wd*2; i <= wd*2; i+=spacing) {
			Line nl = new Line(this, i,0,-2*wd,i,0,2*wd);
			nl.posa = pos;
			lines.add(nl);
		}
		for(float i = -wd*2; i <= wd*2; i+= spacing) {
			Line nl = new Line(this, -2*wd,0,i,2*wd,0,i);
			nl.posa = pos;
			lines.add(nl);
		}  
		stroke(0,0,0,255);
	}

	public void make_normals() {
		back_heading = new PVector(-forward_heading.x, -forward_heading.y, -forward_heading.z);
		up_heading = new PVector(0,speed,0);
		down_heading = new PVector(0,-speed,0);
		right_heading = forward_heading.cross(up_heading);
		left_heading = forward_heading.cross(down_heading);
	}



	@Override
	public void mouseClicked() {
		if(mouseEvent.getClickCount() == 3) { 
			dvals();
		}
	}
	@Override
	public void mouseDragged(){}

	@Override
	public void mouseMoved() { 
		if(matrix_panel_active==1)
			check_matrices_pass();
		//check_nodes_pass();
	}

	@Override
	public void mousePressed() {
		if(matrix_panel_active==1)
			check_matrices_press();
		//check_nodes_press(); 
	}

	@Override
	public void mouseReleased() { }

	public void plot_v(ArrayList y, float xc, float yc, int sz, int highlight, boolean anchored) {
		int ys = y.size();
		int xstep = sz/ys;
		if(xstep < 1)
			xstep = 1;
		float  lx = xc;
		float  ly = yc;
		if(anchored) {
			strokeWeight(1);
			stroke(360,255,250,255);
			line(lx,ly,lx+ys*xstep,ly);
		}
		for(int i = 0; i < ys;i++) {
			float vi = (Float)y.get(i);
			float xi = xc+i*xstep;
			float yi = yc-vi;
			strokeWeight(1);
			stroke(255,255,255,255);    
			if(!anchored) {
				if(i > 0)
					line(lx,ly,xi,yi);
			}
			else
				line(lx,ly,xi,yi);
			if(i == highlight)
				rect(xi,yi,7,7);
			lx = xi;
			ly = yi;
		}
		if(anchored) {
			line(lx,ly,xc+ys*xstep,yc);
		}
	}

	@Override
	public void setup() {
		// For applet
		face = createFont("MingLiU",40);
		textSize(20);
		th_constant = textAscent();
		textFont(face);
		int xw = (screenWidth*15)/16;
		int yh = (screenHeight*7)/8;
		//applet
		//m_sz = 100;
		//m_sz_actv = 110;
		//m_actv_delta = (m_sz_actv-m_sz)/2.0f;
		//if(!online)
		//size(xw,yh, P3D);
		//else
		size(900,568,P3D);
		sbxs = width/50;
		sbys = -height/50;
		rectMode(CENTER);
		colorMode(HSB,360,255,255,255);
		longest = width < height ? 1 : 0;
		scalar = width < height ? width : height;
		//tthea = tan(scalar);
		g3 = (PGraphics3D)g;
		background(0);
		cs = new Visionary(this);
		cs.randomGraph();
		cs.randomizeLayoutColorSize();
		g1 = cs.g1;
		dvals();
		attach_graph_to_local_coordinate_master();
		make_normals();
		make_grid();
		update_heading();
		// not for applet
	}

	public void setup_matrix_panel() {
		mpa_sz = m_sz;
		mpa_sd = 0;
		mpl_sz = m_sz;
		mpl_sd = 0;
		mpd_sz = m_sz_actv;
		mpd_sd = 1;
	}

	public void show_matrix_f( int index, float xc, float yc, int sz ) {
		PGraphics m=null; // = (PGraphics)g1.matrix_cache.get(index);
		if(m != null) {
			if(sz == m_sz) {
				//			image((PImage)g1.matrix_cache.get(index), xc, yc, sz, sz);
			}
			else if(sz == m_sz_actv) {
				//			image((PImage)g1.matrix_cache.get(index), xc, yc, sz, sz);
			}
		}
		else {
			rectMode(CORNER);
			fill(0,0,128,128);
			rect(xc+m_sz/5,yc+m_sz/5,m_sz/2,m_sz/2);
			rectMode(CENTER);
		}
	}

	public void silence_grid( ) {

	}

	public void update_heading() {
		input();
		check_nodes_pass();
		update_v();
	}
	public void update_normals() {
		back_heading.set(-forward_heading.x, -forward_heading.y, -forward_heading.z); 
		up_heading.set(0,speed,0);
		down_heading.set(0,-speed,0);
		straight_left.set(-1,0,0);
		straight_right.set(1,0,0);
		straight_forward.set(0,0,1);
		straight_back.set(0,0,-1);
		straight_up.set(0,1,0);
		straight_down.set(0,-1,0);
		right_heading = forward_heading.cross(straight_up);
		left_heading = forward_heading.cross(straight_down);
		back_heading.mult(speed);
		right_heading.mult(speed);
		left_heading.mult(speed);
	}

	/**
  public void check_nodes_pass() {
    boolean result = false;
    float cz = (float)cam.getDistance();
    float[] la = cam.getLookAt();
    float[] cam_rots = cam.getRotations();
    float mx = mouseX-width/2.0f;
    float my = mouseY-height/2.0f;
    peasy.org.apache.commons.math.geometry.PVector cl = new peasy.org.apache.commons.math.geometry.PVector(la[0],la[1],la[2]);
    Rotation inv_cam = new Rotation(RotationOrder.ZYX, -cam_rots[2], -cam_rots[1], -cam_rots[0]);
    int ns = g1.nodes.size();
    float ndx,ndy,ndz,dx,dy,shifts;
    Node lh = null;
    float lhndz = -height;
    if(g1.hoverNode != null) {
      if(g1.allPaths == 1 || (g1.pointToPointOn == 1))
        g1.silence_all_paths(1);
      retract_activity(g1.hoverNode);  
      g1.hoverNode.hovered = 0;
      g1.hoverNode = null;
    }
    for(int i = 0; i < ns;i++) {
      Node nd = (Node)g1.nodes.get(i);
      if(nd == null)
        continue;
      peasy.org.apache.commons.math.geometry.PVector ndp = new peasy.org.apache.commons.math.geometry.PVector(nd.position.x,nd.position.y,nd.position.z);
      ndp = inv_cam.applyTo(ndp.subtract(cl));
      ndx = (float)(((ndp.x))); ndy = (float)(((ndp.y))); ndz = (float)(((ndp.z)));
      //applet
      //float shifts = (float)((height-64)/(cz-ndp.z));
      shifts = (float)((height-100)/(cz-ndz));
      dx = (float)mx - ndx*shifts;
      dy = (float)my - ndy*shifts;
      if(g1.hoverOnlyOn == 0)
        result = shifts*shifts*nd.radiusq > dx*dx + dy*dy;
      else
        result = 0.5*shifts*nd.radiusq > dx*dx + dy*dy;
      if(result) {
        if(ndz > lhndz && (cz-ndz) > 40) {
          lh = nd;  
          lhndz = ndz;
        }
      }
    }
    if(lh != null) {
       Node nd = lh;
       if(g1.hoverNode != null) {
          if(g1.allPaths == 1 || g1.pointToPointOn == 1)
            g1.silence_all_paths(1);
          g1.hoverNode.hovered = 0;
          retract_activity(g1.hoverNode);
          g1.hoverNode = null;
        }
        if(g1.allPaths == 1 || (g1.pointToPointOn == 1 && g1.directedOn == 0))
          g1.select_all_paths(nd,1);
        else if(g1.selectNode != null && g1.pointToPointOn == 1 && g1.directedOn == 1) {
          g1.select_all_paths(g1.selectNode, nd,1);
        }
        g1.hoverNode = nd;
        nd.hovered = 1;
        spread_activity(nd);
    }
  }
	 **/

	public void update_selected() {
		int ns = g1.v.length;
		for(int i = 0; i < ns;i++) {
			Graph.Vertex n = g1.v[i];
			if(n == null)
				continue;
		}
	}

	public void update_v() {
		switch(direction) {
		case FORWARD:pos.sub(forward_heading); break;
		case BACKWARD:pos.sub(back_heading); break;
		case LEFTWARD:pos.sub(left_heading); break;
		case RIGHTWARD:pos.sub(right_heading); break;
		case UPWARD:pos.add(up_heading); break;
		case DOWNWARD:pos.add(down_heading); break;
		case FORWARD|LEFTWARD:pos.sub(forward_heading);pos.sub(left_heading); break;
		case FORWARD|RIGHTWARD:pos.sub(forward_heading);pos.sub(right_heading);  break;
		case BACKWARD|LEFTWARD:pos.sub(back_heading);pos.sub(left_heading);  break;
		case BACKWARD|RIGHTWARD:pos.sub(back_heading);pos.sub(right_heading);  break;
		case UPWARD|FORWARD:pos.sub(forward_heading);pos.add(up_heading);  break;
		case UPWARD|BACKWARD:pos.sub(back_heading);pos.add(up_heading);  break;
		case UPWARD|LEFTWARD:pos.add(up_heading);pos.sub(left_heading);  break;
		case UPWARD|RIGHTWARD:pos.add(up_heading);pos.sub(right_heading);  break;
		case DOWNWARD|FORWARD:pos.sub(forward_heading);pos.add(down_heading);  break;
		case DOWNWARD|BACKWARD:pos.sub(back_heading);pos.add(down_heading);  break;
		case DOWNWARD|LEFTWARD:pos.add(down_heading);pos.sub(left_heading);  break;
		case DOWNWARD|RIGHTWARD:pos.add(down_heading);pos.sub(right_heading);  break;
		case UPWARD|FORWARD|RIGHTWARD:pos.sub(forward_heading);pos.add(up_heading);pos.sub(right_heading);break;
		case UPWARD|FORWARD|LEFTWARD:pos.sub(forward_heading);pos.add(up_heading);pos.sub(left_heading);  break;
		case UPWARD|BACKWARD|RIGHTWARD:pos.sub(back_heading);pos.add(up_heading);pos.sub(right_heading);  break;
		case UPWARD|BACKWARD|LEFTWARD:pos.sub(back_heading);pos.add(up_heading);pos.sub(left_heading);  break;
		case DOWNWARD|FORWARD|LEFTWARD:pos.sub(forward_heading);pos.add(down_heading);pos.sub(left_heading);  break;
		case DOWNWARD|FORWARD|RIGHTWARD:pos.sub(forward_heading);pos.add(down_heading);pos.sub(right_heading);  break;
		case DOWNWARD|BACKWARD|LEFTWARD:pos.sub(back_heading);pos.add(down_heading);pos.sub(left_heading); break;
		case DOWNWARD|BACKWARD|RIGHTWARD:pos.sub(back_heading);pos.add(down_heading);pos.sub(right_heading);  break;
		case LEFTWARD|RIGHTWARD: 
			if(last_direction==LEFTWARD)
				pos.sub(left_heading);
			else
				pos.sub(right_heading);  
			break;
		case FORWARD|BACKWARD: 
			if(last_direction==FORWARD)
				pos.sub(forward_heading);
			else
				pos.sub(back_heading);
			break;
		case UPWARD|DOWNWARD:
			if(last_direction==UPWARD)
				pos.add(up_heading);
			else
				pos.add(down_heading);
			break;
		} 
		if(pos.y < height/4)
			pos.y = height/4;
	}
}
import java.util.HashMap;
import java.util.HashSet;
import java.util.Iterator;
import java.util.Set;

import processing.core.PVector;

// The local image of the global concept
public class Graph { // implements Comparable<Graph> .... just joking
	public class MatrixTypeIndex {
		public static final int ADJACENCY = 0;
		public static final int DEGREE = 1;
		public static final int ADJACENCY_2STEP = 2;
		public static final int DEGREE_1STEP = 3;
		public static final int DEGREE_2STEP = 4;
		public static final int SEIDELADJACENCY = 5;
		public static final int AUGMENTEDADJACENCY = 6;
		public static final int SHAREOFDEGREE = 7;
		public static final int MINIMUMDISTANCEMATRIX = 8;
		public static final int GENERALIZEDADJACENCY = 9;
		public static final int WEIGHTEDDEGREE = 10;
		public static final int WEIGHTEDADJACENCY = 11;
		public static final int SPIELMANLAZYWALK = 12;
		public static final int LAZYWALKSOFTINVERSE = 13;
		public static final int NSTEPWALKS = 14;
		public static final int TRANSITIONMATRIXRANDOMWALK = 15;
		public static final int LAPLACIAN = 16;
		public static final int SOFTNORMALIZEDLAPLACIAN = 17;
		public static final int HARDNORMALIZEDLAPLACIAN = 18;
		public static final int COMMUTETIME = 19;
		public static final int DEFORMEDLAPLACIAN = 20;
		public static final int PROJECTIONMATRIX = 21;
		public static final int MAXIMUMDISTANCEMATRIX = 22;
		public static final int MATRIX_FIELD_COUNT = 23;
	}
	public class EdgeBattle { 
		public static final int ESCHER = -8;
		public static final int ONEFROMTHREE = -5;
		public static final int TOTALFAIL = -2;
		public static final int DEPENDSWHICHWAYYOULOOKATIT = 1;
		public static final int PARALLEL = 4;
	}
	public class Vertex {
		public int id; public PVector pos, posa; public int color; public float radius; public int hover;
		public Vertex() {
			id = 0; 
			hover = 0;
			pos = new PVector(0,0,0);
		}
		public void draw() {
			aspektWorld.pushMatrix();
			aspektWorld.fill(color);
			aspektWorld.stroke(color);
			aspektWorld.strokeWeight(aspektWorld.SMALL_STROKE);
			aspektWorld.translate(posa.x + pos.x,posa.y + pos.y, posa.z + pos.z);
			aspektWorld.box(radius);
			aspektWorld.popMatrix();
		}
	}
	public class Edge {
		public final int edgeColorB = aspektWorld.color(141,250,250,250); 
		public Vertex a; public Vertex b; public int[] strength; 
		public Edge( Vertex av, Vertex bv ) {
			a = av; b = bv; strength = new int[6];
		}
		public void draw() {
			if(a.hover==1) {
				aspektWorld.strokeWeight(aspektWorld.BIG_STROKE);
				aspektWorld.stroke(edgeColorB);
			}
			else {
				aspektWorld.strokeWeight(aspektWorld.SMALL_STROKE);
				aspektWorld.stroke(a.color);
			}
			aspektWorld.line(a.posa.x + a.pos.x,a.posa.y + a.pos.y,a.posa.z + a.pos.z,
					a.posa.x + b.pos.x,a.posa.y + b.pos.y,a.posa.z + b.pos.z);
		}
	}
	public HashMap<Vertex, Set<Edge>> edgesByVertexA, edgesByVertexB;
	public HashMap<Vertex,Integer> numbering;
	public HashMap<Edge,Integer> edgeNumbering;
	public Matrix[] matrix;
	public Vertex[] v;
	public Edge[] e;
	public float volume;
	private final AspektWorld aspektWorld;

	public Graph(AspektWorld aspektWorld) {
		this.aspektWorld = aspektWorld;
		edgesByVertexA = new HashMap<Vertex, Set<Edge>>();
		edgesByVertexB = new HashMap<Vertex, Set<Edge>>();
		matrix = new Matrix[MatrixTypeIndex.MATRIX_FIELD_COUNT];
	}

	public HashMap<Graph.Vertex,Integer> numbering() {
		if(numbering != null)
			return numbering;
		int vertices = edgesByVertexA.size();
		Set<Graph.Vertex> verticeSet = edgesByVertexA.keySet();
		Iterator<Graph.Vertex> vertexIterator = verticeSet.iterator();
		numbering = new HashMap<Graph.Vertex,Integer> ();
		int i = vertices;
		while(--i >= 0)
			numbering.put(vertexIterator.next(), new Integer(i));
		return numbering;
	}

	public HashMap<Graph.Edge, Integer> edgeNumbering() {
		if(edgeNumbering != null)
			return edgeNumbering;
		Set<Graph.Vertex> verticeSet = edgesByVertexA.keySet();
		Iterator<Graph.Vertex> vertexIterator = verticeSet.iterator();
		edgeNumbering = new HashMap<Edge,Integer> ();
		int edges = 0;
		while(vertexIterator.hasNext()) {
			Graph.Vertex stupidFuckFace = vertexIterator.next();
			Set<Graph.Edge> fuckingMoron = edgesByVertexA.get(stupidFuckFace);
			Iterator<Graph.Edge> edgeIterator = fuckingMoron.iterator();
			while(edgeIterator.hasNext()) {
				Edge e = edgeIterator.next();
				if(!edgeNumbering.containsKey(e))
					edgeNumbering.put(e,edges++);
			}
		}
		return edgeNumbering;
	}

	public void addVertex( Vertex  a ) {
		if(edgesByVertexA.containsKey(a))
			return;
		edgesByVertexA.put( a, new HashSet<Edge>() );
		edgesByVertexB.put( a, new HashSet<Edge>() );
	}

	public void cutVertex( Vertex a ) {
		edgesByVertexA.remove( a );
		edgesByVertexB.remove( a );
	}

	public void addEdge( Edge e ) {
		edgesByVertexA.get(e.a).add(e);
		edgesByVertexB.get(e.b).add(e);
	}

	public void cutEdge( Edge e ) {
		edgesByVertexA.get(e.a).remove(e);
		edgesByVertexB.get(e.b).remove(e);
	}

	public float degree ( Vertex v ) {
		int index = numbering.get(v);
		return matrix[Graph.MatrixTypeIndex.DEGREE].entry[index][index];  
	}

	public void crystallizeGraph() {
		numbering();
		edgeNumbering();
		Set<Vertex> vertices = numbering.keySet();
		v = new Vertex[vertices.size()];
		vertices.toArray(v);
		Set<Edge> edges = edgeNumbering.keySet();
		e = new Edge[edges.size()];
		edges.toArray(e);
	}

	public void draw() {
		// TODO Auto-generated method stub
		int i = e.length;
		if(aspektWorld.edgesOn==1)
			while(--i >= 0)
				e[i].draw();
		else 
			while(--i >= 0)
				if(e[i].a.hover==1)
					e[i].draw();

		i = v.length;
		if(aspektWorld.verticesOn==1)
			while(--i >= 0)
				v[i].draw();
	}

}
import processing.core.PVector;


class Line {
	/**
	 * 
	 */
	private final AspektWorld aspektWorld;
	public PVector another, posa;
	public PVector yetanother;
	public int lc;

	public Line( AspektWorld aspektWorld, float ax, float ay, float az, float tx, float ty, float tz ) {
		this.aspektWorld = aspektWorld;
		another = new PVector(ax,ay,az);
		yetanother = new PVector(tx,ty,tz);
		lc = this.aspektWorld.color (180,200,200,200);
	}  

	public void draw() {
		this.aspektWorld.stroke(lc);
		this.aspektWorld.strokeWeight(aspektWorld.MIDDLE_STROKE);
		this.aspektWorld.line(posa.x+another.x,posa.y+another.y,posa.z+another.z,posa.x+yetanother.x,posa.y+yetanother.y,posa.z+yetanother.z);  
	}
}
import java.util.Arrays;

public class Matrix {
	public int rows, columns, shortest, symmetric;
	public int i, j;
	public static final int NONE = -1;
	public static final int ZERO = 0;
	public static final int IDENTITY = 1;
	public static final int COUNTERIDENTITY = 2;
	public static final int ONES = 3;
	public static final int entry_changed = 1;
	public static final int det_changed = 2;
	public static final int eigen_changed = 4;
	public static final int norms_changed = 8;
	public static final int lu_changed = 16;
	public static final int inverse_changed = 32;
	public static final int minmax_changed = 64;
	public static final int all_changed = 127;
	public static final float ZEROF = 0.0f;
	public static final float IDENTITYF = 1.0f;
	public float eij;
	public float entry[][];
	public float trace;
	public float det;
	public float[] minmax;
	public float eigenvector[][];
	public float eigenvalue[];
	public float krylovs[][];
	public float norms[];
	public float lu_entry[][];
	public float inverse_entry[][];
	public int updates;
	public int eigens_calculated;
	public int norms_calculated;

	public Matrix(int TYPE, int rows_num, int columns_num) {
		rows = rows_num;
		columns = columns_num;
		symmetric = rows == columns ? 1 : 0;
		shortest = rows < columns ? rows : columns;
		updates = all_changed;
		switch (TYPE) {
		case (NONE):
		case (ZERO):
			make_zero();
		break;
		case (IDENTITY):
			make_identity();
		break;
		case (COUNTERIDENTITY):
			make_counterIdentity();
		break;
		case (ONES):
			make_ones();
		break;
		}
		eigenvector = new float[rows][columns];
		eigenvalue = new float[rows];
	}

	public void show() {
		int i = rows;
		while(--i >= 0)
			show_v(entry[i]);
		System.out.println(rows + "," + columns);
	}

	public static void show_v (float[] v) {
		int i = v.length;
		while(--i >= 0)
			System.out.print(v[i] + " ");
		System.out.println();
	}

	public float[] minmax() {
		if (minmax!=null && (updates & minmax_changed) == 0)
			return minmax;
		minmax = new float[2];
		if (rows == 0 || columns == 0)
			return minmax;
		minmax[0] = minmax[1] = entry[0][0];
		i = rows;
		while (--i >= 0) {
			j = columns;
			while (--j >= 0) {
				eij = entry[i][j];
				if (eij > minmax[1])
					minmax[1] = eij;
				else if (eij < minmax[0])
					minmax[0] = eij;
			}
		}
		updates &= (all_changed ^ minmax_changed);
		return minmax;
	}

	public float[] minmaxRow(int row) {
		float[] minmaxRow = new float[2];
		if (rows < row || columns == 0)
			return minmaxRow;
		minmaxRow[0] = minmaxRow[1] = entry[row][0];
		j = columns;
		while (--j >= 0) {
			eij = entry[row][j];
			if (eij > minmaxRow[1])
				minmaxRow[1] = eij;
			else if (eij < minmaxRow[0])
				minmaxRow[0] = eij;
		}
		return minmaxRow;
	}

	public float trace() {
		i = shortest;
		trace = 0.0f;
		while (--i >= 0)
			trace += entry[i][i];
		return trace;
	}

	public void make_zero() {
		entry = new float[rows][columns];
		updates &= (all_changed ^ entry_changed);
	}

	public void make_identity() {
		entry = new float[rows][columns];
		i = shortest;
		while (--i >= 0)
			entry[i][i] = IDENTITYF;
		updates &= (all_changed ^ entry_changed);
	}

	public void make_counterIdentity() {
		entry = new float[rows][columns];
		i = shortest;
		while (--i >= 0)
			entry[i][shortest - i - 1] = IDENTITYF;
		updates &= (all_changed ^ entry_changed);
	}

	public void make_ones() {
		entry = new float[rows][columns];
		i = rows;
		while (--i >= 0) {
			Arrays.fill(entry[i], IDENTITYF);
		}
	}

	public static void copy(float[][] source, float[][] dest) {
		int si = dest.length;
		int sj = dest[0].length;
		while (--si >= 0) {
			System.arraycopy(source[si], 0, dest[si], 0, sj);
		}
	}

	public Matrix(float[][] entries) {
		rows = entries.length;
		if (rows != 0)
			columns = entries[0].length;
		else
			columns = 0;
		symmetric = rows == columns ? 1 : 0;
		shortest = rows < columns ? rows : columns;
		entry = new float[rows][columns];
		Matrix.copy(entries,entry);
		eigenvector = new float[rows][columns];
		eigenvalue = new float[rows];
		updates = all_changed;
	}

	public Matrix(Matrix expired) {
		rows = expired.rows;
		columns = expired.columns;
		entry = new float[rows][columns];
		eigenvector = new float[rows][columns];
		eigenvalue = new float[rows];
		Matrix.copy(expired.entry, entry);
		Matrix.copy(expired.eigenvector, eigenvector);
		System.arraycopy(expired.eigenvalue, 0, eigenvalue, 0,
				eigenvalue.length);
		shortest = expired.shortest;
	}

	public Matrix transpose() {
		int j = columns;
		int i = rows;
		Matrix transpose = new Matrix(Matrix.NONE, j, i);
		while (--j >= 0) {
			while (--i >= 0) {
				transpose.entry[j][i] = entry[i][j];
			}
			i = rows;
		}
		return transpose;
	}

}
import java.util.ArrayList;
import processing.core.PVector;

public class Visionary {
	public final static int VERTICESONSCREEN = 60;
	public final static int NOMOD = 0;
	public final static int MOD = 1;
	public int[] gridsz;
	public int mc = 0;
	public Matrix lrf;
	public Matrix gridv;
	public int etype = 0;
	public int offset = 0;

	public final AspektWorld aspektWorld;
	public Graph g1;
	public Adventurer bg;

	public Visionary(AspektWorld aw) {
		aspektWorld = aw;
		gridsz = new int[] { 70,140,210,280,350,420,490,560,630,700 };
		lrf = new Matrix(Matrix.IDENTITY,VERTICESONSCREEN, VERTICESONSCREEN);
	}

	public Graph randomGraph() {
		Graph g = new Graph(aspektWorld);
		ArrayList<Graph.Vertex> knownVertices = new ArrayList<Graph.Vertex>();
		int i = VERTICESONSCREEN;
		while(--i >= 0) {
			Graph.Vertex v = g.new Vertex();
			g.addVertex(v);
			knownVertices.add(v);
			int k,j;
			k = VERTICESONSCREEN-i;
			j = (int) aspektWorld.random(k);
			while(--j >= 0) {
				Graph.Vertex randomNeighbour = knownVertices.get((int)(aspektWorld.random(k)));
				Graph.Edge newEdge = g.new Edge(v, randomNeighbour);
				g.addEdge(newEdge);
			}
		}
		g.crystallizeGraph();
		g1 = g;
		bg = new Adventurer(g1);
		g1 = bg.g;
		return g;
	}

	public int vertexToGrid(Graph.Vertex v) {
		int vn = g1.numbering.get(v);
		if(gridv == null) {
			float maxGrid = gridsz.length;
			Matrix nd = new Matrix(g1.matrix[Graph.MatrixTypeIndex.DEGREE]);
			Arithmetician.scaleToRange(nd, 0.0f, maxGrid-1 ); 
			gridv = nd;
		}
		return (int)gridv.entry[vn][vn];
	}

	public void boxit() {
		int i = g1.v.length;
		while(--i >= 0) {
			Graph.Vertex v = g1.v[i];
			int grid = vertexToGrid(v);
			v.pos.set(boxCoords(nextBox(v.pos, grid), grid)); 
		}
	}

	public int[] nextBox( PVector start, int grid ) {
		int sz = gridsz[grid];
		int x = (int)(start.x/sz); int y = (int)(start.y/sz); int z = (int)(start.z/sz);
		return new int[] { x,y,z };
	}

	public PVector boxCoords( int[] boxp, int grid ) {
		int sz = gridsz[grid];
		PVector j = new PVector(boxp[0]*sz,boxp[1]*sz,boxp[2]*sz);
		return j;
	}

	public Graph randomizeLayoutColorSize( ) {
		int i = g1.v.length;
		while(--i >= 0) {
			Graph.Vertex v = g1.v[i];
			v.color = aspektWorld.color(aspektWorld.random(360),200,200,200);
			v.radius = 3 + aspektWorld.random(40);
			v.pos.set(aspektWorld.random(aspektWorld.width)-aspektWorld.width/2, -aspektWorld.random(aspektWorld.height), aspektWorld.random(aspektWorld.height)-aspektWorld.height/2);
		}
		i = g1.e.length;
		while(--i >= 0) {
			Graph.Edge e = g1.e[i];
		}
		return g1;
	}

	public int[] randomi( int sz, int max ) {
		int[] ri = new int[sz];
		while(--sz >= 0)
			ri[sz] = (int)aspektWorld.random(max);
		return ri;
	}

	public int min(int a, int b) {
		return a < b ? a : b;
	}

	public void embedPlain(Matrix f) {
		Matrix ef = new Matrix(f);
		Arithmetician.scaleToRangeRow(ef, offset+0, -aspektWorld.height*2, aspektWorld.height*2);
		Arithmetician.scaleToRangeRow(ef, offset+1, -aspektWorld.height*2, aspektWorld.height*2);
		Arithmetician.vectorScalarAdd(ef.entry[offset+1], -aspektWorld.height*2);
		Arithmetician.scaleToRangeRow(ef, offset+2, -aspektWorld.height, aspektWorld.height*2);
		Arithmetician.scaleToRangeRow(ef, offset+3, 0, 360);
		Arithmetician.scaleToRangeRow(ef, offset+4, 10, 70);
		int count = min(ef.shortest,g1.v.length);
		while(--count >= 0 ) {
			Graph.Vertex v = g1.v[count];
			v.pos.x = ef.entry[offset+0][count];
			v.pos.y = ef.entry[offset+1][count];
			v.pos.z = ef.entry[offset+2][count];
			v.color = aspektWorld.color(ef.entry[offset+3][count],200,200,200);
			v.radius = ef.entry[offset+4][count];
		}
	}

	public void embedProjected(Matrix f) {
		Matrix ef;
		ef = Algebraicist.projectedMatrix(f, 0.5f);
		//ef = Algebraicist.eigenvectors(ef);
		Arithmetician.scaleToRangeRow(ef, offset+0, -aspektWorld.height*2, aspektWorld.height*2);
		Arithmetician.scaleToRangeRow(ef, offset+1, -aspektWorld.height*2, aspektWorld.height*2);
		Arithmetician.vectorScalarAdd(ef.entry[offset+1], -aspektWorld.height*2);
		Arithmetician.scaleToRangeRow(ef, offset+2, -aspektWorld.height, aspektWorld.height*2);
		Arithmetician.scaleToRangeRow(ef, offset+3, 0, 360);
		Arithmetician.scaleToRangeRow(ef, offset+4, 10, 70);
		int count = min(ef.shortest,g1.v.length);
		while(--count >= 0 ) {
			Graph.Vertex v = g1.v[count];
			v.pos.x = ef.entry[offset+0][count];
			v.pos.y = ef.entry[offset+1][count];
			v.pos.z = ef.entry[offset+2][count];
			v.color = aspektWorld.color(ef.entry[offset+3][count],200,200,200);
			v.radius = ef.entry[offset+4][count];
		}
	}

	public void embedProjectedEigeneered(Matrix f) {
		Matrix ef;
		ef = Algebraicist.projectedMatrix(f, 0.5f);
		Arithmetician.PowerIterationDeflate(ef, 5+offset, 0.00000001, 0);
		ef = new Matrix(ef.eigenvector);
		Arithmetician.scaleToRangeRow(ef, offset+0, -aspektWorld.height*2, aspektWorld.height*2);
		Arithmetician.scaleToRangeRow(ef, offset+1, -aspektWorld.height*2, aspektWorld.height*2);
		Arithmetician.vectorScalarAdd(ef.entry[offset+1], -aspektWorld.height*2);
		Arithmetician.scaleToRangeRow(ef, offset+2, -aspektWorld.height, aspektWorld.height*2);
		Arithmetician.scaleToRangeRow(ef, offset+3, 0, 360);
		Arithmetician.scaleToRangeRow(ef, offset+4, 10, 70);
		int count = min(ef.shortest,g1.v.length);
		while(--count >= 0 ) {
			Graph.Vertex v = g1.v[count];
			v.pos.x = ef.entry[offset+0][count];
			v.pos.y = ef.entry[offset+1][count];
			v.pos.z = ef.entry[offset+2][count];
			v.color = aspektWorld.color(ef.entry[offset+3][count],200,200,200);
			v.radius = ef.entry[offset+4][count];
		}
	}

	public static void show_v (float[] v) {
		int i = v.length;
		while(--i >= 0)
			System.out.print(v[i] + " ");
		System.out.println();
	}

	public Graph randomFilterEmbed() {
		offset = 0;
		if(bg == null)
			bg = new Adventurer(g1);
		if(g1.matrix[0] == null)
			bg.make_all_basic_matrices();
		// Better random filter code
		Matrix rf = bg.make_random_filter((int)aspektWorld.random(Graph.MatrixTypeIndex.MATRIX_FIELD_COUNT), 
				(int)aspektWorld.random(Graph.MatrixTypeIndex.MATRIX_FIELD_COUNT), 
				(int)aspektWorld.random(13), 
				(int)aspektWorld.random(Arithmetician.MatrixOperation.MatrixBinary.BINARY_FIELD_COUNT));
		// Older rf code
		/**
		 * Matrix rf = bg.make_random_Matrix_filter(NOMOD, 

				randomi(5, Graph.MatrixTypeIndex.MATRIX_FIELD_COUNT), 
				randomi(5, Arithmetician.MatrixOperation.MatrixUnary.UNARY_FIELD_COUNT), 
				randomi(4, Arithmetician.MatrixOperation.MatrixBinary.BINARY_FIELD_COUNT), 
				new int[] { 0, 0 }, new float[] { aspektWorld.random(10), 10+aspektWorld.random(20), 20+aspektWorld.random(30) });
		 **/
		etype = 4;
		lrf = rf;
		embedProjectedEigeneered(rf);
		return g1;
	}

	public Graph basicMatrixEmbed() {
		offset = 0;
		if(bg == null)
			bg = new Adventurer(g1);
		if(g1.matrix[0] == null)
			bg.make_all_basic_matrices();
		//System.out.println(mc);
		Matrix rf = g1.matrix[mc];
		mc ++;
		if(mc >= Graph.MatrixTypeIndex.MATRIX_FIELD_COUNT)
			mc = 0;
		lrf = rf;
		etype = 1;
		embedPlain(rf);
		return g1;
	}

	public Graph basicMatrixProjectedEmbed() {
		offset = 0;
		if(bg == null)
			bg = new Adventurer(g1);
		if(g1.matrix[0] == null)
			bg.make_all_basic_matrices();
		//System.out.println(mc);
		Matrix rf = g1.matrix[mc];
		mc ++;
		if(mc >= Graph.MatrixTypeIndex.MATRIX_FIELD_COUNT)
			mc = 0;
		lrf = rf;
		etype = 2;
		embedProjected(rf);
		return g1;
	}

	public Graph eigeneeringEmbed() {
		offset = 0;
		if(bg == null)
			bg = new Adventurer(g1);
		if(g1.matrix[0] == null)
			bg.make_all_basic_matrices();
		//System.out.println(mc);
		Matrix rf = g1.matrix[mc];
		mc ++;
		if(mc >= Graph.MatrixTypeIndex.MATRIX_FIELD_COUNT)
			mc = 0;
		lrf = rf;
		etype = 3;
		embedProjectedEigeneered(rf);
		return g1;
	}

	public void reEmbedWithSameFilter() {
		// TODO Auto-generated method stub
		switch(etype) {
		case(1):
			embedPlain(lrf);
		break;
		case(2):
			embedProjected(lrf);
		break;
		case(3):case(4):default:
			embedProjectedEigeneered(lrf);
		}
	}
}