List of usage examples for org.apache.commons.math3.linear ArrayRealVector getDimension
@Override public int getDimension()
From source file:edu.stanford.cfuller.imageanalysistools.filter.LocalMaximumSeparabilityThresholdingFilter.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>// www. j a v a 2 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()); ArrayRealVector eta_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); eta_v.setEntry(c, eta); if (eta >= best_eta) { best_eta = eta; best_index = k; } } else { eta_v.setEntry(c, 0); } c++; } int orig_method_best_index = best_index; c = 1; ArrayList<Integer> maxima = new ArrayList<Integer>(); Map<Integer, Integer> k_by_c = new HashMap<Integer, Integer>(); Map<Integer, Integer> c_by_k = new HashMap<Integer, Integer>(); for (int k = h.getMinValue() + 1; k < h.getMaxValue(); k += increment) { //detect if this is a local maximum k_by_c.put(c, k); c_by_k.put(k, c); int lastEntryNotEqual = c - 1; int nextEntryNotEqual = c + 1; while (lastEntryNotEqual > 0 && eta_v.getEntry(lastEntryNotEqual) == eta_v.getEntry(c)) { --lastEntryNotEqual; } while (nextEntryNotEqual < (eta_v.getDimension() - 1) && eta_v.getEntry(nextEntryNotEqual) == eta_v.getEntry(c)) { ++nextEntryNotEqual; } if (eta_v.getEntry(c) > eta_v.getEntry(lastEntryNotEqual) && eta_v.getEntry(c) > eta_v.getEntry(nextEntryNotEqual)) { maxima.add(k); } c++; } //now that we have maxima, try doing a gaussian fit to find the positions. If there's only one, we need to guess at a second RealVector parameters = new ArrayRealVector(6, 0.0); int position0 = 0; int position1 = h.getMaxValue(); if (maxima.size() > 1) { double best_max = 0; double second_best_max = 0; int best_pos = 0; int second_best_pos = 0; for (Integer k : maxima) { int ck = c_by_k.get(k); double eta_k = eta_v.getEntry(ck); if (eta_k > best_max) { second_best_max = best_max; second_best_pos = best_pos; best_max = eta_k; best_pos = ck; } else if (eta_k > second_best_max) { second_best_max = eta_k; second_best_pos = ck; } } position0 = best_pos; position1 = second_best_pos; } else { position0 = c_by_k.get(maxima.get(0)); position1 = (eta_v.getDimension() - position0) / 2 + position0; } //make sure that position 1 is larger than position 0 if (position1 < position0) { int temp = position0; position0 = position1; position1 = temp; } double s = (position1 - position0) / 4.0; parameters.setEntry(0, eta_v.getEntry(position0));//*Math.sqrt(2*Math.PI)*s); parameters.setEntry(1, position0); parameters.setEntry(2, s); parameters.setEntry(3, eta_v.getEntry(position1));//*Math.sqrt(2*Math.PI)*s); parameters.setEntry(4, position1); parameters.setEntry(5, s); DoubleGaussianObjectiveFunction dgof = new DoubleGaussianObjectiveFunction(); dgof.setEta(eta_v); NelderMeadMinimizer nmm = new NelderMeadMinimizer(); RealVector result = nmm.optimize(dgof, parameters); best_index = (int) result.getEntry(4); if (k_by_c.containsKey(best_index)) { best_index = k_by_c.get(best_index); } else { //fall back to the normal global maximum if the fitting seems to have found an invalid value. best_index = orig_method_best_index; } thresholdValue = best_index; if (thresholdValue == Integer.MAX_VALUE) { thresholdValue = 0; } for (ImageCoordinate coord : im) { if (im.getValue(coord) < thresholdValue) im.setValue(coord, 0); } }
From source file:gamlss.algorithm.Gamlss.java
/** This is to emulate the Gamlss algorithm where * the user specifies response variable vector * and design matrix.// w ww.ja v a 2 s .c o m * * @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: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// ww w . j av a 2 s. c o 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:gamlss.algorithm.AdditiveFit.java
/** * Performs a smoothing process - creates an approximating function * that attempts to capture important patterns in the data, while * leaving out noise or other fine-scale structures/rapid phenomena. * @param xMat - design matrix//from ww w. j a va 2 s. c om * @param y - response variable * @param sMat - smoother matrix * @param distParameter - distribution parameter * @return matrix of smoothed values */ public BlockRealMatrix fitSmoother(final ArrayRealVector y, final ArrayRealVector weights, final BlockRealMatrix xMat, final BlockRealMatrix sMat, final int distParameter) { this.s = sMat.copy(); this.w = weights; this.whichDistParameter = distParameter; //residuals <- as.vector(y - s %*% rep(1, ncol(s))) tempV = new ArrayRealVector(s.getColumnDimension()); tempV.set(1.0); tempV = (ArrayRealVector) s.operate(tempV); residuals = MatrixFunctions.vecSub(y, tempV); tempV = null; //fit <- list(fitted.values = 0) fittedValues = new ArrayRealVector(residuals.getDimension()); //rss <- weighted.mean(residuals^2, w) //rss = calculateRss(residuals, w); //tempArr = null; //rssold <- rss * 10 //rssOld = rss*10; //nit <- 0 nit = 0; //df <- rep(NA, length(who)) //lambda <- rep(NA, length(who)) //double[] df = new double[s.getColumnDimension()]; //var <- s var = s.copy(); //if(trace) cat("\nADDITIVE iter rss/n term\n") if (Controls.BACKFIT_TRACE) { System.out.println("ADDITIVE iter rss/n term"); } //ndig <- -log10(tol) + 1 //double ndig = -FastMath.log10(tol)+1; //RATIO <- tol + 1 ratio = Controls.BACKFIT_TOL + 1; //while(RATIO > tol & nit < maxit) while (ratio > Controls.BACKFIT_TOL & nit < Controls.BACKFIT_CYCLES) { //rssold <- rss //rssOld = rss; //nit <- nit + 1 nit++; //z <- residuals + fit$fitted.values z = residuals.add(fittedValues); //org.apache.commons.lang3.time.StopWatch watch = new org.apache.commons.lang3.time.StopWatch(); //watch.start(); //fit <- lm.wfit(x, z, w, method="qr", ...) wls.setNoIntercept(Controls.NO_INTERCEPT[whichDistParameter - 1]); wls.newSampleData(z, xMat.copy(), w.copy()); wls.setNoIntercept(false); //residuals <- fit$residuals fittedValues = (ArrayRealVector) wls.calculateFittedValues(Controls.IS_SVD); //watch.stop(); //System.out.println(watch.getNanoTime()/(1e+9)); //watch.reset(); //residuals = z.subtract(fittedValues); //residuals = (ArrayRealVector) wls.calculateResiduals(); residuals = z.subtract(fittedValues); //rss <- weighted.mean(residuals^2, w) //rss = calculateRss(residuals, w); //if(trace) cat(nit, " ", //format(round(rss, ndig)), " Parametric -- lm.wfit\n", sep = "") if (Controls.BACKFIT_TRACE) { //System.out.println(" " + nit + " " + //rss + " Parametric -- lm.wfit"); } //deltaf <- 0 deltaf = 0; //for(j in seq(names.calls)) for (colNum = 0; colNum < s.getColumnDimension(); colNum++) { //old <- s[, j] old = (ArrayRealVector) s.getColumnVector(colNum); //z <- residuals + s[, j] z = residuals.add(old); //fit.call <- eval(smooth.calls[[j]]) //residuals <- as.double(fit.call$residuals) residuals = smoother.solve(this); //if(length(residuals) != n) //stop(paste(names.calls[j], "returns a vector //of the wrong length")) if (residuals.getDimension() != y.getDimension()) { System.err.println(colNum + " column of matrix s has a" + " vector of the wrong length"); } //s[, j] <- z - residual s.setColumnVector(colNum, z.subtract(residuals)); //if (length(fit.call$lambda)>1) //{for cases where there are multiple lambdas //ambda[j] <- fit.call$lambda[1] //coefSmo[[j]] <- if(is.null(fit.call$coefSmo)) //0 else fit.call$coefSmo //deltaf <- deltaf + weighted.mean((s[, j] - old)^2, w) tempV = MatrixFunctions.vecSub((ArrayRealVector) s.getColumnVector(colNum), old); deltaf = deltaf + mean.evaluate(tempV.ebeMultiply(tempV).getDataRef(), w.getDataRef()); tempV = null; //rss <- weighted.mean(residuals^2, w) //rss = calculateRss(residuals, w); // if(trace) if (Controls.BACKFIT_TRACE) { //cat(" ", nit, " ", format(round(rss, ndig)), //" Nonparametric -- ", names.calls[j], "\n", sep = "") //System.out.println(" " + nit +" " + rss + " //Nonparametric " + "pb(column "+ i+")"); } //df[j] <- fit.call$nl.df //df[i] = pb.getEdf(); //if(se) if (Controls.IS_SE) { //var[, j] <- fit.call$var var.setColumnVector(colNum, smoother.getVar()); } } //RATIO <- sqrt(deltaf/sum(w * apply(s, 1, sum)^2)) tempD = 0.0; tempM = new BlockRealMatrix(s.getRowDimension(), 1); for (int j = 0; j < s.getRowDimension(); j++) { for (int k = 0; k < s.getColumnDimension(); k++) { tempD = tempD + s.getEntry(j, k); } tempM.setEntry(j, 0, tempD); tempD = 0.0; } size = tempM.getRowDimension(); for (int j = 0; j < size; j++) { tempD = tempD + tempM.getEntry(j, 0) * tempM.getEntry(j, 0) * w.getEntry(j); } ratio = FastMath.sqrt(deltaf / tempD); tempD = null; tempM = null; //if(trace) //cat("Relative change in functions:", //format(round(RATIO, ndig)), "\n") if (Controls.BACKFIT_TRACE) { System.out.println("Relative change in functions: " + ratio); } } //if(nit == maxit && maxit > 1) //warning(paste("additive.fit convergence not //obtained in ", maxit, " iterations")) if (ratio > Controls.BACKFIT_TOL) { System.out.println( "AdditiveFit convergence is not obtained in " + Controls.BACKFIT_CYCLES + " iterations"); } //fit$fitted.values <- y - residuals fittedValues = y.subtract(residuals); //rl <- c(fit, list(smooth = s, nl.df = //sum(df), lambda=lambda, coefSmo=coefSmo)) return s; }
From source file:syncleus.gremlann.model.Autoencoder.java
public void get_reconstructed_input(List<Vertex> y, List<Vertex> z) { ArrayRealVector vy = signals(y); for (int i = 0; i < getInputCount(); i++) { double zi = 0; for (int j = 0; j < vy.getDimension(); j++) { zi += weight(i, j) * vy.getEntry(j); }/* w w w . j av a 2 s . com*/ zi += inputBias[i]; Neuron.signal(z.get(i), sigmoid(zi)); } }