List of usage examples for org.apache.commons.math3.linear RealMatrix getEntry
double getEntry(int row, int column) throws OutOfRangeException;
From source file:put.ci.cevo.framework.algorithms.ApacheCMAES.java
/** * @param mat Input matrix./* ww w . ja va 2 s .com*/ * @param n Number of row replicates. * @param m Number of column replicates. * @return a matrix which replicates the input matrix in both directions. */ private static RealMatrix repmat(final RealMatrix mat, int n, int m) { final int rd = mat.getRowDimension(); final int cd = mat.getColumnDimension(); final double[][] d = new double[n * rd][m * cd]; for (int r = 0; r < n * rd; r++) { for (int c = 0; c < m * cd; c++) { d[r][c] = mat.getEntry(r % rd, c % cd); } } return new Array2DRowRealMatrix(d, false); }
From source file:put.ci.cevo.framework.algorithms.ApacheCMAES.java
/** * @param m Input matrix.//from ww w. ja v a2 s . c om * @return the maximum of the matrix element values. */ private static double max(final RealMatrix m) { double max = -Double.MAX_VALUE; for (int r = 0; r < m.getRowDimension(); r++) { for (int c = 0; c < m.getColumnDimension(); c++) { double e = m.getEntry(r, c); if (max < e) { max = e; } } } return max; }
From source file:put.ci.cevo.framework.algorithms.ApacheCMAES.java
/** * @param m Input matrix.// w w w. j a va 2 s . c o m * @return the minimum of the matrix element values. */ private static double min(final RealMatrix m) { double min = Double.MAX_VALUE; for (int r = 0; r < m.getRowDimension(); r++) { for (int c = 0; c < m.getColumnDimension(); c++) { double e = m.getEntry(r, c); if (min > e) { min = e; } } } return min; }
From source file:RBFNN.MainFrame.java
public static int[] compare(RealMatrix YactM, RealMatrix YexpM) { int[] confusion = new int[4]; int len = YactM.getRowDimension(); for (int i = 0; i < len; i++) {//classify the result into the category if (YactM.getEntry(i, 0) > 0.5) YactM.setEntry(i, 0, 1);/*from w w w . j a v a 2 s . co m*/ else YactM.setEntry(i, 0, 0); if (YactM.getEntry(i, 0) == 1 && YexpM.getEntry(i, 0) == 1) confusion[0]++; if (YactM.getEntry(i, 0) == 1 && YexpM.getEntry(i, 0) == 0) confusion[1]++; if (YactM.getEntry(i, 0) == 0 && YexpM.getEntry(i, 0) == 1) confusion[2]++; if (YactM.getEntry(i, 0) == 0 && YexpM.getEntry(i, 0) == 0) confusion[3]++; } return confusion; }
From source file:restclient.service.DailyRecordFacadeREST.java
@GET @Path("findCorrelationByStartAndEndDateAndWeatherVariable/{startDate}/{endDate}/{attribute}") @Produces({ "application/json" }) public String findCorrelationByStartAndEndDateAndWeatherVariable(@PathParam("startDate") Date startDate, @PathParam("endDate") Date endDate, @PathParam("attribute") String attribute) { String attributeInLowerCase = attribute.toLowerCase(); switch (attributeInLowerCase) { case "windspeed": case "wind speed": attributeInLowerCase = "windSpeed"; break;/*w w w . j a v a 2 s . c o m*/ case "atmosphericpressure": case "atmospheric pressure": attributeInLowerCase = "atmosphericPressure"; break; default: break; } String sth = "NEW restclient.Result2(d.painLevel, d." + attributeInLowerCase + ")"; String jpql = "SELECT " + sth + " FROM restclient.DailyRecord d WHERE d.recordDate >= :startDate AND d.recordDate <= :endDate"; TypedQuery<Result2> q = em.createQuery(jpql, Result2.class); q.setParameter("startDate", startDate); q.setParameter("endDate", endDate); List<Result2> result = q.getResultList(); double data[][] = new double[result.size()][]; for (int i = 0; i < result.size(); i++) { data[i] = new double[] { result.get(i).painLevel, result.get(i).weather }; } RealMatrix m = MatrixUtils.createRealMatrix(data); String first = ""; for (int i = 0; i < m.getColumnDimension(); i++) for (int j = 0; j < m.getColumnDimension(); j++) { PearsonsCorrelation pc = new PearsonsCorrelation(); double cor = pc.correlation(m.getColumn(i), m.getColumn(j)); first += (i + "," + j + "=[" + String.format(".%2f", cor) + "," + "]" + "; "); } PearsonsCorrelation pc = new PearsonsCorrelation(m); RealMatrix corM = pc.getCorrelationMatrix(); String second = ("!correlation:" + corM.getEntry(0, 1) + " "); RealMatrix pM = pc.getCorrelationPValues(); String third = ("!p value:" + pM.getEntry(0, 1)); return first + second + third; }
From source file:restclient.service.RecordFacadeREST.java
@GET @Path("findCorrelation/{uid}/{sdate}/{edate}/{wvariable}") @Produces({ "application/json" }) public List<Correlation> findCorrelation(@PathParam("uid") Integer uid, @PathParam("sdate") String date1, @PathParam("edate") String date2, @PathParam("wvariable") String wv) throws ParseException { SimpleDateFormat sdf1 = new SimpleDateFormat("yyyy-MM-dd"); SimpleDateFormat sdf2 = new SimpleDateFormat("dd/MM/yyyy"); Date sdate = sdf1.parse(date1); Date edate = sdf1.parse(date2); TypedQuery<Record> q = em.createQuery( "SELECT r FROM Record r WHERE r.date >= :sdate AND r.date <= :edate AND r.uid.uid = :uid order by r.date ASC", Record.class); q.setParameter("uid", uid); q.setParameter("sdate", sdate); q.setParameter("edate", edate); List<Record> qr = q.getResultList(); List<Correlation> re = new ArrayList<Correlation>(); Correlation cl = new Correlation(); double data[][] = new double[qr.size()][2]; for (int i = 0; i < qr.size(); ++i) { Record r = qr.get(i);/* w w w . ja va2 s . c o m*/ data[i][0] = r.getPlevel(); if (wv.equals("temperature")) { data[i][1] = r.getTemp(); } else if (wv.equals("humidity")) { data[i][1] = r.getHumidity(); } else if (wv.equals("windspeed")) { data[i][1] = r.getWindspeed(); } else { data[i][1] = r.getPressure(); } } RealMatrix m = MatrixUtils.createRealMatrix(data); PearsonsCorrelation pc = new PearsonsCorrelation(m); RealMatrix corM = pc.getCorrelationMatrix(); cl.setRvalue(corM.getEntry(0, 1)); RealMatrix pM = pc.getCorrelationPValues(); cl.setSvalue(pM.getEntry(0, 1)); re.add(cl); return re; }
From source file:Rotationforest.rotationalgo.java
public static void main(String[] args) { // TODO code application logic here double[][] a = new double[303][7]; double[][] b = new double[303][7]; double[][] c = new double[225][7]; double[][] d = new double[225][7]; double sum[] = new double[7]; double mean[]; double[][] inputs = new double[303][14]; final RealMatrix cov; final RealMatrix cov1; read obj = new read(); ArrayList<values> ar = new ArrayList<values>(); ar = obj.run();//from w w w . j a v a2 s. com // System.out.println("records= " + ar.size()); Iterator<values> arIterator = ar.iterator(); int i = 0; while (arIterator.hasNext()) { values values = arIterator.next(); //inputs= new double[][]{{values.age,values.sex,values.cp,values.restbps,values.chol,values.fbs,values.restecg,values.thalach,values.exang,values.oldpeak,values.slope,values.ca,values.thal,values.clas }}; inputs[i][0] = values.age; inputs[i][1] = values.sex; inputs[i][2] = values.cp; inputs[i][3] = values.restbps; inputs[i][4] = values.chol; inputs[i][5] = values.fbs; inputs[i][6] = values.restecg; inputs[i][7] = values.thalach; inputs[i][8] = values.exang; inputs[i][9] = values.oldpeak; inputs[i][10] = values.slope; inputs[i][11] = values.ca; inputs[i][12] = values.thal; inputs[i][13] = values.clas; i++; /*for (int m = 0; m <inputs.length; m++) { for (int n = 0; n <inputs[0].length; n++) { //System.out.print(inputs[m][n] + " "); } //System.out.print("\n\n\n"); } */ } /* final double [][]in =new double[303][14]; for (int r=0; r<inputs.length; r++) { for (int c=0; c<inputs[r].length; c++) { in[r][c]=inputs[r][c]; } }*/ for (int m = 0; m < inputs.length; m++) { for (int n = 0; n <= 6; n++) { a[m][n] = inputs[m][n]; } } // System.out.println(" main inputs = " + inputs); // Corresponding outputs, xor training data /* for (int m = 0; m <inputs.length; m++) { for (int n = 0; n <=6; n++) { System.out.print(a[m][n] + " "); } System.out.print("\n\n\n"); } */ int k = 7; /* for (int m = 0; m <inputs.length; m++) { for (int n = 0; n <=6 && k<=13; n++,k++) { b[m][n]=inputs[m][k]; System.out.print(inputs[m][k] + ""); } System.out.println("\n"); } */ for (int m = 0; m < inputs.length; m++) { for (int n = 0; n < 7; n++) { b[m][n] = inputs[m][k]; k++; } k = 7; } /* for (int m = 0; m <b.length; m++) { for (int n = 0; n <=6; n++) { System.out.print(b[m][n] + " "); } System.out.print("\n\n"); } */ for (int m = 0; m < 225; m++) { for (int n = 0; n <= 6; n++) { c[m][n] = a[m][n]; } } //printing c /* for (int m = 0; m <225; m++) { for (int n = 0; n <=6; n++) { System.out.print(c[m][n] + " "); } System.out.print("\n\n"); } */ for (int m = 0; m < 225; m++) { for (int n = 0; n <= 6; n++) { d[m][n] = b[m][n]; } } Mean ob = new Mean(); PCAMEAN ob1 = new PCAMEAN(); Covariance b1 = new Covariance(); double[] meanc = ob.calmean(c); double[] meand = ob.calmean(d); double[][] adc = ob1.getMeanAdjusted(c, meanc); double[][] add = ob1.getMeanAdjusted(d, meand); cov = b1.computeCovarianceMatrix(adc); cov1 = b1.computeCovarianceMatrix(add); // double[][]digi=ob1.covar(adc); int p = cov.getRowDimension(); int r = cov.getColumnDimension(); double[][] cc = new double[p][r]; double[][] cd = new double[p][r]; //print realmatrix for (int m = 0; m < p; m++) { for (int n = 0; n < r; n++) { cc[m][n] = cov.getEntry(m, n); } } for (int m = 0; m < p; m++) { for (int n = 0; n < r; n++) { cd[m][n] = cov1.getEntry(m, n); } } eigencal e = new eigencal(); double[][] evectorc = e.eigenv(cc); double[][] evectord = e.eigenv(cd); int r1 = evectord.length + evectorc.length; int r2 = evectord[0].length + evectorc[0].length; rot = new double[r1][r2]; for (int m = 0; m < r1; m++) { for (int n = 0; n < r2; n++) { rot[m][n] = 0; } } for (int m = 0; m < evectorc.length; m++) { for (int n = 0; n < evectorc[0].length; n++) { rot[m][n] = evectorc[m][n]; } for (int m1 = 7, m3 = 0; m < evectord.length; m1++, m3++) { for (int n = 7, m4 = 0; n < evectord[0].length; n++, m4++) { rot[m][n] = evectord[m3][m4]; } } } for (int m = 0; m < r1; m++) { for (int n = 0; n < r2; n++) { System.out.print(rot[m][n] + " "); } System.out.print("\n\n"); } System.out.println("Splited the data set into two arrays a & b"); System.out.println("Readed" + "\t" + c.length + "\t" + "records and 7 attributes in array a"); System.out.println("Readed" + "\t" + d.length + "\t" + "records and 7 attributes in array b"); }
From source file:scorePairing.ExtendPairing.java
private ResultsFromEvaluateCost generatePairingAndNullSpaceUsingMatchingPointWithPropertiesNostaticlist( ResultsFromEvaluateCost resultSeed, final CollectionOfPointsWithPropertiesIfc shape1, final CollectionOfPointsWithPropertiesIfc shape2, AlgoParameters algoParameters) { // Both points has StrikingProperties not none and they all match (also if only one of course) Map<Integer, List<PairPointWithDistance>> mapforAllMatchingStrikingProperties = new LinkedHashMap<>(); // Both points has more than one striking properties and at least one of them match List<PairPointWithDistance> listPairsMatchingOneAmongOthersStrikingPropertiesWithShortDistance = new ArrayList<>(); // Both points has None Striking Properties List<PairPointWithDistance> listPairsMatchingOnlyNoneStrikingPropertiesWithShortDistance = new ArrayList<>(); // One point has None striking property and the others has striking properties List<PairPointWithDistance> listPairsNoneMatchingWithoutStrikingPropertiesWithShortDistance = new ArrayList<>(); // Both points has Striking properties but not a single one is matching List<PairPointWithDistance> listPairsNonMatchingWithoutStrikingPropertiesWithShortDistance = new ArrayList<>(); for (int i = 0; i < shape1.getSize(); i++) { for (int j = 0; j < shape2.getSize(); j++) { // Only Pairs of matching properties point are kept PointWithPropertiesIfc pointWithProperties1 = shape1.getPointFromId(i); PointWithPropertiesIfc pointWithProperties2 = shape2.getPointFromId(j); RealMatrix matrix = resultSeed.getRotationMatrix(); float[] vector1 = new float[3]; for (int k = 0; k < 3; k++) { vector1[k] = (float) (pointWithProperties2.getCoords().getCoords()[k] + resultSeed.getTranslationVectorToTranslateShape2ToOrigin().getEntry(k)); }//from w ww. ja v a2 s . c om double[] vector2 = new double[3]; for (int k = 0; k < 3; k++) { for (int l = 0; l < 3; l++) { vector2[k] += matrix.getEntry(k, l) * vector1[l]; } } for (int k = 0; k < 3; k++) { vector1[k] = (float) (vector2[k] - resultSeed.getTranslationVectorToTranslateShape2ToOrigin().getEntry(k) + resultSeed.getTranslationVector().getEntry(k)); } double distance = MathTools.computeDistance(pointWithProperties1.getCoords().getCoords(), vector1); if (distance < algoParameters.getDISTANCE_MIN_FOR_EXTENDED_PAIRING_FROM_SEED()) { PairPointWithDistance pairPointIdsFromMinishape1andMinishape2WithMatchingStrikingPropertiesWithDistance = new PairPointWithDistance( i, j, distance); int countOfMatchingStrikingPropertiesWhenAllAreMatching = StrikingPropertiesTools .evaluatePointsMatchingAllNotNoneProperties(pointWithProperties1, pointWithProperties2); if (countOfMatchingStrikingPropertiesWhenAllAreMatching > 0) { AddToMap.addElementToAMapOfList(mapforAllMatchingStrikingProperties, countOfMatchingStrikingPropertiesWhenAllAreMatching, pairPointIdsFromMinishape1andMinishape2WithMatchingStrikingPropertiesWithDistance); continue; } if (StrikingPropertiesTools.evaluatePointsMatchingOneButNotAllNotNoneProperties( pointWithProperties1, pointWithProperties2)) { listPairsMatchingOneAmongOthersStrikingPropertiesWithShortDistance.add( pairPointIdsFromMinishape1andMinishape2WithMatchingStrikingPropertiesWithDistance); continue; } if (StrikingPropertiesTools.evaluatePointsMatchingOnlyNoneProperties(pointWithProperties1, pointWithProperties2)) { listPairsMatchingOnlyNoneStrikingPropertiesWithShortDistance.add( pairPointIdsFromMinishape1andMinishape2WithMatchingStrikingPropertiesWithDistance); continue; } if (StrikingPropertiesTools.evaluatePointsNoneMatchingANotNONEtoOnlyNoneProperties( pointWithProperties1, pointWithProperties2)) { listPairsNoneMatchingWithoutStrikingPropertiesWithShortDistance.add( pairPointIdsFromMinishape1andMinishape2WithMatchingStrikingPropertiesWithDistance); continue; } if (StrikingPropertiesTools.evaluatePointsNoneMatchingANotNONEtoNotNoneProperties( pointWithProperties1, pointWithProperties2)) { listPairsNonMatchingWithoutStrikingPropertiesWithShortDistance.add( pairPointIdsFromMinishape1andMinishape2WithMatchingStrikingPropertiesWithDistance); continue; } } } } // for (Entry<Integer, List<PairPointWithDistance>> entry: mapforAllMatchingStrikingProperties.entrySet()){ // cleanListPairsFromPointsInvolvedInSeedPairings(resultSeed.getPairingAndNullSpaces(), entry.getValue()); // } // cleanListPairsFromPointsInvolvedInSeedPairings(resultSeed.getPairingAndNullSpaces(), listPairsMatchingOneAmongOthersStrikingPropertiesWithShortDistance); // cleanListPairsFromPointsInvolvedInSeedPairings(resultSeed.getPairingAndNullSpaces(), listPairsMatchingOnlyNoneStrikingPropertiesWithShortDistance); // cleanListPairsFromPointsInvolvedInSeedPairings(resultSeed.getPairingAndNullSpaces(), listPairsNoneMatchingWithoutStrikingPropertiesWithShortDistance); // cleanListPairsFromPointsInvolvedInSeedPairings(resultSeed.getPairingAndNullSpaces(), listPairsNonMatchingWithoutStrikingPropertiesWithShortDistance); for (Map.Entry<Integer, List<PairPointWithDistance>> entry : mapforAllMatchingStrikingProperties .entrySet()) { Collections.sort(entry.getValue(), new PairPointWithDistance.LowestDistancePairPointWithDistance()); } Collections.sort(listPairsMatchingOneAmongOthersStrikingPropertiesWithShortDistance, new PairPointWithDistance.LowestDistancePairPointWithDistance()); Collections.sort(listPairsMatchingOnlyNoneStrikingPropertiesWithShortDistance, new PairPointWithDistance.LowestDistancePairPointWithDistance()); Collections.sort(listPairsNoneMatchingWithoutStrikingPropertiesWithShortDistance, new PairPointWithDistance.LowestDistancePairPointWithDistance()); Collections.sort(listPairsNonMatchingWithoutStrikingPropertiesWithShortDistance, new PairPointWithDistance.LowestDistancePairPointWithDistance()); // loop on the list // when a point already found is found then skip the pair // for (Entry<Integer, List<PairPointWithDistance>> entry: mapforAllMatchingStrikingProperties.entrySet()){ // cleanOfDuplicatePoints(entry.getValue()); // } // cleanOfDuplicatePoints(listPairsMatchingOneAmongOthersStrikingPropertiesWithShortDistance); // cleanOfDuplicatePoints(listPairsMatchingOnlyNoneStrikingPropertiesWithShortDistance); // cleanOfDuplicatePoints(listPairsNoneMatchingWithoutStrikingPropertiesWithShortDistance); // cleanOfDuplicatePoints(listPairsNonMatchingWithoutStrikingPropertiesWithShortDistance); // Now I should try to pair unpaired point so far // clean the list of pairs matching distance but without matching properties OF pairs having points in matching properties list List<PairPointWithDistance> globalList = new ArrayList<>(); List<Integer> keys = new ArrayList<>(); keys.addAll(mapforAllMatchingStrikingProperties.keySet()); Collections.sort(keys, Collections.reverseOrder()); for (Integer key : keys) { globalList.addAll(mapforAllMatchingStrikingProperties.get(key)); } globalList.addAll(listPairsMatchingOneAmongOthersStrikingPropertiesWithShortDistance); globalList.addAll(listPairsMatchingOnlyNoneStrikingPropertiesWithShortDistance); globalList.addAll(listPairsNoneMatchingWithoutStrikingPropertiesWithShortDistance); globalList.addAll(listPairsNonMatchingWithoutStrikingPropertiesWithShortDistance); cleanOfDuplicatePoints(globalList); // now i can build the pairing // initial pairing might be on minishape and shape1 and aligned shape on a larger set like Shape //PairingAndNullSpaces newPairingAndNewNullSpaces = deepCopyNewPairingAndNewNullSpacesAndExtendIfNeeded(resultSeed.getPairingAndNullSpaces(), shape1, shape2); double newcost = 0.0; double newdistanceResidual = 0.0; RealMatrix newRotationMatrix = resultSeed.getRotationMatrix().copy(); RealVector newtranslationVector = resultSeed.getTranslationVector().copy(); RealVector newtranslationVectorToTranslateShape2ToOrigin = resultSeed .getTranslationVectorToTranslateShape2ToOrigin().copy(); Map<Integer, Integer> newpairing = new LinkedHashMap<>(); List<Integer> newnullSpaceOfMap1 = new ArrayList<>(); List<Integer> newnullSpaceOfMap2 = new ArrayList<>(); PairingAndNullSpaces newPairingAndNullSpaces = new PairingAndNullSpaces(newpairing, newnullSpaceOfMap1, newnullSpaceOfMap2); float newratioPairedPointInQuery = 0.0f; ResultsFromEvaluateCost extendedResult = new ResultsFromEvaluateCost(newcost, newdistanceResidual, newRotationMatrix, newtranslationVector, newtranslationVectorToTranslateShape2ToOrigin, newPairingAndNullSpaces, newratioPairedPointInQuery, algoParameters); fillNullSpaceMap1(shape1, extendedResult.getPairingAndNullSpaces()); fillNullSpaceMap2(shape2, extendedResult.getPairingAndNullSpaces()); // loop on new pairs and add Iterator<PairPointWithDistance> itr = globalList.iterator(); while (itr.hasNext()) { PairPointWithDistance pairPointWithDistance = itr.next(); if (pairPointWithDistance.getDistance() < 0.7) { // We really use the sorting for very short distances if (extendedResult.getPairingAndNullSpaces().getNullSpaceOfMap1() .contains(Integer.valueOf(pairPointWithDistance.getPairInteger().point1)) && extendedResult.getPairingAndNullSpaces().getNullSpaceOfMap2() .contains(Integer.valueOf(pairPointWithDistance.getPairInteger().point2))) { extendedResult.getPairingAndNullSpaces().getPairing().put( pairPointWithDistance.getPairInteger().point1, pairPointWithDistance.getPairInteger().point2); boolean toto = extendedResult.getPairingAndNullSpaces().getNullSpaceOfMap1() .remove(Integer.valueOf(pairPointWithDistance.getPairInteger().point1)); boolean toto2 = extendedResult.getPairingAndNullSpaces().getNullSpaceOfMap2() .remove(Integer.valueOf(pairPointWithDistance.getPairInteger().point2)); itr.remove(); } } } Collections.sort(globalList, new PairPointWithDistance.LowestDistancePairPointWithDistance()); Iterator<PairPointWithDistance> itr2 = globalList.iterator(); while (itr2.hasNext()) { PairPointWithDistance pairPointWithDistance = itr2.next(); if (extendedResult.getPairingAndNullSpaces().getNullSpaceOfMap1() .contains(Integer.valueOf(pairPointWithDistance.getPairInteger().point1)) && extendedResult.getPairingAndNullSpaces().getNullSpaceOfMap2() .contains(Integer.valueOf(pairPointWithDistance.getPairInteger().point2))) { extendedResult.getPairingAndNullSpaces().getPairing().put( pairPointWithDistance.getPairInteger().point1, pairPointWithDistance.getPairInteger().point2); boolean toto = extendedResult.getPairingAndNullSpaces().getNullSpaceOfMap1() .remove(Integer.valueOf(pairPointWithDistance.getPairInteger().point1)); boolean toto2 = extendedResult.getPairingAndNullSpaces().getNullSpaceOfMap2() .remove(Integer.valueOf(pairPointWithDistance.getPairInteger().point2)); } } if (debug == true) { boolean validPairing = PairingTools.validate(extendedResult.getPairingAndNullSpaces()); if (validPairing == false) { System.out.println("Extended pairing is not valid"); } // for debug I recheck if there are remaining short distance List<PairPointWithDistance> remainingShortDistance = new ArrayList<>(); for (Integer point1 : extendedResult.getPairingAndNullSpaces().getNullSpaceOfMap1()) { PointWithPropertiesIfc pointWithProperties1 = shape1.getPointFromId(point1); for (Integer point2 : extendedResult.getPairingAndNullSpaces().getNullSpaceOfMap2()) { PointWithPropertiesIfc pointWithProperties2 = shape2.getPointFromId(point2); RealMatrix matrix = extendedResult.getRotationMatrix(); float[] vector1 = new float[3]; for (int k = 0; k < 3; k++) { vector1[k] = (float) (pointWithProperties2.getCoords().getCoords()[k] + extendedResult.getTranslationVectorToTranslateShape2ToOrigin().getEntry(k)); } double[] vector2 = new double[3]; for (int k = 0; k < 3; k++) { for (int l = 0; l < 3; l++) { vector2[k] += matrix.getEntry(k, l) * vector1[l]; } } for (int k = 0; k < 3; k++) { vector1[k] = (float) (vector2[k] - extendedResult.getTranslationVectorToTranslateShape2ToOrigin().getEntry(k) + extendedResult.getTranslationVector().getEntry(k)); } double distance = MathTools.computeDistance(pointWithProperties1.getCoords().getCoords(), vector1); if (distance < algoParameters.getDISTANCE_MIN_FOR_EXTENDED_PAIRING_FROM_SEED()) { PairPointWithDistance pairPointIdsFromMinishape1andMinishape2WithMatchingStrikingPropertiesWithDistance = new PairPointWithDistance( point1, point2, distance); remainingShortDistance.add( pairPointIdsFromMinishape1andMinishape2WithMatchingStrikingPropertiesWithDistance); //System.out.println(distance); } } } System.out.println("remainingShortDistance = " + remainingShortDistance.size()); if (remainingShortDistance.size() > 10) { System.out.println(); } } return extendedResult; }
From source file:tinfour.gwr.SurfaceGwr.java
/** * Computes the elevation for a point at the specified query * coordinates, by performing a multiple-linear regression using the * observed values. A number of statistics are simultaneously computed * and may be obtained using access methods. The more heavy weight * computations (such as the computation of a variance and covariance * matrix) are deferred until actually requested by the calling application. * <p>/* w w w. j ava 2 s .c o m*/ * <strong>Note:</strong> For efficiency purposes, the arrays for * samples and weights passed to this method are stored in the class * directly. * <p> * The sample weights matrix is a two dimensional array giving * weights based on the distance between samples. It is used when performing * calculations for general statistics such as standard deviation, * confidence intervals, etc. Because of the high cost of initializing this * array, it can be treated as optional in cases where only the regression * coefficients are required. * <p> * A convenience routine for populating the sample weights matrix is * supplied by the initWeightsUsingGaussianKernal method. * * @param model the model to be used for the regression * @param xQuery x coordinate of the query position * @param yQuery y coordinate of the query position * @param nSamples the number of sample points to be used for regression * @param samples an array of dimension [n][3] giving at least nSamples * points with the x, y, and z values for the regression. * @param weights an array of weighting factors for samples * @param sampleWeightsMatrix an optional array of weights based on the * distances between different samples; if general statistics are * not required, pass a null value for this argument. * @return an array of regression coefficients, or null if the * computation failed. */ @SuppressWarnings({ "PMD.ArrayIsStoredDirectly", "PMD.MethodReturnsInternalArray" }) public double[] computeRegression(SurfaceModel model, double xQuery, double yQuery, int nSamples, double[][] samples, double[] weights, double[][] sampleWeightsMatrix) { // clear out previous solutions areVarianceAndHatPrepped = false; this.model = model; this.sigma2 = Double.NaN; this.rss = Double.NaN; this.beta = null; this.hat = null; if (nSamples < model.getCoefficientCount()) { throw new IllegalArgumentException("Insufficient number of samples for regression: found " + nSamples + ", need " + model.getCoefficientCount()); } this.nSamples = nSamples; this.samples = samples; this.weights = weights; this.sampleWeightsMatrix = sampleWeightsMatrix; this.xOffset = xQuery; this.yOffset = yQuery; double[][] g; double[][] input; // In the following expressions, the layout of the regression // variables is organized to simplify the computation of quantities // such as slope and curvature. Recall that the samples are always // mapped so that the query position (xQUery, yQuery) is treated // as the origin. Therefore, then the derivatives are evaluated at // the query position (adjusted origin), many terms drop out and // the following relationships apply // Z = b[0] // Zx = b[1] // partial of z(x,y) with respect to x, etc. // Zy = b[2] // Zxx = 2*b[3] // 2nd partial of z(x,y) with respect to x, etc. // Zyy = 2*b[4] // Zxy = b[5] if (model == SurfaceModel.CubicWithCrossTerms) { nVariables = 9; g = new double[10][1]; input = new double[10][10]; for (int i = 0; i < nSamples; i++) { double x = samples[i][0] - xQuery; double y = samples[i][1] - yQuery; double z = samples[i][2]; double w = weights[i]; double x2 = x * x; double y2 = y * y; double x3 = x * x2; double y3 = y * y2; double x4 = x2 * x2; double y4 = y2 * y2; double xy = x * y; input[0][0] += w; input[0][1] += w * x; input[0][2] += w * y; input[0][3] += w * x2; input[0][4] += w * y2; input[0][5] += w * xy; input[0][6] += w * x2 * y; input[0][7] += w * x * y2; input[0][8] += w * x3; input[0][9] += w * y3; input[1][1] += w * x2; input[1][2] += w * xy; input[1][3] += w * x3; input[1][4] += w * x * y2; input[1][5] += w * x2 * y; input[1][6] += w * x * x2 * y; input[1][7] += w * x * x * y2; input[1][8] += w * x * x3; input[1][9] += w * x * y3; input[2][2] += w * y2; input[2][3] += w * x2 * y; input[2][4] += w * y3; input[2][5] += w * x * y2; input[2][6] += w * y * x2 * y; input[2][7] += w * y * x * y2; input[2][8] += w * y * x3; input[2][9] += w * y * y3; input[3][3] += w * x4; input[3][4] += w * x2 * y2; input[3][5] += w * x3 * y; input[3][6] += w * x2 * x2 * y; input[3][7] += w * x2 * x * y2; input[3][8] += w * x2 * x3; input[3][9] += w * x2 * y3; input[4][4] += w * y4; input[4][5] += w * x * y3; input[4][6] += w * y2 * x2 * y; input[4][7] += w * y2 * x * y2; input[4][8] += w * y2 * x3; input[4][9] += w * y2 * y3; input[5][5] += w * x2 * y2; input[5][6] += w * xy * x2 * y; input[5][7] += w * xy * x * y2; input[5][8] += w * xy * x3; input[5][9] += w * xy * y3; input[6][6] += w * x2 * y * x2 * y; input[6][7] += w * x2 * y * x * y2; input[6][8] += w * x2 * y * x3; input[6][9] += w * x2 * y * y3; input[7][7] += w * y2 * x * x * y2; input[7][8] += w * y2 * x * x3; input[7][9] += w * y2 * x * y3; input[8][8] += w * x3 * x3; input[8][9] += w * x3 * y3; input[9][9] += w * y3 * y3; g[0][0] += w * z; g[1][0] += w * x * z; g[2][0] += w * y * z; g[3][0] += w * x2 * z; g[4][0] += w * y2 * z; g[5][0] += w * xy * z; g[6][0] += w * x2 * y * z; g[7][0] += w * x * y2 * z; g[8][0] += w * x3 * z; g[9][0] += w * y3 * z; } // the input for a least-squares fit is a real-symmetric // matrix. So here the code assigns the symmetric terms. input[1][0] = input[0][1]; input[2][0] = input[0][2]; input[2][1] = input[1][2]; input[3][0] = input[0][3]; input[3][1] = input[1][3]; input[3][2] = input[2][3]; input[4][0] = input[0][4]; input[4][1] = input[1][4]; input[4][2] = input[2][4]; input[4][3] = input[3][4]; input[5][0] = input[0][5]; input[5][1] = input[1][5]; input[5][2] = input[2][5]; input[5][3] = input[3][5]; input[5][4] = input[4][5]; input[6][0] = input[0][6]; input[6][1] = input[1][6]; input[6][2] = input[2][6]; input[6][3] = input[3][6]; input[6][4] = input[4][6]; input[6][5] = input[5][6]; input[7][0] = input[0][7]; input[7][1] = input[1][7]; input[7][2] = input[2][7]; input[7][3] = input[3][7]; input[7][4] = input[4][7]; input[7][5] = input[5][7]; input[7][6] = input[6][7]; input[8][0] = input[0][8]; input[8][1] = input[1][8]; input[8][2] = input[2][8]; input[8][3] = input[3][8]; input[8][4] = input[4][8]; input[8][5] = input[5][8]; input[8][6] = input[6][8]; input[8][7] = input[7][8]; input[9][0] = input[0][9]; input[9][1] = input[1][9]; input[9][2] = input[2][9]; input[9][3] = input[3][9]; input[9][4] = input[4][9]; input[9][5] = input[5][9]; input[9][6] = input[6][9]; input[9][7] = input[7][9]; input[9][8] = input[8][9]; } else if (model == SurfaceModel.QuadraticWithCrossTerms) { // z(x, y) = b0 + b1*x + b2*y +b3*x^2 +b4*y^2+b5*x*y nVariables = 5; g = new double[6][1]; input = new double[6][6]; for (int i = 0; i < nSamples; i++) { double x = samples[i][0] - xQuery; double y = samples[i][1] - yQuery; double z = samples[i][2]; double w = weights[i]; double x2 = x * x; double y2 = y * y; double x3 = x * x2; double y3 = y * y2; double x4 = x2 * x2; double y4 = y2 * y2; double xy = x * y; input[0][0] += w; input[0][1] += w * x; input[0][2] += w * y; input[0][3] += w * x2; input[0][4] += w * y2; input[0][5] += w * xy; input[1][1] += w * x2; input[1][2] += w * xy; input[1][3] += w * x * x2; input[1][4] += w * x * y2; input[1][5] += w * x * xy; input[2][2] += w * y2; input[2][3] += w * x2 * y; input[2][4] += w * y3; input[2][5] += w * x * y2; input[3][3] += w * x4; input[3][4] += w * x2 * y2; input[3][5] += w * x3 * y; input[4][4] += w * y4; input[4][5] += w * x * y3; input[5][5] += w * x2 * y2; g[0][0] += w * z; g[1][0] += w * x * z; g[2][0] += w * y * z; g[3][0] += w * x2 * z; g[4][0] += w * y2 * z; g[5][0] += w * xy * z; } // the input for a least-squares fit is a real-symmetric matrix. // So here the code assigns the symmetric terms. input[1][0] = input[0][1]; input[2][0] = input[0][2]; input[2][1] = input[1][2]; input[3][0] = input[0][3]; input[3][1] = input[1][3]; input[3][2] = input[2][3]; input[4][0] = input[0][4]; input[4][1] = input[1][4]; input[4][2] = input[2][4]; input[4][3] = input[3][4]; input[5][0] = input[0][5]; input[5][1] = input[1][5]; input[5][2] = input[2][5]; input[5][3] = input[3][5]; input[5][4] = input[4][5]; } else if (model == SurfaceModel.Quadratic) { // z(x, y) = b0 + b1*x + b2*y +b3*x^2 +b4*y^2+b5*x*y nVariables = 4; g = new double[5][1]; input = new double[5][5]; for (int i = 0; i < nSamples; i++) { double x = samples[i][0] - xQuery; double y = samples[i][1] - yQuery; double z = samples[i][2]; double w = weights[i]; double x2 = x * x; double y2 = y * y; double x3 = x * x2; double y3 = y * y2; double x4 = x2 * x2; double y4 = y2 * y2; double xy = x * y; input[0][0] += w; input[0][1] += w * x; input[0][2] += w * y; input[0][3] += w * x2; input[0][4] += w * y2; input[1][1] += w * x2; input[1][2] += w * xy; input[1][3] += w * x3; input[1][4] += w * x * y2; input[2][2] += w * y2; input[2][3] += w * x2 * y; input[2][4] += w * y3; input[3][3] += w * x4; input[3][4] += w * x2 * y2; input[4][4] += w * y4; g[0][0] += w * z; g[1][0] += w * x * z; g[2][0] += w * y * z; g[3][0] += w * x2 * z; g[4][0] += w * y2 * z; } // the input for a least-squares fit is a real-symmetric matrix. // So here the code assigns the symmetric terms. input[1][0] = input[0][1]; input[2][0] = input[0][2]; input[2][1] = input[1][2]; input[3][0] = input[0][3]; input[3][1] = input[1][3]; input[3][2] = input[2][3]; input[4][0] = input[0][4]; input[4][1] = input[1][4]; input[4][2] = input[2][4]; input[4][3] = input[3][4]; } else if (model == SurfaceModel.Planar) { // z(x, y) = b0 + b1*x + b2*y; nVariables = 2; g = new double[3][1]; input = new double[3][3]; for (int i = 0; i < nSamples; i++) { double x = samples[i][0] - xQuery; double y = samples[i][1] - yQuery; double z = samples[i][2]; double w = weights[i]; double x2 = x * x; double y2 = y * y; input[0][0] += w; input[0][1] += w * x; input[0][2] += w * y; input[1][1] += w * x2; input[1][2] += w * x * y; input[2][2] += w * y2; g[0][0] += w * z; g[1][0] += w * x * z; g[2][0] += w * y * z; } // the input for a least-squares fit is a real-symmetric matrix. // So here the code assigns the symmetric terms. input[1][0] = input[0][1]; input[2][0] = input[0][2]; input[2][1] = input[1][2]; } else if (model == SurfaceModel.PlanarWithCrossTerms) { // z(x, y) = b0 + b1*x + b2*y + b3*x*y; nVariables = 3; g = new double[4][1]; input = new double[4][4]; for (int i = 0; i < nSamples; i++) { double x = samples[i][0] - xQuery; double y = samples[i][1] - yQuery; double z = samples[i][2]; double w = weights[i]; double x2 = x * x; double y2 = y * y; double xy = x * y; input[0][0] += w; input[0][1] += w * x; input[0][2] += w * y; input[0][3] += w * xy; input[1][1] += w * x2; // x*x input[1][2] += w * xy; // y*x input[1][3] += w * xy * x; // xy*x input[2][2] += w * y2; // y*y input[2][3] += w * xy * y; // xy*y input[3][3] += w * xy * xy; //xy*xy g[0][0] += w * z; g[1][0] += w * x * z; g[2][0] += w * y * z; g[3][0] += w * xy * z; } // the input for a least-squares fit is a real-symmetric matrix. // So here the code assigns the symmetric terms. input[1][0] = input[0][1]; input[2][0] = input[0][2]; input[2][1] = input[1][2]; input[3][0] = input[0][3]; input[3][1] = input[1][3]; input[3][2] = input[2][3]; } else { // if(model==SurfaceModel.Cubic){ nVariables = 6; g = new double[7][1]; input = new double[7][7]; for (int i = 0; i < nSamples; i++) { double x = samples[i][0] - xQuery; double y = samples[i][1] - yQuery; double z = samples[i][2]; double w = weights[i]; double x2 = x * x; double y2 = y * y; double xy = x * y; double x3 = x * x2; double y3 = y * y2; double x4 = x2 * x2; double y4 = y2 * y2; input[0][0] += w; input[0][1] += w * x; input[0][2] += w * y; input[0][3] += w * x2; input[0][4] += w * y2; input[0][5] += w * x3; input[0][6] += w * y3; input[1][1] += w * x2; input[1][2] += w * xy; input[1][3] += w * x3; input[1][4] += w * y2 * x; input[1][5] += w * x4; input[1][6] += w * y3 * x; input[2][2] += w * y2; input[2][3] += w * x2 * y; input[2][4] += w * y3; input[2][5] += w * x3 * y; input[2][6] += w * y4; input[3][3] += w * x4; input[3][4] += w * y2 * x2; input[3][5] += w * x3 * x2; // x5 input[3][6] += w * y3 * x2; input[4][4] += w * y4; input[4][5] += w * x3 * y2; input[4][6] += w * y3 * y2; // y5 input[5][5] += w * x3 * x3; // x6 input[5][6] += w * y3 * x3; input[6][6] += w * y3 * y3; // y6 g[0][0] += w * z; g[1][0] += w * x * z; g[2][0] += w * y * z; g[3][0] += w * x2 * z; g[4][0] += w * y2 * z; g[5][0] += w * x3 * z; g[6][0] += w * y3 * z; } // the input for a least-squares fit is a real-symmetric matrix. // So here the code assigns the symmetric terms. input[1][0] = input[0][1]; input[2][0] = input[0][2]; input[2][1] = input[1][2]; input[3][0] = input[0][3]; input[3][1] = input[1][3]; input[3][2] = input[2][3]; input[4][0] = input[0][4]; input[4][1] = input[1][4]; input[4][2] = input[2][4]; input[4][3] = input[3][4]; input[5][0] = input[0][5]; input[5][1] = input[1][5]; input[5][2] = input[2][5]; input[5][3] = input[3][5]; input[5][4] = input[4][5]; input[6][0] = input[0][6]; input[6][1] = input[1][6]; input[6][2] = input[2][6]; input[6][3] = input[3][6]; input[6][4] = input[4][6]; input[6][5] = input[5][6]; } nDegOfFreedom = nSamples - nVariables - 1; if (nDegOfFreedom < 1) { throw new IllegalArgumentException( "Inadequate sample size " + nSamples + " for " + nDegOfFreedom + " degrees of freedom"); } RealMatrix matrixG = new BlockRealMatrix(g); RealMatrix matrixA = new BlockRealMatrix(input); // The Apache Commons Math MultipleLinearRegression implementation // uses the QRDecomposition, and we follow their lead. // When I first implemented this, I thought that the input matrix would be // a real symmetric and positive-definite matrix. If that were the case, // it could be solved using a Cholesky decomposition. But, even // when evaluating ordinary least squares, I ran into numeric // issues that led to the matrix violating the positive-definite criterion. try { QRDecomposition cd = new QRDecomposition(matrixA); DecompositionSolver solver = cd.getSolver(); RealMatrix solution = solver.solve(matrixG); beta = new double[nVariables + 1]; for (int i = 0; i < beta.length; i++) { beta[i] = solution.getEntry(i, 0); } } catch (SingularMatrixException npex) { return null; } return beta; }
From source file:tinfour.gwr.SurfaceGwr.java
public void computeVarianceAndHat() { if (areVarianceAndHatPrepped) { return;//from w w w .j a va 2 s .c o m } areVarianceAndHatPrepped = true; if (sampleWeightsMatrix == null) { throw new NullPointerException("Null specification for sampleWeightsMatrix"); } else if (sampleWeightsMatrix.length != nSamples) { throw new IllegalArgumentException("Incorrectly specified sampleWeightsMatrix"); } double[][] bigS = new double[nSamples][nSamples]; double[][] bigW = sampleWeightsMatrix; double[][] input = computeDesignMatrix(model, xOffset, yOffset, nSamples, samples); RealMatrix mX = new BlockRealMatrix(input); RealMatrix mXT = mX.transpose(); // in the loop below, we compute // Tr(hat) and Tr(Hat' x Hat) // this second term is actually the square of the // Frobenius Norm. Internally, the Apache Commons classes // may provide a more numerically stable implementation of this operation. // This may be worth future investigation. double sTrace = 0; double sTrace2 = 0; for (int i = 0; i < nSamples; i++) { DiagonalMatrix mW = new DiagonalMatrix(bigW[i]); //NOPMD RealMatrix mXTW = mXT.multiply(mW); RealMatrix rx = mX.getRowMatrix(i); RealMatrix c = mXTW.multiply(mX); QRDecomposition cd = new QRDecomposition(c); // NOPMD DecompositionSolver cdSolver = cd.getSolver(); RealMatrix cInv = cdSolver.getInverse(); RealMatrix r = rx.multiply(cInv).multiply(mXTW); double[] row = r.getRow(0); sTrace += row[i]; System.arraycopy(row, 0, bigS[i], 0, nSamples); for (int j = 0; j < nSamples; j++) { sTrace2 += row[j] * row[j]; } } hat = new BlockRealMatrix(bigS); traceHat = sTrace; traceHat2 = sTrace2; double[][] zArray = new double[nSamples][1]; for (int i = 0; i < nSamples; i++) { zArray[i][0] = samples[i][2]; } RealMatrix mY = new BlockRealMatrix(zArray); RealMatrix mYH = hat.multiply(mY); double sse = 0; for (int i = 0; i < nSamples; i++) { double yHat = mYH.getEntry(i, 0); double e = zArray[i][0] - yHat; sse += e * e; } rss = sse; double d1 = nSamples - (2 * traceHat - sTrace2); sigma2 = rss / d1; mlSigma2 = rss / nSamples; RealMatrix mIL = hat.copy(); for (int i = 0; i < nSamples; i++) { double c = 1.0 - mIL.getEntry(i, i); mIL.setEntry(i, i, c); } RealMatrix mILT = mIL.transpose().multiply(mIL); delta1 = mILT.getTrace(); delta2 = (mILT.multiply(mILT)).getTrace(); }