List of usage examples for org.apache.commons.math3.linear EigenDecomposition getDeterminant
public double getDeterminant()
From source file:lanchester.TestEigen.java
public static void main(String[] args) { mat = getMat(4);/*w w w . ja v a 2 s .c om*/ // fillMat(mat); eigen = new EigenDecomposition(mat); double det = eigen.getDeterminant(); double[] eVals = eigen.getRealEigenvalues(); if (eigen.hasComplexEigenvalues()) { System.out.println("Complex eigenvalues"); } Array2DRowRealMatrix eVectors = (Array2DRowRealMatrix) eigen.getV(); EigenDecomposition eigen2 = new EigenDecomposition(eVectors); double det2 = eigen2.getDeterminant(); System.out.println("Det 1 = " + det + " and Det 2 = " + det2); }
From source file:lanchester.MultiArena.java
public void step() { int numFoes = forces.size(); if (isMyTurn == null) { isMyTurn = new boolean[numFoes][numFoes]; stanceArray = new int[numFoes][numFoes]; currentFloor = new double[numFoes][numFoes]; for (int i1 = 0; i1 < numFoes; i1++) { int ind1 = forceMap.get(forces.get(i1)); for (int i2 = 0; i2 < numFoes; i2++) { int ind2 = forceMap.get(forces.get(i2)); isMyTurn[i1][i2] = true; if (i1 == i2) { stanceArray[i1][i2] = AthenaConstants.ALLIED_POSTURE; currentFloor[i1][i2] = 0.; } else { stanceArray[i1][i2] = initializeStance(forces.get(i1), forces.get(i2)); setFloor(i1, i2);/* ww w .ja v a 2 s. c om*/ } } } } Array2DRowRealMatrix mat = getMat(); EigenDecomposition eigen = new EigenDecomposition(mat); double det = eigen.getDeterminant(); double[] eVals = eigen.getRealEigenvalues(); // for(int i1=0;i1<eVals.length;i1++){ // System.out.println("eVals["+i1+"] = "+eVals[i1]); // } if (eigen.hasComplexEigenvalues()) { System.out.println("Complex eigenvalues"); for (int i1 = 0; i1 < forces.size(); i1++) { AthenaForce f = forces.get(i1); System.out.println(f.getName() + " has " + f.getForceSize() + " forces remaining"); } } double[] initialNums = getInitialNumbers(forces); Array2DRowRealMatrix eVectors = (Array2DRowRealMatrix) eigen.getV(); LUDecomposition lu = new LUDecomposition(eVectors); double det2 = lu.getDeterminant(); double[] coeffs = new double[numFoes]; for (int i1 = 0; i1 < numFoes; i1++) { Array2DRowRealMatrix tmpMat = (Array2DRowRealMatrix) eVectors.copy(); tmpMat.setColumn(i1, initialNums); LUDecomposition tmpLU = new LUDecomposition(tmpMat); double tmpDet = tmpLU.getDeterminant(); coeffs[i1] = tmpDet / det2; } MultiTimeStep currentStep = new MultiTimeStep(numFoes); currentTime += timeStep; currentStep.setTime(currentTime); for (int i1 = 0; i1 < numFoes; i1++) { double updatedForce = 0.; for (int i2 = 0; i2 < numFoes; i2++) { updatedForce += coeffs[i2] * eVectors.getEntry(i1, i2) * Math.exp(eVals[i2] * timeStep); // updatedForce+=coeffs[i2]*eVectors.getEntry(i2, i1)*Math.exp(eVals[i2]*timeStep); // updatedForce+=coeffs[i1]*eVectors.getEntry(i2, i1)*Math.exp(eVals[i1]*timeStep); } forces.get(i1).updateForce(updatedForce); currentStep.setForceNumber(updatedForce, i1); } history.add(currentStep); // eVectors. // this.currentTime++; // Truncator truncator = new Truncator(); if (true) { // System.out.println("time = " + time); } }
From source file:lanchester.MultiArena2.java
public void step() { boolean aboveFloor = true; double currentCycle = 0.; int numFoes = forces.size(); if (isMyTurn == null) { isMyTurn = new boolean[numFoes][numFoes]; stanceArray = new int[numFoes][numFoes]; currentFloor = new double[numFoes][numFoes]; for (int i1 = 0; i1 < numFoes; i1++) { int ind1 = forceMap.get(forces.get(i1)); for (int i2 = 0; i2 < numFoes; i2++) { int ind2 = forceMap.get(forces.get(i2)); isMyTurn[i1][i2] = true; if (i1 == i2) { stanceArray[i1][i2] = AthenaConstants.ALLIED_POSTURE; currentFloor[i1][i2] = 0.; } else { stanceArray[i1][i2] = initializeStance(forces.get(i1), forces.get(i2)); setFloor(i1, i2);// www . j a v a 2 s.c o m } } } } Array2DRowRealMatrix mat = getMat(); EigenDecomposition eigen = new EigenDecomposition(mat); double det = eigen.getDeterminant(); double[] eVals = eigen.getRealEigenvalues(); // for(int i1=0;i1<eVals.length;i1++){ // System.out.println("eVals["+i1+"] = "+eVals[i1]); // } if (eigen.hasComplexEigenvalues()) { System.out.println("Complex eigenvalues"); for (int i1 = 0; i1 < forces.size(); i1++) { AthenaForce f = forces.get(i1); System.out.println(f.getName() + " has " + f.getForceSize() + " forces remaining"); } } double[] initialNums = getInitialNumbers(forces); Array2DRowRealMatrix eVectors = (Array2DRowRealMatrix) eigen.getV(); LUDecomposition lu = new LUDecomposition(eVectors); double det2 = lu.getDeterminant(); double[] coeffs = new double[numFoes]; for (int i1 = 0; i1 < numFoes; i1++) { Array2DRowRealMatrix tmpMat = (Array2DRowRealMatrix) eVectors.copy(); tmpMat.setColumn(i1, initialNums); LUDecomposition tmpLU = new LUDecomposition(tmpMat); double tmpDet = tmpLU.getDeterminant(); coeffs[i1] = tmpDet / det2; } aboveFloor = true; int cntr = 0; int numGone; do { MultiTimeStep currentStep = new MultiTimeStep(numFoes); currentTime += timeStep; currentCycle += timeStep; currentStep.setTime(currentTime); numGone = 0; for (int i1 = 0; i1 < numFoes; i1++) { double updatedForce = 0.; if (forces.get(i1).getForceSize() > lb) { for (int i2 = 0; i2 < numFoes; i2++) { // updatedForce += coeffs[i2] * eVectors.getEntry(i1, i2) * Math.exp(eVals[i2] * timeStep); updatedForce += coeffs[i2] * eVectors.getEntry(i1, i2) * Math.exp(eVals[i2] * currentCycle); if (updatedForce < 1.) { updatedForce = 0.; numGone++; } // updatedForce+=coeffs[i2]*eVectors.getEntry(i2, i1)*Math.exp(eVals[i2]*timeStep); // updatedForce+=coeffs[i1]*eVectors.getEntry(i2, i1)*Math.exp(eVals[i1]*timeStep); } } else { updatedForce = lb / 2.; numGone++; } forces.get(i1).updateForce(updatedForce); currentStep.setForceNumber(updatedForce, i1); } history.add(currentStep); aboveFloor = checkAboveFloors(); cntr++; } while (aboveFloor && cntr < 2000 && (numFoes - numGone) > 1); for (int i1 = 0; i1 < numFoes; i1++) { for (int i2 = 0; i2 < numFoes; i2++) { if (i1 == i2) { stanceArray[i1][i2] = AthenaConstants.ALLIED_POSTURE; currentFloor[i1][i2] = 0.; } else { stanceArray[i1][i2] = initializeStance(forces.get(i1), forces.get(i2)); setFloor(i1, i2); } } } // eVectors. // this.currentTime++; // Truncator truncator = new Truncator(); if (numFoes - numGone == 1) { loneSurvivor = true; // System.out.println("time = " + time); } }
From source file:MultivariateNormalDistribution.java
/** * Creates a multivariate normal distribution with the given mean vector and * covariance matrix.//from ww w .j ava2 s . com * <br/> * The number of dimensions is equal to the length of the mean vector * and to the number of rows and columns of the covariance matrix. * It is frequently written as "p" in formulae. * * @param rng Random Number Generator. * @param means Vector of means. * @param covariances Covariance matrix. * @throws DimensionMismatchException if the arrays length are * inconsistent. * @throws SingularMatrixException if the eigenvalue decomposition cannot * be performed on the provided covariance matrix. * @throws NonPositiveDefiniteMatrixException if any of the eigenvalues is * negative. */ public MultivariateNormalDistribution(RandomGenerator rng, final double[] means, final double[][] covariances) throws SingularMatrixException, DimensionMismatchException, NonPositiveDefiniteMatrixException { super(rng, means.length); final int dim = means.length; if (covariances.length != dim) { throw new DimensionMismatchException(covariances.length, dim); } for (int i = 0; i < dim; i++) { if (dim != covariances[i].length) { throw new DimensionMismatchException(covariances[i].length, dim); } } this.means = MathArrays.copyOf(means); covarianceMatrix = new Array2DRowRealMatrix(covariances); // Covariance matrix eigen decomposition. final EigenDecomposition covMatDec = new EigenDecomposition(covarianceMatrix); // Compute and store the inverse. covarianceMatrixInverse = covMatDec.getSolver().getInverse(); // Compute and store the determinant. covarianceMatrixDeterminant = covMatDec.getDeterminant(); // Eigenvalues of the covariance matrix. final double[] covMatEigenvalues = covMatDec.getRealEigenvalues(); for (int i = 0; i < covMatEigenvalues.length; i++) { if (covMatEigenvalues[i] < 0) { throw new NonPositiveDefiniteMatrixException(covMatEigenvalues[i], i, 0); } } // Matrix where each column is an eigenvector of the covariance matrix. final Array2DRowRealMatrix covMatEigenvectors = new Array2DRowRealMatrix(dim, dim); for (int v = 0; v < dim; v++) { final double[] evec = covMatDec.getEigenvector(v).toArray(); covMatEigenvectors.setColumn(v, evec); } final RealMatrix tmpMatrix = covMatEigenvectors.transpose(); // Scale each eigenvector by the square root of its eigenvalue. for (int row = 0; row < dim; row++) { final double factor = FastMath.sqrt(covMatEigenvalues[row]); for (int col = 0; col < dim; col++) { tmpMatrix.multiplyEntry(row, col, factor); } } samplingMatrix = covMatEigenvectors.multiply(tmpMatrix); }
From source file:lanchester.MultiArena3.java
public void step() { boolean aboveFloor = true; double currentCycle = 0.; int numFoes = forces.size(); if (isMyTurn == null) { isMyTurn = new boolean[numFoes][numFoes]; stanceArray = new int[numFoes][numFoes]; currentFloor = new double[numFoes][numFoes]; currentCeiling = new double[numFoes][numFoes]; for (int i1 = 0; i1 < numFoes; i1++) { int ind1 = forceMap.get(forces.get(i1)); for (int i2 = 0; i2 < numFoes; i2++) { int ind2 = forceMap.get(forces.get(i2)); isMyTurn[i1][i2] = true; if (i1 == i2) { stanceArray[i1][i2] = AthenaConstants.ALLIED_POSTURE; currentFloor[i1][i2] = 0.; currentCeiling[i1][i2] = 100.; } else { stanceArray[i1][i2] = initializeStance(forces.get(i1), forces.get(i2)); setFloor(i1, i2);/* w w w. j a v a 2s . co m*/ setCeiling(i1, i2); } } } } Array2DRowRealMatrix mat = getMat(); EigenDecomposition eigen = new EigenDecomposition(mat); double det = eigen.getDeterminant(); double[] eVals = eigen.getRealEigenvalues(); if (eigen.hasComplexEigenvalues()) { System.out.println("Complex eigenvalues"); for (int i1 = 0; i1 < forces.size(); i1++) { MultiForce f = forces.get(i1); System.out.println(f.getName() + " has " + f.getNumber() + " forces remaining"); } } double[] initialNums = getInitialNumbers(forces); Array2DRowRealMatrix eVectors = (Array2DRowRealMatrix) eigen.getV(); LUDecomposition lu = new LUDecomposition(eVectors); double det2 = lu.getDeterminant(); double[] coeffs = new double[numFoes]; for (int i1 = 0; i1 < numFoes; i1++) { Array2DRowRealMatrix tmpMat = (Array2DRowRealMatrix) eVectors.copy(); tmpMat.setColumn(i1, initialNums); LUDecomposition tmpLU = new LUDecomposition(tmpMat); double tmpDet = tmpLU.getDeterminant(); coeffs[i1] = tmpDet / det2; } aboveFloor = true; boolean belowCeiling = true; int cntr = 0; int numGone; do { timeStep = determineTimeStep(); MultiTimeStep currentStep = new MultiTimeStep(numFoes); currentTime += timeStep; currentCycle += timeStep; currentStep.setTime(currentTime); numGone = 0; for (int i1 = 0; i1 < numFoes; i1++) { double updatedForce = 0.; if (stillAlive[i1]) { for (int i2 = 0; i2 < numFoes; i2++) { updatedForce += coeffs[i2] * eVectors.getEntry(i1, i2) * Math.exp(eVals[i2] * currentCycle); } if (updatedForce < 1.) { updatedForce = lb; stillAlive[i1] = false; numGone++; } } else { numGone++; updatedForce = lb; } forces.get(i1).updateForce(updatedForce); currentStep.setForceNumber(updatedForce, i1); } history.add(currentStep); aboveFloor = checkAboveFloors(); belowCeiling = checkBelowCeilings(); cntr++; } while (aboveFloor && belowCeiling && cntr < 2000 && (numFoes - numGone) > 1); for (int i1 = 0; i1 < numFoes; i1++) { for (int i2 = 0; i2 < numFoes; i2++) { if (i1 == i2) { stanceArray[i1][i2] = AthenaConstants.ALLIED_POSTURE; currentFloor[i1][i2] = 0.; } else { stanceArray[i1][i2] = initializeStance(forces.get(i1), forces.get(i2)); setFloor(i1, i2); } } } // eVectors. // this.currentTime++; // Truncator truncator = new Truncator(); if (numFoes - numGone == 1) { loneSurvivor = true; // System.out.println("time = " + time); } }
From source file:lanchester.MultiArena4.java
public void step() { boolean aboveFloor = true; double currentCycle = 0.; int numFoes = forces.size(); if (isMyTurn == null) { isMyTurn = new boolean[numFoes][numFoes]; stanceArray = new int[numFoes][numFoes]; currentFloor = new double[numFoes][numFoes]; currentCeiling = new double[numFoes][numFoes]; for (int i1 = 0; i1 < numFoes; i1++) { // int ind1 = forceMap.get(forces.get(i1)); for (int i2 = 0; i2 < numFoes; i2++) { // int ind2 = forceMap.get(forces.get(i2)); isMyTurn[i1][i2] = true; if (i1 == i2) { stanceArray[i1][i2] = AthenaConstants.ALLIED_POSTURE; currentFloor[i1][i2] = 0.; currentCeiling[i1][i2] = 1000000.; } else { stanceArray[i1][i2] = initializeStance(forces.get(i1), forces.get(i2)); // setStances(i1,i2); setFloor(i1, i2);//from www . j ava 2 s. c om setCeiling(i1, i2); } } } } Array2DRowRealMatrix mat = getMat(); EigenDecomposition eigen = new EigenDecomposition(mat); double det = eigen.getDeterminant(); double[] eVals = eigen.getRealEigenvalues(); if (eigen.hasComplexEigenvalues()) { System.out.println("Complex eigenvalues"); for (int i1 = 0; i1 < forces.size(); i1++) { MultiForce f = forces.get(i1); System.out.println(f.getName() + " has " + f.getNumber() + " forces remaining"); } } double[] initialNums = getInitialNumbers(forces); Array2DRowRealMatrix eVectors = (Array2DRowRealMatrix) eigen.getV(); LUDecomposition lu = new LUDecomposition(eVectors); double det2 = lu.getDeterminant(); double[] coeffs = new double[numFoes]; for (int i1 = 0; i1 < numFoes; i1++) { Array2DRowRealMatrix tmpMat = (Array2DRowRealMatrix) eVectors.copy(); tmpMat.setColumn(i1, initialNums); LUDecomposition tmpLU = new LUDecomposition(tmpMat); double tmpDet = tmpLU.getDeterminant(); coeffs[i1] = tmpDet / det2; } aboveFloor = true; boolean belowCeiling = true; int cntr = 0; int numGone; do { timeStep = determineTimeStep(); MultiTimeStep currentStep = new MultiTimeStep(numFoes); currentTime += timeStep; currentCycle += timeStep; currentStep.setTime(currentTime); numGone = 0; for (int i1 = 0; i1 < numFoes; i1++) { double updatedForce = 0.; // if (forces.get(i1).getForceSize() > lb) { if (stillAlive[i1]) { for (int i2 = 0; i2 < numFoes; i2++) { // updatedForce += coeffs[i2] * eVectors.getEntry(i1, i2) * Math.exp(eVals[i2] * timeStep); updatedForce += coeffs[i2] * eVectors.getEntry(i1, i2) * Math.exp(eVals[i2] * currentCycle); // updatedForce+=coeffs[i2]*eVectors.getEntry(i2, i1)*Math.exp(eVals[i2]*timeStep); // updatedForce+=coeffs[i1]*eVectors.getEntry(i2, i1)*Math.exp(eVals[i1]*timeStep); } if (updatedForce < 1.) { updatedForce = lb; stillAlive[i1] = false; numGone++; } } else { updatedForce = lb; numGone++; } forces.get(i1).updateForce(updatedForce); currentStep.setForceNumber(updatedForce, i1); // for (int i2 = 0; i2 < numFoes; i1++) { // if (i1 != i2) { // // } // } } history.add(currentStep); aboveFloor = checkAboveFloors(); belowCeiling = checkBelowCeilings(); cntr++; } while (aboveFloor && belowCeiling && cntr < 2000 && (numFoes - numGone) > 1); for (int i1 = 0; i1 < numFoes; i1++) { for (int i2 = 0; i2 < numFoes; i2++) { if (i1 == i2) { stanceArray[i1][i2] = AthenaConstants.ALLIED_POSTURE; currentFloor[i1][i2] = 0.; } else { // stanceArray[i1][i2] = initializeStance(forces.get(i1), forces.get(i2)); setStances(i1, i2); setFloor(i1, i2); setCeiling(i1, i2); } } } // eVectors. // this.currentTime++; // Truncator truncator = new Truncator(); if (numFoes - numGone == 1) { loneSurvivor = true; // System.out.println("time = " + time); } }
From source file:net2.N2MultiArena.java
public void step() { boolean aboveFloor = true; double currentCycle = 0.; int numFoes = forces.size(); System.out.println("Num foes = " + numFoes); if (isMyTurn == null) { isMyTurn = new boolean[numFoes][numFoes]; stanceArray = new int[numFoes][numFoes]; currentFloor = new double[numFoes][numFoes]; currentCeiling = new double[numFoes][numFoes]; for (int i1 = 0; i1 < numFoes; i1++) { // int ind1 = forceMap.get(forces.get(i1)); for (int i2 = 0; i2 < numFoes; i2++) { // int ind2 = forceMap.get(forces.get(i2)); isMyTurn[i1][i2] = true; if (i1 == i2) { stanceArray[i1][i2] = AthenaConstants.ALLIED_POSTURE; currentFloor[i1][i2] = 0.; currentCeiling[i1][i2] = 1000000.; } else { stanceArray[i1][i2] = initializeStance(forces.get(i1), forces.get(i2)); setFloor(i1, i2);/*from w w w. j a v a 2 s. c o m*/ setCeiling(i1, i2); } } } } Array2DRowRealMatrix mat = getMat(); EigenDecomposition eigen = new EigenDecomposition(mat); double det = eigen.getDeterminant(); double[] eVals = eigen.getRealEigenvalues(); if (eigen.hasComplexEigenvalues()) { System.out.println("Complex eigenvalues"); for (int i1 = 0; i1 < forces.size(); i1++) { N2ForceUnit f = forces.get(i1); System.out.println(f.getName() + " has " + f.getNumber() + " forces remaining"); } } double[] initialNums = getInitialNumbers(forces); Array2DRowRealMatrix eVectors = (Array2DRowRealMatrix) eigen.getV(); LUDecomposition lu = new LUDecomposition(eVectors); double det2 = lu.getDeterminant(); double[] coeffs = new double[numFoes]; for (int i1 = 0; i1 < numFoes; i1++) { Array2DRowRealMatrix tmpMat = (Array2DRowRealMatrix) eVectors.copy(); tmpMat.setColumn(i1, initialNums); LUDecomposition tmpLU = new LUDecomposition(tmpMat); double tmpDet = tmpLU.getDeterminant(); coeffs[i1] = tmpDet / det2; } aboveFloor = true; boolean belowCeiling = true; // int cntr = 0; int numGone; // do { timeStep = determineTimeStep(); MultiTimeStep currentStep = new MultiTimeStep(numFoes); currentTime += timeStep; currentCycle += timeStep; currentStep.setTime(currentTime); numGone = 0; for (int i1 = 0; i1 < numFoes; i1++) { double updatedForce = 0.; // if (forces.get(i1).getForceSize() > lb) { if (stillAlive[i1]) { for (int i2 = 0; i2 < numFoes; i2++) { // updatedForce += coeffs[i2] * eVectors.getEntry(i1, i2) * Math.exp(eVals[i2] * timeStep); updatedForce += coeffs[i2] * eVectors.getEntry(i1, i2) * Math.exp(eVals[i2] * currentCycle); } if (updatedForce < 1.) { updatedForce = lb; stillAlive[i1] = false; numGone++; } } else { updatedForce = lb; numGone++; } forces.get(i1).updateForce(updatedForce); currentStep.setForceNumber(updatedForce, i1); } history.add(currentStep); aboveFloor = checkAboveFloors(); belowCeiling = checkBelowCeilings(); // cntr++; // } while (aboveFloor && belowCeiling && cntr < 2000 && (numFoes - numGone) > 1); for (int i1 = 0; i1 < numFoes; i1++) { for (int i2 = 0; i2 < numFoes; i2++) { if (i1 == i2) { stanceArray[i1][i2] = AthenaConstants.ALLIED_POSTURE; currentFloor[i1][i2] = 0.; } else { // stanceArray[i1][i2] = initializeStance(forces.get(i1), forces.get(i2)); setStances(i1, i2); setFloor(i1, i2); setCeiling(i1, i2); } } } if (numFoes - numGone == 1) { loneSurvivor = true; // System.out.println("time = " + time); } }
From source file:org.eclipse.dataset.LinearAlgebra.java
/** * @param a/*w w w . j ava 2 s . com*/ * @return determinant of dataset */ public static double calcDeterminant(Dataset a) { EigenDecomposition evd = new EigenDecomposition(createRealMatrix(a)); return evd.getDeterminant(); }
From source file:org.rhwlab.dispim.nucleus.NucleusData.java
public Ellipse2d ellipse(int xi, int yi, int zi, Coeff coef) { Array2DRowRealMatrix Q = new Array2DRowRealMatrix(3, 3); Q.setEntry(0, 0, coef.A);//from w ww . ja v a 2s . c o m Q.setEntry(1, 0, coef.B / 2); Q.setEntry(0, 1, coef.B / 2); Q.setEntry(1, 1, coef.C); Q.setEntry(2, 0, coef.D / 2); Q.setEntry(0, 2, coef.D / 2); Q.setEntry(2, 1, coef.E / 2); Q.setEntry(1, 2, coef.E / 2); Q.setEntry(2, 2, coef.F); EigenDecomposition ed = new EigenDecomposition(Q); double detQ = ed.getDeterminant(); RealMatrix rm = new Array2DRowRealMatrix(2, 2); rm.setEntry(0, 0, coef.A); rm.setEntry(1, 1, coef.C); rm.setEntry(0, 1, coef.B / 2.0); rm.setEntry(1, 0, coef.B / 2.0); EigenDecomposition eigenDecomp = new EigenDecomposition(rm); double detA33 = eigenDecomp.getDeterminant(); double[] eigenValues = eigenDecomp.getRealEigenvalues(); //System.out.printf("Eigenvalues: %f,%f\n",eigenValues[0],eigenValues[1]) ; RealVector eigenvector0 = eigenDecomp.getEigenvector(0); Ellipse2d e = new Ellipse2d(); double cot2theta = (coef.A - coef.C) / coef.B; if (Double.isFinite(cot2theta)) { double d = Math.sqrt(1.0 + cot2theta * cot2theta); double cos2theta = cot2theta / d; e.cosine = Math.sqrt((1.0 + cos2theta) / 2.0); e.sine = Math.sqrt((1.0 - cos2theta) / 2.0); } else { e.cosine = 1.0; e.sine = 0.0; } double dd = coef.B * coef.B - 4.0 * coef.A * coef.C; double xc = (2.0 * coef.C * coef.D - coef.B * coef.E) / dd; double yc = (2.0 * coef.A * coef.E - coef.B * coef.D) / dd; // System.out.printf("dd=%f,xc=%f,xcn=%f,yc=%f,ycn%f\n",dd,xc,xcn,yc,ycn); double[] ce = this.getCenter(); e.x = ce[xi] + xc; e.y = ce[yi] + yc; // double f = -detQ/detA33 * 2; double f = -detQ / detA33; double a = eigenValues[0] / f; double b = eigenValues[1] / f; if (a <= 0.0 || b <= 0.0) { return null; } e.a = 1.0 / Math.sqrt(a); e.b = 1.0 / Math.sqrt(b); //System.out.printf("detQ=%e,detA33=%e,f=%f,a=%e,b=%e\n",detQ,detA33,f,e.a,e.b); //System.out.printf("eigenValues (%f,%f)\n",eigenValues[0],eigenValues[1]); e.cosine = eigenvector0.getEntry(0); e.sine = eigenvector0.getEntry(1); e.low[xi] = (long) (e.x - e.a); e.low[yi] = (long) (e.y - e.b); e.low[zi] = 0; e.high[xi] = (long) (e.x + e.a); e.high[yi] = (long) (e.y + e.b); e.high[zi] = 0; //System.out.printf("Ellipse: a=%f,b=%f,x=%f,y=%f\n", e.a,e.b,e.x,e.y); return e; }
From source file:xyz.lejon.sampling.FastMultivariateNormalDistribution.java
/** * Creates a multivariate normal distribution with the given mean vector and * covariance matrix./*from w w w.j av a2 s . com*/ * <br/> * The number of dimensions is equal to the length of the mean vector * and to the number of rows and columns of the covariance matrix. * It is frequently written as "p" in formulae. * * @param rng Random Number Generator. * @param means Vector of means. * @param covariances Covariance matrix. * @throws DimensionMismatchException if the arrays length are * inconsistent. * @throws SingularMatrixException if the eigenvalue decomposition cannot * be performed on the provided covariance matrix. * @throws NonPositiveDefiniteMatrixException if any of the eigenvalues is * negative. */ public FastMultivariateNormalDistribution(RandomGenerator rng, final double[] means, final double[][] covariances) throws SingularMatrixException, DimensionMismatchException, NonPositiveDefiniteMatrixException { super(rng, means.length); final int dim = means.length; if (covariances.length != dim) { throw new DimensionMismatchException(covariances.length, dim); } for (int i = 0; i < dim; i++) { if (dim != covariances[i].length) { throw new DimensionMismatchException(covariances[i].length, dim); } } this.means = MathArrays.copyOf(means); covarianceMatrix = new Array2DRowRealMatrix(covariances); // Covariance matrix eigen decomposition. final EigenDecomposition covMatDec = new EigenDecomposition(covarianceMatrix); // Compute and store the inverse. covarianceMatrixInverse = covMatDec.getSolver().getInverse(); // Compute and store the determinant. covarianceMatrixDeterminant = covMatDec.getDeterminant(); // Eigenvalues of the covariance matrix. final double[] covMatEigenvalues = covMatDec.getRealEigenvalues(); for (int i = 0; i < covMatEigenvalues.length; i++) { if (covMatEigenvalues[i] < 0) { throw new NonPositiveDefiniteMatrixException(covMatEigenvalues[i], i, 0); } } // Matrix where each column is an eigenvector of the covariance matrix. final Array2DRowRealMatrix covMatEigenvectors = new Array2DRowRealMatrix(dim, dim); final Array2DRowRealMatrix tmpMatrix = new Array2DRowRealMatrix(dim, dim); for (int v = 0; v < dim; v++) { final double factor = FastMath.sqrt(covMatEigenvalues[v]); final double[] evec = covMatDec.getEigenvector(v).toArray(); covMatEigenvectors.setColumn(v, evec); tmpMatrix.setRow(v, evec); for (int col = 0; col < dim; col++) { tmpMatrix.multiplyEntry(v, col, factor); } } samplingMatrix = covMatEigenvectors.multiply(tmpMatrix); }