List of usage examples for org.apache.commons.math3.linear ArrayRealVector setEntry
@Override public void setEntry(int index, double value) throws OutOfRangeException
From source file:gamlss.algorithm.Gamlss.java
/** * Main method.//from w w w . j a va2s . c o m * @param args - command-line arguments */ public static void main(final String[] args) { //String fileName = "Data/dataReduced.csv"; String fileName = "Data/oil.csv"; // String fileName = "Data/sp.csv"; //String fileName = "Data/dataReduced.csv"; CSVFileReader readData = new CSVFileReader(fileName); readData.readFile(); ArrayList<String> data = readData.storeValues; ArrayRealVector y = new ArrayRealVector(data.size()); BlockRealMatrix muX = new BlockRealMatrix(data.size(), 1); BlockRealMatrix sigmaX = new BlockRealMatrix(data.size(), 1); BlockRealMatrix nuX = new BlockRealMatrix(data.size(), 1); BlockRealMatrix tauX = new BlockRealMatrix(data.size(), 1); ArrayRealVector w = new ArrayRealVector(data.size()); BlockRealMatrix muS = new BlockRealMatrix(data.size(), 1); BlockRealMatrix sigmaS = new BlockRealMatrix(data.size(), 1); BlockRealMatrix nuS = new BlockRealMatrix(data.size(), 1); BlockRealMatrix tauS = new BlockRealMatrix(data.size(), 1); for (int i = 0; i < data.size(); i++) { String[] line = data.get(i).split(","); y.setEntry(i, Double.parseDouble(line[0])); muX.setEntry(i, 0, Double.parseDouble(line[1])); muS.setEntry(i, 0, Double.parseDouble(line[1])); sigmaX.setEntry(i, 0, Double.parseDouble(line[1])); sigmaS.setEntry(i, 0, Double.parseDouble(line[1])); nuX.setEntry(i, 0, Double.parseDouble(line[1])); nuS.setEntry(i, 0, Double.parseDouble(line[1])); tauX.setEntry(i, 0, Double.parseDouble(line[1])); tauS.setEntry(i, 0, Double.parseDouble(line[1])); } Hashtable<Integer, BlockRealMatrix> designMatrices = new Hashtable<Integer, BlockRealMatrix>(); designMatrices.put(DistributionSettings.MU, muX); designMatrices.put(DistributionSettings.SIGMA, sigmaX); designMatrices.put(DistributionSettings.NU, nuX); designMatrices.put(DistributionSettings.TAU, tauX); HashMap<Integer, BlockRealMatrix> smoothMatrices = new HashMap<Integer, BlockRealMatrix>(); smoothMatrices.put(DistributionSettings.MU, muS); smoothMatrices.put(DistributionSettings.SIGMA, sigmaS); smoothMatrices.put(DistributionSettings.NU, nuS); smoothMatrices.put(DistributionSettings.TAU, tauS); //smoothMatrices.put(DistributionSettings.MU, null); //smoothMatrices.put(DistributionSettings.SIGMA, null); //smoothMatrices.put(DistributionSettings.NU, null); //smoothMatrices.put(DistributionSettings.TAU, null); DistributionSettings.DISTR = DistributionSettings.SST; Controls.GLOB_DEVIANCE_TOL = 5500; Controls.INTER = 50;//only for the PB smoother Controls.SMOOTHER = Controls.PB; //or PB Controls.IS_SVD = true; Controls.BIG_DATA = true; Controls.JAVA_OPTIMIZATION = false; Controls.GAMLSS_NUM_CYCLES = 50; //Gamlss gamlss = new Gamlss(y, designMatrices, null); Gamlss gamlss = new Gamlss(y, designMatrices, null); //Gamlss gamlss = new Gamlss(y, null, smoothMatrices); gamlss.saveFittedDistributionParameters("Data/oilresults.csv"); }
From source file:ImageProc.java
public static float[][] colormapJet(int nClass) { int n;//from ww w .j a v a 2 s .com n = (int) Math.ceil((float) nClass / 4); ArrayRealVector u, r, g, b; RealVector R, G, B; u = new ArrayRealVector(3 * n - 1); for (int i = 0; i < n; i++) { u.setEntry(i, (i + 1.0) / n); u.setEntry(u.getDimension() - i - 1, (i + 1.0) / n); } u.setSubVector(n, new ArrayRealVector(n, 1)); g = new ArrayRealVector(u.getDimension()); float m; m = (float) Math.ceil((float) n / 2); if (nClass % 4 == 1) m = m - 1; for (int i = 0; i < g.getDimension(); i++) { g.setEntry(i, (i + 1 + m)); } r = g.add(new ArrayRealVector(g.getDimension(), n)); b = g.subtract(new ArrayRealVector(g.getDimension(), n)); R = new ArrayRealVector(); G = new ArrayRealVector(); B = new ArrayRealVector(); for (int i = 0; i < r.getDimension(); i++) { if (r.getEntry(i) <= nClass) R = R.append(r.getEntry(i)); } for (int i = 0; i < g.getDimension(); i++) { if (g.getEntry(i) <= nClass) G = G.append(g.getEntry(i)); } for (int i = 0; i < b.getDimension(); i++) { if (b.getEntry(i) >= 1) B = B.append(b.getEntry(i)); } float[][] J = new float[nClass][3]; int index; for (int i = 0; i < R.getDimension(); i++) { index = (int) R.getEntry(i); J[index - 1][0] = (float) u.getEntry(i); } for (int i = 0; i < G.getDimension(); i++) { index = (int) G.getEntry(i); J[index - 1][1] = (float) u.getEntry(i); } for (int i = u.getDimension() - B.getDimension(), j = 0; i < u.getDimension(); i++, j++) { index = (int) B.getEntry(j); J[index - 1][2] = (float) u.getEntry(i); } /* System.out.println("\nRGB Matrix:"); for(int i=0;i<nClass;i++){ for(int j=0;j<3;j++){ System.out.print(J[i][j]+"\t"); } System.out.println(); } */ return J; }
From source file:net.myrrix.common.math.CommonsMathSolver.java
@Override public double[] solveFToD(float[] b) { ArrayRealVector bVec = new ArrayRealVector(b.length); for (int i = 0; i < b.length; i++) { bVec.setEntry(i, b[i]); }/* ww w . j a v a2 s .co m*/ RealVector vec = solver.solve(bVec); double[] result = new double[b.length]; for (int i = 0; i < result.length; i++) { result[i] = vec.getEntry(i); } return result; }
From source file:lambertmrev.Lambert.java
public ArrayRealVector dTdx(double x, double T, double m_lambda) { double l2 = m_lambda * m_lambda; double l3 = l2 * m_lambda; double umx2 = 1.0 - x * x; double y = FastMath.sqrt(1.0 - l2 * umx2); double y2 = y * y; double y3 = y2 * y; ArrayRealVector out = new ArrayRealVector(3); out.setEntry(0, 1.0 / umx2 * (3.0 * T * x - 2.0 + 2.0 * l3 * x / y)); out.setEntry(1, 1.0 / umx2 * (3.0 * T + 5.0 * x * out.getEntry(0) + 2.0 * (1.0 - l2) * l3 / y3)); out.setEntry(2, 1.0 / umx2//from www .ja va2 s . co m * (7.0 * x * out.getEntry(1) + 8.0 * out.getEntry(0) - 6.0 * (1.0 - l2) * l2 * l3 * x / y3 / y2)); return out; }
From source file:gamlss.algorithm.Gamlss.java
/** This is to emulate the Gamlss algorithm where * the user specifies response variable vector * and design matrix./*from w w w . j av a 2 s. c om*/ * * @param y - vector of response variable values * @param designMatrices - design matrices for each * of the distribution parameters * @param smoothMatrices - smoother matrices for each * of the distribution parameters */ public Gamlss(final ArrayRealVector y, Hashtable<Integer, BlockRealMatrix> designMatrices, final HashMap<Integer, BlockRealMatrix> smoothMatrices) { GAMLSSFamilyDistribution distr = null; switch (DistributionSettings.DISTR) { case DistributionSettings.NO: distr = new NO(); distr.initialiseDistributionParameters(y); break; case DistributionSettings.TF: distr = new TF(); distr.initialiseDistributionParameters(y); break; case DistributionSettings.GA: distr = new GA(); distr.initialiseDistributionParameters(y); break; case DistributionSettings.GT: distr = new GT(); distr.initialiseDistributionParameters(y); break; case DistributionSettings.ST3: distr = new ST3(); distr.initialiseDistributionParameters(y); break; case DistributionSettings.ST4: distr = new ST4(); distr.initialiseDistributionParameters(y); break; case DistributionSettings.JSUo: distr = new JSUo(); distr.initialiseDistributionParameters(y); break; case DistributionSettings.TF2: distr = new TF2(); distr.initialiseDistributionParameters(y); break; case DistributionSettings.SST: distr = new SST(); distr.initialiseDistributionParameters(y); break; case DistributionSettings.BCPE: distr = new BCPE(); distr.initialiseDistributionParameters(y); break; case DistributionSettings.ST1: distr = new ST1(); distr.initialiseDistributionParameters(y); break; case DistributionSettings.PE: distr = new PE(); distr.initialiseDistributionParameters(y); break; default: System.err.println("The specific distribution " + "has not been implemented yet in Gamlss!"); } if (smoothMatrices != null) { Controls.SMOOTHING = true; } if (designMatrices == null) { designMatrices = new Hashtable<Integer, BlockRealMatrix>(); for (int i = 1; i < distr.getNumberOfDistribtionParameters() + 1; i++) { designMatrices.put(i, MatrixFunctions.buildInterceptMatrix(y.getDimension())); Controls.NO_INTERCEPT[i - 1] = true; } } ArrayRealVector w = new ArrayRealVector(y.getDimension()); for (int i = 0; i < w.getDimension(); i++) { w.setEntry(i, Controls.DEFAULT_WEIGHT); } switch (DistributionSettings.FITTING_ALG) { case DistributionSettings.RS: rs = new RSAlgorithm(distr, y, designMatrices, smoothMatrices, w); rs.functionRS(distr, y, w); break; case DistributionSettings.RS20CG20: rs = new RSAlgorithm(distr, y, designMatrices, smoothMatrices, w); rs.functionRS(distr, y, w); rs = null; cg = new CGAlgorithm(y.getDimension()); cg.setnCyc(Controls.GAMLSS_NUM_CYCLES); cg.CGfunction(distr, y, designMatrices, w); break; case DistributionSettings.CG: cg = new CGAlgorithm(y.getDimension()); cg.CGfunction(distr, y, designMatrices, w); break; case DistributionSettings.GO: break; default: System.err.println(" Cannot recognise the " + "fitting algorithm"); } int parNum = distr.getNumberOfDistribtionParameters(); resultsMatrix = new BlockRealMatrix(y.getDimension(), parNum); for (int i = 0; i < y.getDimension(); i++) { for (int j = 0; j < parNum; j++) { resultsMatrix.setEntry(i, j, distr.getDistributionParameter(j + 1).getEntry(i)); } } }
From source file:edu.utexas.cs.tactex.tariffoptimization.TariffOptimizierTOUFixedMargin.java
private ArrayRealVector computeEnergyUnitCosts(CostCurvesPredictor costCurvesPredictor, ArrayRealVector predictedMyCustomersEnergy, ArrayRealVector predictedCompetitorsEnergy, int currentTimeslot) { ArrayRealVector energyUnitCosts = new ArrayRealVector(predictedMyCustomersEnergy.getDimension()); for (int i = 0; i < energyUnitCosts.getDimension(); ++i) { int futureTimeslot = currentTimeslot + i + 1; // +1 since energy record starts next timeslot energyUnitCosts.setEntry(i, costCurvesPredictor.predictUnitCostKwh(currentTimeslot, futureTimeslot, predictedMyCustomersEnergy.getEntry(i), predictedCompetitorsEnergy.getEntry(i)) + costCurvesPredictor.getFudgeFactorKwh(currentTimeslot)); }/*from w w w. j a v a 2 s .c o m*/ //log.debug("energyUnitCosts =" + energyUnitCosts.toString()); return energyUnitCosts; }
From source file:edu.utexas.cs.tactex.BrokerUtilsTest.java
@Test public void test_rotateWeeklyRecordAndAppendTillEndOfDay() { // this test is very similar and is based on // EnergyPredictionTest.testEnergyPredictionOfAboutOneWeek() // we moved it here after refactoring a method to BrokersUtils ////from w ww. j ava 2 s .co m // initialize an record vector where the value of an // entry is its index ArrayRealVector record = new ArrayRealVector(7 * 24); for (int i = 0; i < record.getDimension(); ++i) { record.setEntry(i, i); } int currentTimeslot; ArrayRealVector expected; ArrayRealVector actual; currentTimeslot = 7 * 24 - 1; expected = new ArrayRealVector(record); actual = BrokerUtils.rotateWeeklyRecordAndAppendTillEndOfDay(record, currentTimeslot); assertArrayEquals("no rotation at the beginning of week", expected.toArray(), actual.toArray(), 1e-6); currentTimeslot = 1 * 24 - 1; // actual actual = BrokerUtils.rotateWeeklyRecordAndAppendTillEndOfDay(record, currentTimeslot); // prepare expected RealVector slice1 = record.getSubVector(1 * 24, 7 * 24 - 24); expected.setSubVector(0, slice1); expected.setSubVector(slice1.getDimension(), record.getSubVector(0, 24)); assertArrayEquals("end of first day results in days 2,...,7,1", expected.toArray(), actual.toArray(), 1e-6); currentTimeslot = 6 * 24 - 1; // actual actual = BrokerUtils.rotateWeeklyRecordAndAppendTillEndOfDay(record, currentTimeslot); // prepare expected slice1 = record.getSubVector(6 * 24, 7 * 24 - 6 * 24); expected.setSubVector(0, slice1); expected.setSubVector(slice1.getDimension(), record.getSubVector(0, 6 * 24)); assertArrayEquals("end of 6th day results in days 7,1,...,6", expected.toArray(), actual.toArray(), 1e-6); currentTimeslot = 0; // predict actual = BrokerUtils.rotateWeeklyRecordAndAppendTillEndOfDay(record, currentTimeslot); // prepare expected slice1 = record.getSubVector(1, 7 * 24 - 1); expected.setSubVector(0, slice1); expected.setSubVector(slice1.getDimension(), record.getSubVector(0, 1)); expected.append(record.getSubVector(1, 24 - 1)); assertArrayEquals("if not at day start, appending until the end of day", expected.toArray(), actual.toArray(), 1e-6); }
From source file:edu.utexas.cs.tactex.LWRCustomerPredictionTests.java
/** * testing LWR predictions - // w ww .j a v a2 s .c o m * based on matlab code for LWR * */ @Test public void testLWR() { // Ejml, with intercept SimpleMatrix X; SimpleMatrix y; SimpleMatrix x0; X = new SimpleMatrix(3, 2); y = new SimpleMatrix(3, 1); x0 = new SimpleMatrix(2, 1); // (with intercept) // initialize examples X.setColumn(0, 0, // X-column 1: (intercept) 1, 1, 1); X.setColumn(1, 0, // X-column 2: -1, 0, 1); y.setColumn(0, 0, // Y-column 3, 2, 4); x0.set(1); // initialize to 1s (will override non-intercept terms later) double tau = 0.5; // testing based on matlab results x0.set(1, 0, -1.5); double expected = 3.43855279053977; double actual = lwrCustNewEjml.LWRPredict(X, y, x0, tau); assertEquals("Matlab based expected value 1", expected, actual, 1e-6); x0.set(1, 0, -0.5); expected = 2.62107459906727; actual = lwrCustNewEjml.LWRPredict(X, y, x0, tau); assertEquals("Matlab based expected value 2", expected, actual, 1e-6); x0.set(1, 0, 0); expected = 2.63582467285126; actual = lwrCustNewEjml.LWRPredict(X, y, x0, tau); assertEquals("Matlab based expected value 3", expected, actual, 1e-6); x0.set(1, 0, 0.5); expected = 3.12107459906727; actual = lwrCustNewEjml.LWRPredict(X, y, x0, tau); assertEquals("Matlab based expected value 4", expected, actual, 1e-6); x0.set(1, 0, 1.5); expected = 4.93855279053977; actual = lwrCustNewEjml.LWRPredict(X, y, x0, tau); assertEquals("Matlab based expected value 5", expected, actual, 1e-6); // No intercept - should be same pred. between Ejml and Appache // build Ejml matrix with no intercept SimpleMatrix Xejml = X.extractMatrix(0, X.numRows(), 1, X.numCols()); // build appache matrix with no intercept ArrayRealVector Xappache = new ArrayRealVector(Xejml.numRows()); for (int i = 0; i < Xejml.numRows(); ++i) { Xappache.setEntry(i, Xejml.get(i, 0)); } ArrayRealVector yappache = new ArrayRealVector(y.numRows()); for (int i = 0; i < y.numRows(); ++i) { yappache.setEntry(i, y.get(i, 0)); } // initialize x0 x0 = new SimpleMatrix(1, 1); double x0val; double actualEjml; double actualAppache; x0val = -1.5; x0.set(0, 0, x0val); expected = 4.47403745685534; actualEjml = lwrCustNewEjml.LWRPredict(Xejml, y, x0, tau); actualAppache = lwrCustOldAppache.LWRPredict(Xappache, yappache, x0val, tau); assertEquals("NoIcpt - Matlab based - value 1, Ejml", expected, actualEjml, 1e-6); assertEquals("NoIcpt - Matlab based - value 1, Appache", expected, actualAppache, 1e-6); x0val = -0.5; x0.set(0, 0, x0val); expected = 1.08278977292259; actualEjml = lwrCustNewEjml.LWRPredict(Xejml, y, x0, tau); actualAppache = lwrCustOldAppache.LWRPredict(Xappache, yappache, x0val, tau); assertEquals("NoIcpt - Matlab based - value 2, Ejml", expected, actualEjml, 1e-6); assertEquals("NoIcpt - Matlab based - value 2, Appache", expected, actualAppache, 1e-6); x0val = 0; x0.set(0, 0, x0val); expected = 0; actualEjml = lwrCustNewEjml.LWRPredict(Xejml, y, x0, tau); actualAppache = lwrCustOldAppache.LWRPredict(Xappache, yappache, x0val, tau); assertEquals("NoIcpt - Matlab based - value 3, Ejml", expected, actualEjml, 1e-6); assertEquals("NoIcpt - Matlab based - value 3, Appache", expected, actualAppache, 1e-6); x0val = 0.5; x0.set(0, 0, x0val); expected = 1.58278977292259; actualEjml = lwrCustNewEjml.LWRPredict(Xejml, y, x0, tau); actualAppache = lwrCustOldAppache.LWRPredict(Xappache, yappache, x0val, tau); assertEquals("NoIcpt - Matlab based - value 4, Ejml", expected, actualEjml, 1e-6); assertEquals("NoIcpt - Matlab based - value 4, Appache", expected, actualAppache, 1e-6); x0val = 1.5; x0.set(0, 0, x0val); expected = 5.97403745685534; actualEjml = lwrCustNewEjml.LWRPredict(Xejml, y, x0, tau); actualAppache = lwrCustOldAppache.LWRPredict(Xappache, yappache, x0val, tau); assertEquals("NoIcpt - Matlab based - value 5, Ejml", expected, actualEjml, 1e-6); assertEquals("NoIcpt - Matlab based - value 5, Appache", expected, actualAppache, 1e-6); }
From source file:lambertmrev.Lambert.java
/** Constructs and solves a Lambert problem. * * \param[in] R1 first cartesian position * \param[in] R2 second cartesian position * \param[in] tof time of flight//w w w . ja v a 2 s . co m * \param[in] mu gravity parameter * \param[in] cw when 1 a retrograde orbit is assumed * \param[in] multi_revs maximum number of multirevolutions to compute */ public void lambert_problem(Vector3D r1, Vector3D r2, double tof, double mu, Boolean cw, int multi_revs) { // sanity checks if (tof <= 0) { System.out.println("ToF is negative! \n"); } if (mu <= 0) { System.out.println("mu is below zero"); } // 1 - getting lambda and T double m_c = FastMath.sqrt((r2.getX() - r1.getX()) * (r2.getX() - r1.getX()) + (r2.getY() - r1.getY()) * (r2.getY() - r1.getY()) + (r2.getZ() - r1.getZ()) * (r2.getZ() - r1.getZ())); double R1 = r1.getNorm(); double R2 = r2.getNorm(); double m_s = (m_c + R1 + R2) / 2.0; Vector3D ir1 = r1.normalize(); Vector3D ir2 = r2.normalize(); Vector3D ih = Vector3D.crossProduct(ir1, ir2); ih = ih.normalize(); if (ih.getZ() == 0) { System.out.println("angular momentum vector has no z component \n"); } double lambda2 = 1.0 - m_c / m_s; double m_lambda = FastMath.sqrt(lambda2); Vector3D it1 = new Vector3D(0.0, 0.0, 0.0); Vector3D it2 = new Vector3D(0.0, 0.0, 0.0); if (ih.getZ() < 0.0) { // Transfer angle is larger than 180 degrees as seen from abive the z axis m_lambda = -m_lambda; it1 = Vector3D.crossProduct(ir1, ih); it2 = Vector3D.crossProduct(ir2, ih); } else { it1 = Vector3D.crossProduct(ih, ir1); it2 = Vector3D.crossProduct(ih, ir2); } it1.normalize(); it2.normalize(); if (cw) { // Retrograde motion m_lambda = -m_lambda; it1.negate(); it2.negate(); } double lambda3 = m_lambda * lambda2; double T = FastMath.sqrt(2.0 * mu / m_s / m_s / m_s) * tof; // 2 - We now hava lambda, T and we will find all x // 2.1 - let us first detect the maximum number of revolutions for which there exists a solution int m_Nmax = FastMath.toIntExact(FastMath.round(T / FastMath.PI)); double T00 = FastMath.acos(m_lambda) + m_lambda * FastMath.sqrt(1.0 - lambda2); double T0 = (T00 + m_Nmax * FastMath.PI); double T1 = 2.0 / 3.0 * (1.0 - lambda3); double DT = 0.0; double DDT = 0.0; double DDDT = 0.0; if (m_Nmax > 0) { if (T < T0) { // We use Halley iterations to find xM and TM int it = 0; double err = 1.0; double T_min = T0; double x_old = 0.0, x_new = 0.0; while (true) { ArrayRealVector deriv = dTdx(x_old, T_min, m_lambda); DT = deriv.getEntry(0); DDT = deriv.getEntry(1); DDDT = deriv.getEntry(2); if (DT != 0.0) { x_new = x_old - DT * DDT / (DDT * DDT - DT * DDDT / 2.0); } err = FastMath.abs(x_old - x_new); if ((err < 1e-13) || (it > 12)) { break; } tof = x2tof(x_new, m_Nmax, m_lambda); x_old = x_new; it++; } if (T_min > T) { m_Nmax -= 1; } } } // We exit this if clause with Mmax being the maximum number of revolutions // for which there exists a solution. We crop it to multi_revs m_Nmax = FastMath.min(multi_revs, m_Nmax); // 2.2 We now allocate the memory for the output variables m_v1 = MatrixUtils.createRealMatrix(m_Nmax * 2 + 1, 3); RealMatrix m_v2 = MatrixUtils.createRealMatrix(m_Nmax * 2 + 1, 3); RealMatrix m_iters = MatrixUtils.createRealMatrix(m_Nmax * 2 + 1, 3); //RealMatrix m_x = MatrixUtils.createRealMatrix(m_Nmax*2+1, 3); ArrayRealVector m_x = new ArrayRealVector(m_Nmax * 2 + 1); // 3 - We may now find all solution in x,y // 3.1 0 rev solution // 3.1.1 initial guess if (T >= T00) { m_x.setEntry(0, -(T - T00) / (T - T00 + 4)); } else if (T <= T1) { m_x.setEntry(0, T1 * (T1 - T) / (2.0 / 5.0 * (1 - lambda2 * lambda3) * T) + 1); } else { m_x.setEntry(0, FastMath.pow((T / T00), 0.69314718055994529 / FastMath.log(T1 / T00)) - 1.0); } // 3.1.2 Householder iterations //m_iters.setEntry(0, 0, housOutTmp.getEntry(0)); m_x.setEntry(0, householder(T, m_x.getEntry(0), 0, 1e-5, 15, m_lambda)); // 3.2 multi rev solutions double tmp; double x0; for (int i = 1; i < m_Nmax + 1; i++) { // 3.2.1 left householder iterations tmp = FastMath.pow((i * FastMath.PI + FastMath.PI) / (8.0 * T), 2.0 / 3.0); m_x.setEntry(2 * i - 1, (tmp - 1) / (tmp + 1)); x0 = householder(T, m_x.getEntry(2 * i - 1), i, 1e-8, 15, m_lambda); m_x.setEntry(2 * i - 1, x0); //m_iters.setEntry(2*i-1, 0, housOutTmp.getEntry(0)); //3.2.1 right Householder iterations tmp = FastMath.pow((8.0 * T) / (i * FastMath.PI), 2.0 / 3.0); m_x.setEntry(2 * i, (tmp - 1) / (tmp + 1)); x0 = householder(T, m_x.getEntry(2 * i), i, 1e-8, 15, m_lambda); m_x.setEntry(2 * i, x0); //m_iters.setEntry(2*i, 0, housOutTmp.getEntry(0)); } // 4 - For each found x value we recontruct the terminal velocities double gamma = FastMath.sqrt(mu * m_s / 2.0); double rho = (R1 - R2) / m_c; double sigma = FastMath.sqrt(1 - rho * rho); double vr1, vt1, vr2, vt2, y; ArrayRealVector ir1_vec = new ArrayRealVector(3); ArrayRealVector ir2_vec = new ArrayRealVector(3); ArrayRealVector it1_vec = new ArrayRealVector(3); ArrayRealVector it2_vec = new ArrayRealVector(3); // set Vector3D values to a mutable type ir1_vec.setEntry(0, ir1.getX()); ir1_vec.setEntry(1, ir1.getY()); ir1_vec.setEntry(2, ir1.getZ()); ir2_vec.setEntry(0, ir2.getX()); ir2_vec.setEntry(1, ir2.getY()); ir2_vec.setEntry(2, ir2.getZ()); it1_vec.setEntry(0, it1.getX()); it1_vec.setEntry(1, it1.getY()); it1_vec.setEntry(2, it1.getZ()); it2_vec.setEntry(0, it2.getX()); it2_vec.setEntry(1, it2.getY()); it2_vec.setEntry(2, it2.getZ()); for (int i = 0; i < m_x.getDimension(); i++) { y = FastMath.sqrt(1.0 - lambda2 + lambda2 * m_x.getEntry(i) * m_x.getEntry(i)); vr1 = gamma * ((m_lambda * y - m_x.getEntry(i)) - rho * (m_lambda * y + m_x.getEntry(i))) / R1; vr2 = -gamma * ((m_lambda * y - m_x.getEntry(i)) + rho * (m_lambda * y + m_x.getEntry(i))) / R2; double vt = gamma * sigma * (y + m_lambda * m_x.getEntry(i)); vt1 = vt / R1; vt2 = vt / R2; for (int j = 0; j < 3; ++j) m_v1.setEntry(i, j, vr1 * ir1_vec.getEntry(j) + vt1 * it1_vec.getEntry(j)); for (int j = 0; j < 3; ++j) m_v2.setEntry(i, j, vr2 * ir2_vec.getEntry(j) + vt2 * it2_vec.getEntry(j)); } }
From source file:edu.stanford.cfuller.imageanalysistools.filter.MaximumSeparabilityThresholdingFilter.java
/** * Applies the filter to an Image, optionally turning on adaptive determination of the increment size used to find the threshold, and specifying a size for the threshold determination increment. * Turning on adaptive determination of the increment will generally make the threshold slightly less optimal, but can sometimes speed up the filtering, especially * for images with a large dynamic range. * <p>// w w w . j av a2 s . c o m * The increment size specified (in greylevels) will be used to determine the threshold only if adaptive determination is off; otherwise this parameter will be ignored. * @param im The Image to be thresholded; this will be overwritten by the thresholded Image. * @param adaptiveincrement true to turn on adaptive determination of the threshold increment; false to turn it off and use the default value * @param increment the increment size (in greylevels) to use for determining the threshold; must be positive. */ public void apply_ext(WritableImage im, boolean adaptiveincrement, int increment) { Histogram h = new Histogram(im); int thresholdValue = 0; final int numSteps = 1000; double best_eta = 0; int best_index = Integer.MAX_VALUE; int nonzerocounts = h.getTotalCounts() - h.getCounts(0); double meannonzero = h.getMeanNonzero(); ArrayRealVector omega_v = new ArrayRealVector(h.getMaxValue()); ArrayRealVector mu_v = new ArrayRealVector(h.getMaxValue()); int c = 0; if (adaptiveincrement) { increment = (int) ((h.getMaxValue() - h.getMinValue() + 1) * 1.0 / numSteps); if (increment < 1) increment = 1; } for (int k = h.getMinValue(); k < h.getMaxValue() + 1; k += increment) { if (k == 0) continue; omega_v.setEntry(c, h.getCumulativeCounts(k) * 1.0 / nonzerocounts); if (c == 0) { mu_v.setEntry(c, k * omega_v.getEntry(c)); } else { mu_v.setEntry(c, mu_v.getEntry(c - 1) + k * h.getCounts(k) * 1.0 / nonzerocounts); for (int i = k - increment + 1; i < k; i++) { mu_v.setEntry(c, mu_v.getEntry(c) + h.getCounts(i) * i * 1.0 / nonzerocounts); } } double omega = omega_v.getEntry(c); double mu = mu_v.getEntry(c); if (omega > 1e-8 && 1 - omega > 1e-8) { double eta = omega * (1 - omega) * Math.pow((meannonzero - mu) / (1 - omega) - mu / omega, 2); if (eta >= best_eta) { best_eta = eta; best_index = k; } //System.out.printf("%d, %f\n", k, eta); } c++; } thresholdValue = best_index; if (thresholdValue == Integer.MAX_VALUE) { thresholdValue = 0; } for (ImageCoordinate coord : im) { if (im.getValue(coord) < thresholdValue) im.setValue(coord, 0); } }