List of usage examples for org.apache.commons.math.analysis MultivariateRealFunction MultivariateRealFunction
MultivariateRealFunction
From source file:com.opengamma.analytics.math.util.wrapper.CommonsMathWrapper.java
/** * @param f An OG 1-D function mapping vectors of doubles onto doubles, not null * @return A Commons multivariate real function *//* w ww .j a v a 2 s . c o m*/ public static MultivariateRealFunction wrapMultivariate(final Function1D<DoubleMatrix1D, Double> f) { Validate.notNull(f); return new MultivariateRealFunction() { @Override public double value(final double[] point) throws FunctionEvaluationException, IllegalArgumentException { return f.evaluate(new DoubleMatrix1D(point)); } }; }
From source file:com.opengamma.analytics.math.util.wrapper.CommonsMathWrapper.java
/** * @param f An OG n-D function mapping doubles onto doubles, not null * @return A Commons multivariate real function *//*from ww w . j a v a2 s . c om*/ public static MultivariateRealFunction wrap(final FunctionND<Double, Double> f) { Validate.notNull(f); return new MultivariateRealFunction() { @Override public double value(final double[] point) throws FunctionEvaluationException, IllegalArgumentException { final int n = point.length; final Double[] coordinate = new Double[n]; for (int i = 0; i < n; i++) { coordinate[i] = point[i]; } return f.evaluate(coordinate); } }; }
From source file:endrov.nucAutoJH.FitGaussian.java
private static double[] fitGaussian2D_(EvPixels p, double sigmaInit, final double midxInit, final double midyInit) { //sigma00, sigma01, sigma11, mu_x, mu_y, c p = p.getReadOnly(EvPixelsType.DOUBLE); final double[] arrPixels = p.getArrayDouble(); final int w = p.getWidth(); final int h = p.getHeight(); int extent = (int) Math.round(3 * sigmaInit); extent = Math.max(extent, 2); final int sx = Math.max(0, (int) (midxInit - extent)); final int ex = Math.min(w, (int) (midxInit + extent + 1)); //+1 to the right? final int sy = Math.max(0, (int) (midyInit - extent)); final int ey = Math.min(h, (int) (midyInit + extent + 1)); double minIntensity = Double.MAX_VALUE; double maxIntensity = Double.MIN_VALUE; for (int y = sy; y < ey; y++) { int base = y * w; double dy2 = y - midyInit; dy2 = dy2 * dy2;//from ww w .j av a 2 s .c o m for (int x = sx; x < ex; x++) { double dx2 = x - midxInit; dx2 = dx2 * dx2; double t = arrPixels[base + x]; //if(dx2+dy2<=extent*extent) { if (t < minIntensity) minIntensity = t; if (t > maxIntensity) maxIntensity = t; } } } //double[] weights=new double[]{1}; double[] startPoint = new double[] { sigmaInit, 0, sigmaInit, midxInit, midyInit, minIntensity, maxIntensity - minIntensity }; //double[] output=new double[startPoint.length]; try { MultivariateRealFunction func = new MultivariateRealFunction() { // opt.optimize( public double value(double[] arg) throws FunctionEvaluationException, IllegalArgumentException { double sigma00 = arg[0]; double sigma01 = arg[1]; double sigma11 = arg[2]; double mu0 = arg[3]; double mu1 = arg[4]; double C = arg[5]; double D = arg[6]; double sumError = 0; Matrix2d sigma = new Matrix2d(sigma00, sigma01, sigma01, sigma11); Matrix2d sigmaInv = new Matrix2d(); sigma.invert(sigmaInv); double sigmaDet = sigma.determinant(); double front = 1.0 / (2 * Math.PI * Math.sqrt(sigmaDet)); //System.out.println("front: "+front); //System.out.println("sigma inv "+sigmaInv); if (mu0 < sx || mu0 > ex) sumError += 1000000; if (mu1 < sy || mu1 > ey) sumError += 1000000; if (sigma00 < 1) sumError += 1000000; //if(sigma01<0) sumError+=1000000; if (sigma11 < 1) sumError += 1000000; if (D <= 0) sumError += 1000000; for (int y = sy; y < ey; y++) { int base = y * w; double dy2 = y - midyInit; dy2 = dy2 * dy2; for (int x = sx; x < ex; x++) { double dx2 = x - midxInit; dx2 = dx2 * dx2; double thisReal = arrPixels[base + x]; // if(dx2+dy2<=extent*extent) { // DoubleMatrix2D sigma=new DenseDoubleMatrix2D(new double[][]{{sigma00,sigma01},{sigma01,sigma11}}); //double sigmaDet=sigma00*sigma11-sigma01*sigma01; double dx0 = x - mu0; double dx1 = y - mu1; //http://en.wikipedia.org/wiki/Multivariate_normal_distribution Vector2d vX = new Vector2d(dx0, dx1); Vector2d op = new Vector2d(vX); sigmaInv.transform(op); double upper = -0.5 * op.dot(vX); double exp = Math.exp(upper); //System.out.println("front "+front+" "+exp+" C "+C+" thisreal"+thisReal+" upper "+upper); if (upper > -0.4) exp = 1; else exp = Math.max(0, 1 + upper + 0.4); /* if(exp<0.7) exp=0; else exp=1; */ double thisExpected = D * front * exp + C; double diff = thisExpected - thisReal; sumError += diff * diff; } } } //System.out.println(sigma00+"\t"+sigma01+"\t"+sigma11+"\tC"+C+"\tmu "+mu0+","+mu1+"\terr "+sumError); return sumError; // return new double[]{sumError}; } }; NelderMead opt = new NelderMead(); //LevenbergMarquardtOptimizer opt=new LevenbergMarquardtOptimizer(); opt.setMaxIterations(10000); RealPointValuePair pair = opt.optimize(func, GoalType.MINIMIZE, startPoint); int numit = opt.getIterations(); System.out.println("#it " + numit); System.out.println("err " + func.value(pair.getPointRef())); return pair.getPointRef(); // for(int i=0;i<startPoint.length;i++) // System.out.println("i: "+i+" "+output[i]); //, output, weights, startPoint); } /* catch (MaxIterationsExceededException e) { System.out.println("max it reached"); }*/ catch (Exception e) { e.printStackTrace(); } //Maybe this is a bad point? System.out.println("max it reached"); return startPoint; // return output; }
From source file:endrov.typeLineageAutoNucJH.FitGaussian.java
private static double[] fitGaussian2D_(EvPixels p, double sigmaInit, final double midxInit, final double midyInit) { //sigma00, sigma01, sigma11, mu_x, mu_y, c p = p.getReadOnly(EvPixelsType.DOUBLE); final double[] arrPixels = p.getArrayDouble(); final int w = p.getWidth(); final int h = p.getHeight(); int extent = (int) Math.round(3 * sigmaInit); extent = Math.max(extent, 2); final int sx = Math.max(0, (int) (midxInit - extent)); final int ex = Math.min(w, (int) (midxInit + extent + 1)); //+1 to the right? final int sy = Math.max(0, (int) (midyInit - extent)); final int ey = Math.min(h, (int) (midyInit + extent + 1)); double minIntensity = Double.MAX_VALUE; double maxIntensity = -Double.MAX_VALUE; for (int y = sy; y < ey; y++) { int base = y * w; double dy2 = y - midyInit; dy2 = dy2 * dy2;/*from w w w .j a v a2 s . co m*/ for (int x = sx; x < ex; x++) { double dx2 = x - midxInit; dx2 = dx2 * dx2; double t = arrPixels[base + x]; //if(dx2+dy2<=extent*extent) { if (t < minIntensity) minIntensity = t; if (t > maxIntensity) maxIntensity = t; } } } //double[] weights=new double[]{1}; double[] startPoint = new double[] { sigmaInit, 0, sigmaInit, midxInit, midyInit, minIntensity, maxIntensity - minIntensity }; //double[] output=new double[startPoint.length]; try { MultivariateRealFunction func = new MultivariateRealFunction() { // opt.optimize( public double value(double[] arg) throws FunctionEvaluationException, IllegalArgumentException { double sigma00 = arg[0]; double sigma01 = arg[1]; double sigma11 = arg[2]; double mu0 = arg[3]; double mu1 = arg[4]; double C = arg[5]; double D = arg[6]; double sumError = 0; Matrix2d sigma = new Matrix2d(sigma00, sigma01, sigma01, sigma11); Matrix2d sigmaInv = new Matrix2d(); sigma.invert(sigmaInv); double sigmaDet = sigma.determinant(); double front = 1.0 / (2 * Math.PI * Math.sqrt(sigmaDet)); //System.out.println("front: "+front); //System.out.println("sigma inv "+sigmaInv); if (mu0 < sx || mu0 > ex) sumError += 1000000; if (mu1 < sy || mu1 > ey) sumError += 1000000; if (sigma00 < 1) sumError += 1000000; //if(sigma01<0) sumError+=1000000; if (sigma11 < 1) sumError += 1000000; if (D <= 0) sumError += 1000000; for (int y = sy; y < ey; y++) { int base = y * w; double dy2 = y - midyInit; dy2 = dy2 * dy2; for (int x = sx; x < ex; x++) { double dx2 = x - midxInit; dx2 = dx2 * dx2; double thisReal = arrPixels[base + x]; // if(dx2+dy2<=extent*extent) { // DoubleMatrix2D sigma=new DenseDoubleMatrix2D(new double[][]{{sigma00,sigma01},{sigma01,sigma11}}); //double sigmaDet=sigma00*sigma11-sigma01*sigma01; double dx0 = x - mu0; double dx1 = y - mu1; //http://en.wikipedia.org/wiki/Multivariate_normal_distribution Vector2d vX = new Vector2d(dx0, dx1); Vector2d op = new Vector2d(vX); sigmaInv.transform(op); double upper = -0.5 * op.dot(vX); double exp = Math.exp(upper); //System.out.println("front "+front+" "+exp+" C "+C+" thisreal"+thisReal+" upper "+upper); if (upper > -0.4) exp = 1; else exp = Math.max(0, 1 + upper + 0.4); /* if(exp<0.7) exp=0; else exp=1; */ double thisExpected = D * front * exp + C; double diff = thisExpected - thisReal; sumError += diff * diff; } } } //System.out.println(sigma00+"\t"+sigma01+"\t"+sigma11+"\tC"+C+"\tmu "+mu0+","+mu1+"\terr "+sumError); return sumError; // return new double[]{sumError}; } }; NelderMead opt = new NelderMead(); //LevenbergMarquardtOptimizer opt=new LevenbergMarquardtOptimizer(); opt.setMaxIterations(10000); RealPointValuePair pair = opt.optimize(func, GoalType.MINIMIZE, startPoint); int numit = opt.getIterations(); System.out.println("#it " + numit); System.out.println("err " + func.value(pair.getPointRef())); return pair.getPointRef(); // for(int i=0;i<startPoint.length;i++) // System.out.println("i: "+i+" "+output[i]); //, output, weights, startPoint); } /* catch (MaxIterationsExceededException e) { System.out.println("max it reached"); }*/ catch (Exception e) { e.printStackTrace(); } //Maybe this is a bad point? System.out.println("max it reached"); return startPoint; // return output; }
From source file:uk.ac.diamond.scisoft.analysis.optimize.ApacheConjugateGradient.java
@Override public void optimize(IDataset[] coords, IDataset data, final IFunction function) throws Exception { // pull out appropriate data from the inputs in a way that they can be used final int numCoords = coords.length; final DoubleDataset[] newCoords = new DoubleDataset[numCoords]; for (int i = 0; i < numCoords; i++) { newCoords[i] = (DoubleDataset) DatasetUtils.convertToAbstractDataset(coords[i]) .cast(AbstractDataset.FLOAT64); }/*from ww w . j a v a2 s . co m*/ final DoubleDataset values = (DoubleDataset) DatasetUtils.convertToAbstractDataset(data) .cast(AbstractDataset.FLOAT64); NonLinearConjugateGradientOptimizer cg = new NonLinearConjugateGradientOptimizer( ConjugateGradientFormula.POLAK_RIBIERE); //FIXME - this does not work very well for many cases. DifferentiableMultivariateRealFunction f = new DifferentiableMultivariateRealFunction() { @Override public double value(double[] arg0) throws FunctionEvaluationException, IllegalArgumentException { function.setParameterValues(arg0); return function.residual(true, values, newCoords); } @Override public MultivariateRealFunction partialDerivative(final int parameter) { MultivariateRealFunction result = new MultivariateRealFunction() { @Override public double value(double[] parameters) throws FunctionEvaluationException, IllegalArgumentException { double step = 0.1; parameters[parameter] -= step; function.setParameterValues(parameters); double min = function.residual(true, values, newCoords); parameters[parameter] += 2.0 * step; function.setParameterValues(parameters); double max = function.residual(true, values, newCoords); return -((max - min) / (2.0 * step)); } }; return result; } @Override public MultivariateVectorialFunction gradient() { MultivariateVectorialFunction result = new MultivariateVectorialFunction() { @Override public double[] value(double[] parameters) throws FunctionEvaluationException, IllegalArgumentException { double[] result = new double[parameters.length]; for (int i = 0; i < parameters.length; i++) { result[i] = partialDerivative(i).value(parameters); } return result; } }; return result; } }; double[] start = function.getParameterValues(); RealPointValuePair result = cg.optimize(f, GoalType.MINIMIZE, start); function.setParameterValues(result.getPoint()); }
From source file:uk.ac.diamond.scisoft.analysis.optimize.ApacheMultiDirectional.java
@Override public void optimize(IDataset[] coords, IDataset data, final IFunction function) throws Exception { // Pull out the data which is required from the inputs final int numCoords = coords.length; final DoubleDataset[] newCoords = new DoubleDataset[numCoords]; for (int i = 0; i < numCoords; i++) { newCoords[i] = (DoubleDataset) DatasetUtils.convertToAbstractDataset(coords[i]) .cast(AbstractDataset.FLOAT64); }// ww w. j a v a2 s . c o m final DoubleDataset values = (DoubleDataset) DatasetUtils.convertToAbstractDataset(data) .cast(AbstractDataset.FLOAT64); // create an instance of the fitter MultiDirectional md = new MultiDirectional(); // provide the fitting function which wrappers all the normal fitting functionality MultivariateRealFunction f1 = new MultivariateRealFunction() { @Override public double value(double[] arg0) throws FunctionEvaluationException, IllegalArgumentException { function.setParameterValues(arg0); return function.residual(true, values, newCoords); } }; double[] start = function.getParameterValues(); // preform the optimisation RealPointValuePair result = md.optimize(f1, GoalType.MINIMIZE, start); // set the input functions parameters to be the result before finishing. function.setParameterValues(result.getPoint()); }
From source file:uk.ac.diamond.scisoft.analysis.optimize.ApacheNelderMead.java
@Override public void optimize(IDataset[] coords, IDataset data, final IFunction function) throws Exception { // Pull out the data which is required from the inputs final int numCoords = coords.length; final DoubleDataset[] newCoords = new DoubleDataset[numCoords]; for (int i = 0; i < numCoords; i++) { newCoords[i] = (DoubleDataset) DatasetUtils.convertToAbstractDataset(coords[i]) .cast(AbstractDataset.FLOAT64); }//from w ww . j a va 2 s .c om final DoubleDataset values = (DoubleDataset) DatasetUtils.convertToAbstractDataset(data) .cast(AbstractDataset.FLOAT64); // create an instance of the fitter NelderMead nm = new NelderMead(); // provide the fitting function which wrappers all the normal fitting functionality MultivariateRealFunction f1 = new MultivariateRealFunction() { @Override public double value(double[] arg0) throws FunctionEvaluationException, IllegalArgumentException { function.setParameterValues(arg0); return function.residual(true, values, newCoords); } }; // preform the optimisation double[] start = function.getParameterValues(); RealPointValuePair result = nm.optimize(f1, GoalType.MINIMIZE, start); // set the input functions parameters to be the result before finishing. function.setParameterValues(result.getPoint()); }