List of usage examples for org.apache.commons.math3.linear ArrayRealVector dotProduct
@Override public double dotProduct(RealVector v) throws DimensionMismatchException
From source file:EndmemberExtraction.java
public static ArrayList<Integer> ORASIS(double[][] data, int nData, int nDim, double threshold, int[] exemplarLabel) { ArrayRealVector vec;/*from w w w. ja v a 2s . c o m*/ ArrayList<ArrayRealVector> X = new ArrayList<>(); ArrayList<ArrayRealVector> E = new ArrayList<>(); ArrayList<Integer> exemplarIndex = new ArrayList<>(); for (int i = 0; i < nData; i++) { vec = new ArrayRealVector(data[i]); vec.unitize(); X.add(vec); } E.add(X.get(0)); exemplarIndex.add(0); double t = Math.sqrt(2 * (1 - threshold)); //Add first element of test spectra to set of exemplar spectra exemplarLabel[0] = 0; boolean flag; double maxCos, sigmaMin, sigmaMax, dotXR, dotER, cosTheta; double[] vecR = new double[nDim]; for (int i = 0; i < nDim; i++) { vecR[i] = 1 / Math.sqrt(nDim); } ArrayRealVector R = new ArrayRealVector(vecR); ArrayRealVector exemplarSpec, testSpec; for (int i = 0; i < X.size(); i++) { if (i == 0 || exemplarLabel[i] == -1) { continue; } flag = false; maxCos = 0; testSpec = X.get(i); dotXR = testSpec.dotProduct(R); sigmaMin = dotXR - t; sigmaMax = dotXR + t; for (int j = 0; j < E.size(); j++) { exemplarSpec = E.get(j); dotER = exemplarSpec.dotProduct(R); if (dotER < sigmaMax && dotER > sigmaMin) { cosTheta = testSpec.dotProduct(exemplarSpec); if (cosTheta > threshold) { //Test spectra is similar to one of the exemplar spectra if (cosTheta > maxCos) { maxCos = cosTheta; exemplarLabel[i] = j; //System.out.println("Count: "+i+"\texemplarLabel: "+exemplarLabel[i]); flag = true; } } } } if (!flag) { //Test spectra is unique, add it to set of exemplars E.add(testSpec); exemplarIndex.add(i); exemplarLabel[i] = E.size() - 1; //System.out.println("Count: "+i+"\texemplarLabel: "+exemplarLabel[i]); } } return exemplarIndex; }
From source file:fp.overlapr.algorithmen.StressMajorization.java
@Deprecated private static ArrayRealVector conjugateGradientsMethod(Array2DRowRealMatrix A, ArrayRealVector b, ArrayRealVector werte) {/*from www . ja va 2s . c om*/ Array2DRowRealMatrix preJacobi = new Array2DRowRealMatrix(A.getRowDimension(), A.getColumnDimension()); // Predconditioner berechnen preJacobi.walkInRowOrder(new DefaultRealMatrixChangingVisitor() { @Override public double visit(int row, int column, double value) { if (row == column) { return 1 / A.getEntry(row, column); } else { return 0.0; } } }); // x_k beliebig whlen ArrayRealVector x_k = new ArrayRealVector(werte); // r_k berechnen ArrayRealVector r_k = b.subtract(A.operate(x_k)); // h_k berechnen ArrayRealVector h_k = (ArrayRealVector) preJacobi.operate(r_k); // d_k = r_k ArrayRealVector d_k = h_k; // x_k+1 und r_k+1 und d_k+1, sowie alpha und beta und z erzeugen ArrayRealVector x_k1; ArrayRealVector r_k1; ArrayRealVector d_k1; ArrayRealVector h_k1; double alpha; double beta; ArrayRealVector z; do { // Speichere Matrix-Vektor-Produkt, um es nur einmal auszurechnen z = (ArrayRealVector) A.operate(d_k); // Finde von x_k in Richtung d_k den Ort x_k1 des Minimums der // Funktion E // und aktualisere den Gradienten bzw. das Residuum alpha = r_k.dotProduct(h_k) / d_k.dotProduct(z); x_k1 = x_k.add(d_k.mapMultiply(alpha)); r_k1 = r_k.subtract(z.mapMultiply(alpha)); h_k1 = (ArrayRealVector) preJacobi.operate(r_k1); // Korrigiere die Suchrichtung d_k1 beta = r_k1.dotProduct(h_k1) / r_k.dotProduct(h_k); d_k1 = h_k1.add(d_k.mapMultiply(beta)); // Werte "eins" weitersetzen x_k = x_k1; r_k = r_k1; d_k = d_k1; h_k = h_k1; } while (r_k1.getNorm() > TOL); return x_k1; }
From source file:gamlss.algorithm.RSAlgorithm.java
/** * The Rigby and Stasinopoulos fitting algorithm. * Uses WLSMultipleLinearRegression to calculate the fitting * coefficents to define the fitted values of the distribution * parameters, also updating them and storing them in the * supplied distribution class (distr). * /*from ww w .ja va 2s . c o m*/ * @param distr - object of the fitted distribution belonging * to the gamlss family * @param response - vector of response variable values * @param w - vector of the weight values * the original data or not * @param designMatrices - design matrices for * each of the distribution parameters * @param smoothMatrices - smoother matrices for each * of the distribution parameters * distribution parameters */ public void functionRS(final GAMLSSFamilyDistribution distr, final ArrayRealVector response, final ArrayRealVector w) { //iter <- control$iter iter = Controls.ITERATIONS; //conv <- FALSE conv = false; //G.dev <- sum(w*G.dev.incr) gDev = w.dotProduct(distr.globalDevianceIncreament(response)); //G.dev.old <- G.dev+1 gDevOld = gDev + 1; //while ( abs(G.dev.old-G.dev) > c.crit && iter < n.cyc ){ while (FastMath.abs(gDevOld - gDev) > Controls.GAMLSS_CONV_CRIT && iter < Controls.GAMLSS_NUM_CYCLES) { //if ("mu"%in%names(family$parameters)){ //if ("sigma"%in%names(family$parameters)){ //if ("nu"%in%names(family$parameters)){ //if ("tau"%in%names(family$parameters)){ for (int i = 1; i < distr.getNumberOfDistribtionParameters() + 1; i++) { switch (i) { case DistributionSettings.MU: whichDistParameter = DistributionSettings.MU; break; case DistributionSettings.SIGMA: whichDistParameter = DistributionSettings.SIGMA; break; case DistributionSettings.NU: whichDistParameter = DistributionSettings.NU; break; case DistributionSettings.TAU: whichDistParameter = DistributionSettings.TAU; break; default: System.err.println(" Distribution parameter " + "is not allowed"); } //if (family$parameter$mu==TRUE & mu.fix==FALSE){ glimfit.setWLSnoIntercept(Controls.NO_INTERCEPT[whichDistParameter - 1]); // mu.fit <<- glim.fit(f = mu.object, X = mu.X, y = y, w = w, //fv = mu, os = mu.offset, step = mu.step,control = i.control, //gd.tol = gd.tol,auto = autostep) glimfit.glimFitFunctionRS(whichDistParameter); } //G.dev.old <- G.dev gDevOld = gDev; //G.dev <- sum(w*G.dev.incr) gDev = w.dotProduct(distr.globalDevianceIncreament(response)); //iter <- iter+1 iter = iter + 1; //fiter <<- iter // if(trace) cat("GAMLSS-RS iteration ", iter, ": Global Deviance //= ", format(round(G.dev, 4)), " \n", sep = "") if (Controls.GAMLSS_TRACE) { System.out.println("GAMLSS-RS iteration " + iter + " : Global Deviance = " + gDev); } //if (G.dev > (G.dev.old+gd.tol) && iter >1 ) stop(paste("The //global deviance is increasing", "\n", "Try different //steps for the parameters or the model maybe inappropriate")) if (gDev > (gDevOld + Controls.GLOB_DEVIANCE_TOL) && iter > 1) { System.err.println("The global deviance is increasing, " + "Try different steps for the parameters or " + "the model maybe inappropriate"); break; } //if ( abs(G.dev.old-G.dev) < c.crit ) conv <- TRUE else FALSE if (FastMath.abs(gDevOld - gDev) < Controls.GAMLSS_CONV_CRIT) { conv = true; } } //if (!conv && no.warn ) if (!conv) { //warning("Algorithm RS has not yet converged") System.err.println("Algorithm RS has not yet converged"); } }
From source file:gamlss.algorithm.CGAlgorithm.java
public Hashtable<Integer, ArrayRealVector> CGfunction(GAMLSSFamilyDistribution distr, ArrayRealVector response, Hashtable<Integer, BlockRealMatrix> X, ArrayRealVector weights) { //iter <- control$iter double iter = Controls.ITERATIONS; //conv <- FALSE boolean conv = false; //G.dev.incr <- eval(G.dev.expr) //G.dev <- sum(w*G.dev.incr) double gDev = weights.dotProduct(distr.globalDevianceIncreament(response)); //G.dev.old <- G.dev+1 double gDevOld = gDev + 1; //while ( abs(G.dev.old-G.dev) > c.crit && iter < n.cyc ) while (FastMath.abs(gDevOld - gDev) > Controls.GAMLSS_CONV_CRIT && iter < nCyc) { for (int i = 1; i < distr.getNumberOfDistribtionParameters() + 1; i++) { switch (i) { case DistributionSettings.MU: whichDistParameter = DistributionSettings.MU; break; case DistributionSettings.SIGMA: whichDistParameter = DistributionSettings.SIGMA; break; case DistributionSettings.NU: whichDistParameter = DistributionSettings.NU; break; case DistributionSettings.TAU: whichDistParameter = DistributionSettings.TAU; break; }//from w ww .j a v a 2 s. co m //eta.mu <- eta.old.mu <- family$mu.linkfun(mu) etaOld = makelink.link(distr.getDistributionParameterLink(whichDistParameter), distr.getDistributionParameter(whichDistParameter)); cgData.put("etaOld" + whichDistParameter, etaOld); eta = etaOld; cgData.put("eta" + whichDistParameter, eta); //u.mu <- mu.object$dldp(mu=mu) ArrayRealVector u = distr.firstDerivative(whichDistParameter, response); //u2.mu <- mu.object$d2ldp2(mu=mu) ArrayRealVector u2 = distr.secondDerivative(whichDistParameter, response); if (whichDistParameter == DistributionSettings.SIGMA) { u2MuSigma = distr.secondCrossDerivative(DistributionSettings.MU, DistributionSettings.SIGMA, response); } //dr.mu <- family$mu.dr(eta.mu) dr = makelink.distParameterEta(distr.getDistributionParameterLink(whichDistParameter), eta); //dr.mu <- 1/dr.mu dr = drInverse(dr); cgData.put("dr" + whichDistParameter, dr); //who.mu <- mu.object$who //smooth.frame.mu <- mu.object$smooth.frame //s.mu <- if(first.iter) mu.object$smooth else s.mu //w.mu <- -u2.mu/(dr.mu*dr.mu) w = wtSet(u2, dr); cgData.put("w" + whichDistParameter, w); if (whichDistParameter == DistributionSettings.SIGMA) { wMuSigma = wMUSIGMAset(u2MuSigma, cgData.get("dr" + DistributionSettings.MU), dr); cgData.put("wMuSigma", wMuSigma); } //z.mu <- (eta.old.mu-mu.offset)+mu.step*u.mu/(dr.mu*w.mu) z = wvSet(etaOld, Controls.STEP[whichDistParameter - 1], Controls.OFFSET[whichDistParameter - 1], u, dr, w); cgData.put("z" + whichDistParameter, z); } for (int i = 1; i < distr.getNumberOfDistribtionParameters() + 1; i++) { switch (i) { case DistributionSettings.MU: whichDistParameter = DistributionSettings.MU; break; case DistributionSettings.SIGMA: whichDistParameter = DistributionSettings.SIGMA; break; case DistributionSettings.NU: whichDistParameter = DistributionSettings.NU; break; case DistributionSettings.TAU: whichDistParameter = DistributionSettings.TAU; break; } //if (family$parameter$mu==TRUE & mu.fix==FALSE){ glimfitCG.setWLSnoIntercept(Controls.NO_INTERCEPT[whichDistParameter - 1]); bettas.put(whichDistParameter, glimfitCG.glimFitFunctionCG(distr, response, X, distr.getDistributionParameter(whichDistParameter), weights, whichDistParameter, Controls.STEP[whichDistParameter - 1], Controls.OFFSET[whichDistParameter - 1], gDev, cgData, makelink)); } gDevOld = gDev; //G.dev.incr <- eval(G.dev.expr) //G.dev <- sum(w*G.dev.incr) gDev = weights.dotProduct(distr.globalDevianceIncreament(response)); //if (G.dev < G.dev.old) break // if (gDev < gDevOld) // { // break; // } //iter <- iter+1 iter = iter + 1; //if(trace) if (Controls.GAMLSS_TRACE) { //cat("GAMLSS-CG iteration ", iter, ": Global Deviance = ", format(round(G.dev, 4)), " \n", sep = "") System.out.println("GAMLSS-CG iteration " + iter + " : Global Deviance = " + gDev); } //if (G.dev > (G.dev.old+gd.tol) && iter >1 ) if (gDev > (gDevOld + Controls.GLOB_DEVIANCE_TOL) && iter > 1) { //stop(paste("The global deviance is increasing in CG-algorithm ", "\n", "Try different steps for the parameters or the model maybe inappropriate")) System.err.println( "The global deviance is increasing in CG-algorithm, Try different steps for the parameters or the model maybe inappropriate"); break; } //if ( abs(G.dev.old-G.dev) < c.crit ) conv <- TRUE else FALSE if (FastMath.abs(gDevOld - gDev) < Controls.GAMLSS_CONV_CRIT) { conv = true; } else { conv = false; } // if (!conv) if (!conv && iter == nCyc) { //warning("Algorithm CG has not yet converged"); System.out.println("Algorithm CG has not yet converged"); } } return bettas; }
From source file:edu.utexas.cs.tactex.subscriptionspredictors.LWRCustOldAppache.java
/** * LWR prediction// w w w . j av a 2 s.c om * * @param X * @param Y * @param x0 * @param tau * @return */ public Double LWRPredict(ArrayRealVector X, ArrayRealVector Y, double x0, final double tau) { ArrayRealVector X0 = new ArrayRealVector(X.getDimension(), x0); ArrayRealVector delta = X.subtract(X0); ArrayRealVector sqDists = delta.ebeMultiply(delta); UnivariateFunction expTau = new UnivariateFunction() { @Override public double value(double arg0) { //log.info(" cp univariate tau " + tau); return Math.pow(Math.E, -arg0 / (2 * tau)); } }; ArrayRealVector W = sqDists.map(expTau); double Xt_W_X = X.dotProduct(W.ebeMultiply(X)); if (Xt_W_X == 0.0) { log.error(" cp LWR cannot predict - 0 denominator returning NULL"); log.error("Xcv is " + X.toString()); log.error("Ycv is " + Y.toString()); log.error("x0 is " + x0); return null; // <==== NOTE: a caller must be prepared for it } double theta = (1.0 / Xt_W_X) * X.ebeMultiply(W).dotProduct(Y); return theta * x0; }
From source file:gamlss.algorithm.GlimFitCG.java
public ArrayRealVector glimFitFunctionCG(GAMLSSFamilyDistribution distr, ArrayRealVector response, Hashtable<Integer, BlockRealMatrix> X, ArrayRealVector fv, ArrayRealVector weights, int whichDistParameter, double step, double offSet, double globalDeviance, Hashtable<String, ArrayRealVector> cgData, MakeLinkFunction makelink) { int itn = 0;//from w w w . j av a2 s . co m double[] temp = new double[0]; //G.dev.in <- G.dev+1 // double gDevIn = globalDeviance +1; //i.G.dev <- G.dev double gDev = globalDeviance; double gDevOld = globalDeviance + 1; //first.iter <- FALSE boolean firstIter = false; //while ( abs(G.dev.in -i.G.dev) > i.c.crit && i.iter < i.n.cyc ) while (FastMath.abs(gDevOld - gDev) > Controls.GLIM_CONV_CRIT && itn < Controls.GLIM_NUM_CYCLES) { //i.iter <- i.iter+1 itn = itn + 1; //for (int i = 1; i < distr.getNumberOfDistribtionParameters()+1; i++ ){ switch (whichDistParameter) { case DistributionSettings.MU: // whichDistParameter = DistributionSettings.MU; //adj.mu <- -(w.mu.sigma*(eta.sigma-eta.old.sigma)+w.mu.nu*(eta.nu-eta.old.nu)+w.mu.tau*(eta.tau-eta.old.tau))/w.mu adj = setAdj(cgData.get("wMuSigma"), cgData.get("eta" + DistributionSettings.SIGMA), cgData.get("etaOld" + DistributionSettings.SIGMA), cgData.get("wMuNu"), cgData.get("eta" + DistributionSettings.NU), cgData.get("etaOld" + DistributionSettings.NU), cgData.get("wMuTau"), cgData.get("eta" + DistributionSettings.TAU), cgData.get("etaOld" + DistributionSettings.TAU), cgData.get("w" + DistributionSettings.MU)); cgData.put("adj" + DistributionSettings.MU, adj); break; case DistributionSettings.SIGMA: // whichDistParameter = DistributionSettings.SIGMA; // adj.sigma <- -(w.mu.sigma*(eta.mu-eta.old.mu)+ w.sigma.nu*(eta.nu-eta.old.nu)+w.sigma.tau*(eta.tau-eta.old.tau))/w.sigma adj = setAdj(cgData.get("wMuSigma"), cgData.get("eta" + DistributionSettings.MU), cgData.get("etaOld" + DistributionSettings.MU), cgData.get("wSigmaNu"), cgData.get("eta" + DistributionSettings.NU), cgData.get("etaOld" + DistributionSettings.NU), cgData.get("wSigmaTau"), cgData.get("eta" + DistributionSettings.TAU), cgData.get("etaOld" + DistributionSettings.TAU), cgData.get("w" + DistributionSettings.SIGMA)); cgData.put("adj" + DistributionSettings.SIGMA, adj); break; case DistributionSettings.NU: // whichDistParameter = DistributionSettings.NU; break; case DistributionSettings.TAU: // whichDistParameter = DistributionSettings.TAU; break; } //wv.mu <- z.mu+adj.mu ArrayRealVector wv = setWV(cgData.get("z" + whichDistParameter), cgData.get("adj" + whichDistParameter)); //mu.fit <<- lm.wfit(x=mu.X,y=wv.mu,w=w.mu*w,method="qr") wls.newSampleData(wv, X.get(whichDistParameter).copy(), cgData.get("w" + whichDistParameter).ebeMultiply(weights)); ArrayRealVector fit = (ArrayRealVector) wls.calculateFittedValues(false); //System.out.println(wls.calculateBeta()); // bettas.put(whichDistParameter, (ArrayRealVector) wls.calculateBeta()); //mu.fit$eta <<- eta.mu <- mu.fit$fitted.values+mu.offset temp = new double[fit.getDimension()]; for (int j = 0; j < temp.length; j++) { temp[j] = fit.getEntry(j) + offSet; } //eta = new ArrayRealVector(temp,false); cgData.put("eta" + whichDistParameter, new ArrayRealVector(temp, false)); temp = null; //mu.fit$fv <<- mu <<- mu.object$linkinv(eta.mu) ArrayRealVector dParam = makelink.linkInv(distr.getDistributionParameterLink(whichDistParameter), cgData.get("eta" + whichDistParameter)); distr.setDistributionParameter(whichDistParameter, dParam); //mu.fit$wv <<- wv.mu //mu.fit$wt <<- w.mu //G.dev.in <- i.G.dev gDevOld = gDev; //G.dev.incr <- eval(G.dev.expr) //i.G.dev <- sum(w*G.dev.incr) gDev = weights.dotProduct(distr.globalDevianceIncreament(response)); if (gDev > gDevOld && itn >= 2 && Controls.AUTO_STEP) { for (int i = 0; i < 5; i++) { //eta.mu <- (eta.mu+eta.old.mu)/2 ArrayRealVector etaM = etaMean(cgData.get("eta" + whichDistParameter), cgData.get("etaOld" + whichDistParameter)); cgData.put("eta" + whichDistParameter, etaM); //mu <<- mu.object$linkinv(eta.mu) distr.setDistributionParameter(whichDistParameter, makelink.linkInv(distr.getDistributionParameterLink(whichDistParameter), etaM)); //if(length(who.mu) > 0) s.mu <- (s.mu+s.mu.old)/2 // } } } //if(i.trace) if (Controls.GLIM_TRACE) { //cat("CG inner iteration ", iter, ": Global Deviance = ",format(round(i.G.dev, 4)), " \n", sep = "") System.err.println("CG inner iteration " + itn + " : Global Deviance = " + gDev); } //if (i.G.dev > (G.dev.in+gd.tol) && iter >1 ) if (gDev > (gDevOld + Controls.GLOB_DEVIANCE_TOL) && itn > 1) { //stop(paste("The global deviance is increasing in the inner CG loop", "\n","Try different steps for the parameters or the model maybe inappropriate")) System.out.println( "The global deviance is increasing in the inner CG loop, try different steps for the parameters or the model maybe inappropriate"); break; } } //G.dev.old <- G.dev // gDevOld = gDev; //G.dev.incr <- eval(G.dev.expr) //G.dev <- sum(w*G.dev.incr) // gDev = weights.dotProduct(distr.globalDevianceIncreament(response)); //if (G.dev > G.dev.old && iter >= 2 && autostep == TRUE) return (ArrayRealVector) wls.calculateBeta(); }