Example usage for org.apache.commons.math3.optim.nonlinear.scalar.noderiv CMAESOptimizer CMAESOptimizer

List of usage examples for org.apache.commons.math3.optim.nonlinear.scalar.noderiv CMAESOptimizer CMAESOptimizer

Introduction

In this page you can find the example usage for org.apache.commons.math3.optim.nonlinear.scalar.noderiv CMAESOptimizer CMAESOptimizer.

Prototype

public CMAESOptimizer(int maxIterations, double stopFitness, boolean isActiveCMA, int diagonalOnly,
        int checkFeasableCount, RandomGenerator random, boolean generateStatistics,
        ConvergenceChecker<PointValuePair> checker) 

Source Link

Usage

From source file:norbert.mynemo.core.selection.UserRecommenderSelector.java

/**
 * Evaluates the given recommender type for several configurations. The number of neighbors for
 * each evaluation is chosen via an optimizer. The number of evaluation performed depends of the
 * number of steps necessary for the evaluation results to converge. The number of evaluations
 * performed is also bounded by an internal maximum.
 *
 * <p>/*from   w ww  .ja v a 2  s.  c o m*/
 * The given recommender type must be part of the user similarity based family.
 *
 * @return all evaluations done during the selection process
 */
public Collection<RecommenderEvaluation> select(RecommenderType type, double minimumCoverage) {
    checkArgument(type.getFamily() == RecommenderFamily.USER_SIMILARITY_BASED);
    checkArgument(0 <= minimumCoverage && minimumCoverage <= 1,
            "The minimum coverage must be" + "between 0 and 1.");

    final UserBasedRecommenderEvaluationFunction function = new UserBasedRecommenderEvaluationFunction(
            configuration, type, minimumCoverage);

    // find a good initial guess: test 10%, 20% 90% of maximum neighbors
    double secondBestInitialValue = 1;
    double bestInitialValue = 0;
    double bestEvalResult = Integer.MAX_VALUE;
    for (double factor = 0.1; factor <= 0.9; factor += 0.1) {
        double currentNumNeighbors = maxNeighbors * factor;
        double evalResult = function.value(currentNumNeighbors);
        if (evalResult < bestEvalResult) {
            bestEvalResult = evalResult;
            secondBestInitialValue = bestInitialValue;
            bestInitialValue = currentNumNeighbors;
        }
    }

    // initialize the quality parameters of the CMAES optimizer
    final double minNeighbors = 1;
    final double initialGuessNeighbors = bestInitialValue;
    final double sigmaNeighbors = Math.min(Math.abs(bestInitialValue - secondBestInitialValue),
            Math.min(maxNeighbors - bestInitialValue, bestInitialValue - 1));
    // not sure about that
    final int populationSize = (int) (Math.log10(maxNeighbors) / Math.log10(2));
    // inhibit the exception throw if the maximum number of evaluation is reached
    final int maxCmaesEvaluations = Integer.MAX_VALUE;
    // not sure about that
    final int maxCmaesIterations = (int) (Math.log10(maxNeighbors) / Math.log10(2));

    // initialize the other parameters
    final ConvergenceChecker<PointValuePair> checker = new MaxIterationChecker<PointValuePair>(20);
    final ObjectiveFunction objectiveFunction = new ObjectiveFunction(function);
    final CMAESOptimizer optimizer = new CMAESOptimizer(maxCmaesIterations, 1.0, true, 0, 0,
            new JDKRandomGenerator(), false, checker);

    final MaxEval maxEval = new MaxEval(maxCmaesEvaluations);
    final GoalType goalType = GoalType.MINIMIZE;
    final SimpleBounds bounds = new SimpleBounds(new double[] { minNeighbors, minNeighbors },
            new double[] { maxNeighbors, maxNeighbors });
    final InitialGuess initialGuess = new InitialGuess(
            new double[] { initialGuessNeighbors, initialGuessNeighbors });
    final CMAESOptimizer.PopulationSize popSize = new CMAESOptimizer.PopulationSize((populationSize));
    final CMAESOptimizer.Sigma sigma = new CMAESOptimizer.Sigma(
            new double[] { sigmaNeighbors, sigmaNeighbors });

    // run the optimizer
    optimizer.optimize(objectiveFunction, goalType, initialGuess, popSize, sigma, bounds, maxEval);

    return function.getEvaluations();
}

From source file:org.compevol.ssgd.MaximumLikelihood.java

public MaximumLikelihood(final MultivariateFunction likelihood, final Parameter variables) {
    this.likelihood = likelihood;
    this.variables = variables;
    initial = new double[variables.getDimension()];
    Arrays.fill(initial, 1);/*from   w  w  w . j  a  va  2 s  .  c om*/
    optimizer = new CMAESOptimizer(Integer.MAX_VALUE, 0.0, true, 0, 8096, new RandomGenerator() {

        @Override
        public void setSeed(int i) {
            throw new UnsupportedOperationException();
        }

        @Override
        public void setSeed(int[] ints) {
            throw new UnsupportedOperationException();
        }

        @Override
        public void setSeed(long l) {
            throw new UnsupportedOperationException();
        }

        @Override
        public void nextBytes(byte[] bytes) {
            MathUtils.nextBytes(bytes);
        }

        @Override
        public int nextInt() {
            return MathUtils.nextInt();
        }

        @Override
        public int nextInt(int i) {
            return MathUtils.nextInt(i);
        }

        @Override
        public long nextLong() {
            return MathUtils.nextLong();
        }

        @Override
        public boolean nextBoolean() {
            return MathUtils.nextBoolean();
        }

        @Override
        public float nextFloat() {
            return MathUtils.nextFloat();
        }

        @Override
        public double nextDouble() {
            return MathUtils.nextDouble();
        }

        @Override
        public double nextGaussian() {
            return MathUtils.nextGaussian();
        }
    }, true, new SimplePointChecker<PointValuePair>(MachineAccuracy.SQRT_EPSILON, MachineAccuracy.EPSILON));
}

From source file:org.dawnsci.plotting.tools.diffraction.BeamCenterRefinement.java

/**
 * Run optimisation of sector region position in a separate job. 
 * /*www .  j av  a  2  s  .c  o  m*/
 * @param startPosition Initial position of sector region
 */
public void optimize(final double[] startPosition) {
    final int cmaesLambda = this.cmaesLambda;
    final double[] cmaesInputSigma = this.cmaesInputSigma;
    final int cmaesMaxIterations = this.cmaesMaxIterations;
    final int cmaesCheckFeasableCount = this.cmaesCheckFeasableCount;
    final ConvergenceChecker<PointValuePair> cmaesChecker = this.cmaesChecker;
    final BeamCenterRefinement function = this;
    Job job = new Job("Beam Position Refinement") {
        @Override
        protected IStatus run(IProgressMonitor monitor) {

            function.setInitPeaks(initPeaks);
            function.setMonitor(monitor);

            final double[] lB = new double[] { startPosition[0] - 20, startPosition[1] - 20 };
            final double[] uB = new double[] { startPosition[0] + 20, startPosition[1] + 20 };
            CMAESOptimizer beamPosOptimizer = new CMAESOptimizer(cmaesMaxIterations, 0.0, true, 0,
                    cmaesCheckFeasableCount, new Well19937a(), false, cmaesChecker);
            final PointValuePair result = beamPosOptimizer.optimize(new MaxEval(cmaesMaxIterations),
                    new ObjectiveFunction(function), GoalType.MAXIMIZE,
                    new CMAESOptimizer.PopulationSize(cmaesLambda), new CMAESOptimizer.Sigma(cmaesInputSigma),
                    new SimpleBounds(lB, uB), new InitialGuess(startPosition));

            final double[] newBeamPosition = result.getPoint();
            logger.info("Optimiser terminated at beam position ({}, {}) with the value {}",
                    new Object[] { newBeamPosition[0], newBeamPosition[1], result.getValue() });
            Display.getDefault().syncExec(new Runnable() {
                @Override
                public void run() {
                    ((IDiffractionMetadata) dataset.getMetadata()).getDetector2DProperties()
                            .setBeamCentreCoords(newBeamPosition);
                }
            });

            return Status.OK_STATUS;
        }
    };
    job.schedule();
}

From source file:playground.thibautd.initialdemandgeneration.socnetgensimulated.framework.ModelIterator.java

public SocialNetwork iterateModelToTarget(final ModelRunner runner, final Thresholds initialThresholds) {
    final MultivariateOptimizer optimizer = new CMAESOptimizer(maxIterations, 1E-9, true, 3, 50,
            new MersenneTwister(42), false, new Convergence());

    final double x = initialThresholds.getPrimaryThreshold();
    final double y = initialThresholds.getSecondaryReduction();

    final PointValuePair result = optimizer.optimize(GoalType.MINIMIZE, new MaxEval(maxIterations),
            new InitialGuess(new double[] { x, y }), new ObjectiveFunction(new Function(runner)),
            new CMAESOptimizer.Sigma(new double[] { 5, 500 }), new CMAESOptimizer.PopulationSize(7),
            new SimpleBounds(new double[] { Double.NEGATIVE_INFINITY, 0 }, // lower bounds: constrain secondary reduction to be positive
                    new double[] { Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY }) // upper bounds
    );//from   www.j  av a2 s  .co  m

    final Thresholds bestThresholds = new Thresholds(result.getPoint()[0], result.getPoint()[1]);
    final SocialNetwork bestSn = generate(runner, bestThresholds);

    log.info("best social network found for thresholds: " + bestThresholds);

    return bestSn;
}

From source file:uk.ac.diamond.scisoft.analysis.diffraction.FittingUtils.java

/**
 * @param opt/*from   w w  w  .  j a v a 2s .co  m*/
 * @param n
 * @return optimizer
 */
public static MultivariateOptimizer createOptimizer(Optimizer opt, int n) {
    switch (opt) {
    case BOBYQA:
        return new BOBYQAOptimizer(n + 2);
    case CMAES:
        return new CMAESOptimizer(MAX_ITER, 0., true, 0, 10,
                seed == null ? new Well19937c() : new Well19937c(seed), false,
                new SimplePointChecker<PointValuePair>(REL_TOL, ABS_TOL));
    case Simplex:
    default:
        return new SimplexOptimizer(new SimplePointChecker<PointValuePair>(REL_TOL * 1e3, ABS_TOL * 1e3));
    }
}

From source file:uk.ac.diamond.scisoft.analysis.optimize.ApacheOptimizer.java

private MultivariateOptimizer createOptimizer() {
    SimplePointChecker<PointValuePair> checker = new SimplePointChecker<PointValuePair>(REL_TOL, ABS_TOL);
    switch (optimizer) {
    case CONJUGATE_GRADIENT:
        return new NonLinearConjugateGradientOptimizer(Formula.POLAK_RIBIERE, checker);
    case BOBYQA:/*from  w  w  w .j  a va2s.  co  m*/
        return new BOBYQAOptimizer(n + 2);
    case CMAES:
        return new CMAESOptimizer(MAX_ITER, 0., true, 0, 10,
                seed == null ? new Well19937c() : new Well19937c(seed), false,
                new SimplePointChecker<PointValuePair>(REL_TOL, ABS_TOL));
    case POWELL:
        return new PowellOptimizer(REL_TOL, ABS_TOL, checker);
    case SIMPLEX_MD:
    case SIMPLEX_NM:
        return new SimplexOptimizer(checker);
    default:
        throw new IllegalStateException("Should not be called");
    }
}

From source file:uk.ac.diamond.scisoft.ncd.calibration.rcp.views.MultivariateFunctionWithMonitor.java

public void optimize(final double[] startPosition) {
    final int cmaesLambda = this.cmaesLambda;
    final double[] cmaesInputSigma = this.cmaesInputSigma;
    final int cmaesMaxIterations = this.cmaesMaxIterations;
    final int cmaesCheckFeasableCount = this.cmaesCheckFeasableCount;
    final ConvergenceChecker<PointValuePair> cmaesChecker = this.cmaesChecker;
    final MultivariateFunctionWithMonitor function = this;
    Job job = new Job(jobName) {
        @Override//from w  w w. j av a 2 s  .c  o m
        protected IStatus run(IProgressMonitor monitor) {

            function.setInitPeaks(initPeaks);
            function.setMonitor(monitor);

            CMAESOptimizer beamPosOptimizer = new CMAESOptimizer(cmaesMaxIterations, 0.0, true, 0,
                    cmaesCheckFeasableCount, new Well19937a(), false, cmaesChecker);
            final PointValuePair res = beamPosOptimizer.optimize(new MaxEval(cmaesMaxIterations),
                    new ObjectiveFunction(function), GoalType.MAXIMIZE,
                    new CMAESOptimizer.PopulationSize(cmaesLambda), new CMAESOptimizer.Sigma(cmaesInputSigma),
                    SimpleBounds.unbounded(2), new InitialGuess(startPosition));
            final double[] newBeamPos = res.getPoint();
            logger.info("Optimiser terminated at beam position ({}, {}) with the value {}",
                    new Object[] { newBeamPos[0], newBeamPos[1], res.getValue() });
            // Run calculation with optimised beam center to update UI 
            function.value(newBeamPos);

            return Status.OK_STATUS;
        }
    };
    job.schedule();
}

From source file:uk.ac.diamond.scisoft.ncd.core.data.plots.GuinierPlotData.java

public Object[] getGuinierPlotParameters(IDataset data, IDataset axis) {
    Dataset guinierData = getSaxsPlotDataset(data, axis);
    Dataset guinierAxis = getSaxsPlotAxis(axis);
    /*      int sliceSize = data.getSize() / 10;
          int[] step = new int[] {sliceSize};
          IndexIterator dataSliceIter = guinierData.getSliceIterator(null, null, null);
          IndexIterator axisSliceIter = guinierAxis.getSliceIterator(null, null, null);
          // Iterate over data slices/*  w  w  w. j  a  va  2s .  c  o m*/
          int[] minPosition = new int[] {-1};
          double minError = Double.MAX_VALUE;
          double slope = Double.NaN;
          double intercept = Double.NaN;
          Map<DoublePoint, Double> clusterInputMap = new HashMap<DoublePoint, Double>();
          while (dataSliceIter.hasNext() && axisSliceIter.hasNext()) {
             SimpleRegression regression = new SimpleRegression();
             int[] pos = dataSliceIter.getPos();
             // Iterate over data values for linear regression
             IndexIterator dataIter = new SliceIterator(
       guinierData.getShape(),
       guinierData.getSize(),
       pos, step);
             pos = axisSliceIter.getPos();
             IndexIterator axisIter = new SliceIterator(
       guinierAxis.getShape(),
       guinierAxis.getSize(),
       pos, step);
             int points = 0;
             while (dataIter.hasNext() && axisIter.hasNext()) {
    double dataVal = guinierData.getDouble(dataIter.getPos());
    double axisVal = guinierAxis.getDouble(axisIter.getPos());
    regression.addData(axisVal, dataVal);
    points++;
             }
             if (points == sliceSize) {
    regression.regress();
    double err = regression.getMeanSquareError();
    if (err < minError) {
       minError = err;
       minPosition = Arrays.copyOf(pos, pos.length);
       slope = regression.getSlope();
       intercept = regression.getIntercept();
       double I0 = Math.exp(intercept);
       double Rg = Math.sqrt(-3.0*slope);
       System.out.println("    Min Pos : " + Arrays.toString(minPosition));
       System.out.println("   Min Error : " + Double.toString(minError));
       System.out.println("   Slope : " + Double.toString(slope) + "   Intercept : " + Double.toString(intercept));
       System.out.println("   I(0) : " + Double.toString(I0) + "   Rg : " + Double.toString(Rg));
       clusterInputMap.put(new DoublePoint(minPosition), minError);
    }
             } else {
    break;
             }
          }
                  
          DBSCANClusterer<DoublePoint> clusterer = new DBSCANClusterer<DoublePoint>(5, 5);
          List<Cluster<DoublePoint>> clusterResults = clusterer.cluster(clusterInputMap.keySet());
            
          // output the clusters
          for (int i = 0; i < clusterResults.size(); i++) {
              System.out.println("Cluster " + i);
             double[] minPoint = null;
             double minVal = Double.MAX_VALUE;
              for (DoublePoint point : clusterResults.get(i).getPoints()) {
      System.out.println(Arrays.toString(point.getPoint()));
      Double val = clusterInputMap.get(point);
      if (val < minVal) {
         minVal = val;
         minPoint = Arrays.copyOf(point.getPoint(), point.getPoint().length);
      }
      minVal = (val < minVal ? val : minVal);
              }
               System.out.println("Min cluster point : " + Arrays.toString(minPoint));
               System.out.println("Min cluster value : " + Double.toString(minVal));
              System.out.println();
          }
                  
    */ ConvergenceChecker<PointValuePair> cmaesChecker = new SimplePointChecker<PointValuePair>(1e-6, 1e-8);
    RandomDataGenerator rnd = new RandomDataGenerator();
    CMAESOptimizer optimizer = new CMAESOptimizer(cmaesMaxIterations, 0.0, true, 0, cmaesCheckFeasableCount,
            rnd.getRandomGenerator(), false, cmaesChecker);
    GuinierLineFitFunction function = new GuinierLineFitFunction(guinierData, guinierAxis);

    Amount<Dimensionless> I0 = Amount.valueOf(Double.NaN, Double.NaN, Dimensionless.UNIT);
    Amount<Dimensionless> Rg = Amount.valueOf(Double.NaN, Double.NaN, Dimensionless.UNIT);
    double[] qvals = new double[] { Double.NaN, Double.NaN };

    double q0 = guinierAxis.getDouble(0);
    double qMin = guinierAxis.getDouble(1);
    double qMax = guinierAxis.getDouble(guinierAxis.getSize() - 1);
    double[] startPosition = new double[] { guinierAxis.getDouble(0),
            guinierAxis.getDouble(GuinierLineFitFunction.MIN_POINTS) };
    double[] cmaesInputSigma = new double[] { (qMin - q0) * 0.1, qMax * 0.1 };
    try {
        final PointValuePair res = optimizer.optimize(new MaxEval(cmaesMaxIterations),
                new ObjectiveFunction(function), GoalType.MAXIMIZE,
                new CMAESOptimizer.PopulationSize(cmaesLambda), new CMAESOptimizer.Sigma(cmaesInputSigma),
                new SimpleBounds(new double[] { q0, q0 }, new double[] { qMin, qMax }),
                new InitialGuess(startPosition));

        qvals = res.getPoint();
        function.value(qvals);
        I0 = getI0(function.regression);
        Rg = getRg(function.regression);

        System.out.println("Final Result");
        String msg = StringUtils.join(new String[] { "   I(0) ", I0.toString(), "   Rg ", Rg.toString() },
                " : ");
        System.out.println(msg);
        msg = StringUtils.join(new String[] { "Slice", ArrayUtils.toString(res.getPoint()), "R",
                Double.toString(function.regression.getR()) }, " : ");
        System.out.println(msg);

        /*         // Run Monte-Carlo simulation to generate error estimates Rg values 
                 //double finalR = function.regression.getR();
                 int maxSample = 10000;
                 int minSample = 10;
                 int totalSample = 100000;
                 int counter = 0;
                 int totalCounter = 0;
                 GuinierLineFitFunction mcFunction = new GuinierLineFitFunction(guinierData, guinierAxis);
                 DescriptiveStatistics statsR = new DescriptiveStatistics();
                 List<Pair<Double, Amount<Dimensionless>>> listI0 = new ArrayList<Pair<Double,Amount<Dimensionless>>>();
                 List<Pair<Double, Amount<Dimensionless>>> listRg = new ArrayList<Pair<Double,Amount<Dimensionless>>>();
                 while ((counter < maxSample && totalCounter < totalSample)
                       || (counter < minSample && totalCounter >= totalSample)) {
                    double q1 = rnd.nextUniform(q0, qMin);
                    double q2 = rnd.nextUniform(q0, qMax);
                    if (!(q2 > q1)) {
                       continue;
                    }
                    totalCounter++;
                     mcFunction.value(new double[] {q1, q2});
                     double tmpR = Math.abs(mcFunction.regression.getR());
                     //boolean equalsR = Precision.equalsWithRelativeTolerance(tmpR, finalR, 0.1); 
                     if (!(Double.isNaN(tmpR) || Double.isInfinite(tmpR))) {
        statsR.addValue(tmpR);
        Amount<Dimensionless> tmpI0 = getI0(mcFunction.regression);
        Amount<Dimensionless> tmpRg = getRg(mcFunction.regression);
         if (Double.isNaN(tmpI0.getEstimatedValue()) || Double.isInfinite(tmpI0.getEstimatedValue()) ||
               Double.isNaN(tmpRg.getEstimatedValue()) || Double.isInfinite(tmpRg.getEstimatedValue())) {
            continue;
         }
        listI0.add(new Pair<Double, Amount<Dimensionless>>(tmpR, tmpI0));
        listRg.add(new Pair<Double, Amount<Dimensionless>>(tmpR, tmpRg));
        counter++;
                     }
                 }
                         
                 double threshold = statsR.getPercentile(90);
                 //double threshold = 0.95*statsR.getMax();
                 SummaryStatistics statsI0 = new SummaryStatistics();
                 SummaryStatistics statsRg = new SummaryStatistics();
                 for (Pair<Double, Amount<Dimensionless>> tmpVal : listRg) {
                    if (tmpVal.getFirst() > threshold) {
        statsRg.addValue(tmpVal.getSecond().getEstimatedValue());
                    }
                 }
                 for (Pair<Double, Amount<Dimensionless>> tmpVal : listI0) {
                    if (tmpVal.getFirst() > threshold) {
        statsI0.addValue(tmpVal.getSecond().getEstimatedValue());
                    }
                 }
                         
                 double meanI0 = statsI0.getMean();
                 double stdI0 = statsI0.getStandardDeviation();
                 I0 = Amount.valueOf(meanI0, stdI0, Dimensionless.UNIT);
                         
                 double meanRg = statsRg.getMean();
                 double stdRg = statsRg.getStandardDeviation();
                 Rg = Amount.valueOf(meanRg, stdRg, Dimensionless.UNIT);
                         
                 String msg = StringUtils.join(new String[] {
                       "Monte-Carlo Rg", Rg.toString()
                       },
                       " : ");
                 System.out.println(msg);
        */
    } catch (MaxCountExceededException e) {
        System.out.println("Maximum counts exceeded");
        return null;
    }
    return new Object[] { I0, Rg, qvals[0], qvals[1] };
}

From source file:uk.ac.diamond.scisoft.ncd.core.data.plots.LogLogPlotData.java

public double[] getPorodPlotParameters(IDataset data, IDataset axis) {
    Dataset loglogData = getSaxsPlotDataset(data, axis);
    Dataset loglogAxis = getSaxsPlotAxis(axis);
    CMAESOptimizer optimizer = new CMAESOptimizer(cmaesMaxIterations, 0.0, true, 0, cmaesCheckFeasableCount,
            new Well19937a(), false, cmaesChecker);
    PorodLineFitFunction function = new PorodLineFitFunction(loglogData, loglogAxis);

    int dataSize = loglogAxis.getSize();
    double qMax = Math.pow(10.0, loglogAxis.getDouble(dataSize - 1));
    double i0 = Math.pow(10.0, loglogData.getDouble(0));
    double[] startPosition = new double[] { i0, 1.0 / qMax, 4.0 };
    double[] cmaesInputSigma = new double[] { 0.1 * i0, 0.1 / qMax, 0.1 };
    double[] lb = new double[] { 0.0, 0.0, 0.0 };
    double[] ub = new double[] { Double.MAX_VALUE, Double.MAX_VALUE, Double.MAX_VALUE };
    double[] params = new double[] { Double.NaN, Double.NaN, Double.NaN };
    try {/*from   w  w  w  . j  a  v a  2s.  c o m*/
        final PointValuePair res = optimizer.optimize(new MaxEval(cmaesMaxIterations),
                new ObjectiveFunction(function), GoalType.MAXIMIZE,
                new CMAESOptimizer.PopulationSize(cmaesLambda), new CMAESOptimizer.Sigma(cmaesInputSigma),
                new SimpleBounds(lb, ub), new InitialGuess(startPosition));

        params = res.getPoint();
        double solvation = params[0];
        double correlation = params[1];
        double exponent = params[2];

        double rms = Math.exp(-function.value(params));

        System.out.println();
        System.out.println("Final Result");
        String msg = StringUtils.join(new String[] { "   solvation ", Double.toString(solvation),
                "   correlation ", Double.toString(correlation), "   exponent ", Double.toString(exponent) },
                " : ");
        System.out.println(msg);
        msg = StringUtils.join(new String[] { "RMS", Double.toString(rms) }, " : ");
        System.out.println(msg);
        System.out.println();

    } catch (MaxCountExceededException e) {
        params = new double[] { Double.NaN, Double.NaN, Double.NaN };
    }
    return params;
}

From source file:uk.ac.diamond.scisoft.ncd.core.data.plots.PorodPlotData.java

public SimpleRegression getPorodPlotParameters(IDataset data, IDataset axis) {
    Dataset porodData = getSaxsPlotDataset(data, axis);
    Dataset porodAxis = getSaxsPlotAxis(axis);
    CMAESOptimizer optimizer = new CMAESOptimizer(cmaesMaxIterations, 0.0, true, 0, cmaesCheckFeasableCount,
            new Well19937a(), false, cmaesChecker);
    PorodLineFitFunction function = new PorodLineFitFunction(porodData, porodAxis);

    int dataSize = porodAxis.getSize();
    double q0 = porodAxis.getDouble(0);
    double qMax = porodAxis.getDouble(dataSize - 1);
    double[] startPosition = new double[] { porodAxis.getDouble(dataSize - PorodLineFitFunction.MIN_POINTS - 1),
            porodAxis.getDouble(dataSize - 1) };
    double[] cmaesInputSigma = new double[] { qMax * 0.1, qMax * 0.1 };
    try {/*from  ww  w  .  j  ava2 s  .  co  m*/
        final PointValuePair res = optimizer.optimize(new MaxEval(cmaesMaxIterations),
                new ObjectiveFunction(function), GoalType.MAXIMIZE,
                new CMAESOptimizer.PopulationSize(cmaesLambda), new CMAESOptimizer.Sigma(cmaesInputSigma),
                new SimpleBounds(new double[] { q0, q0 }, new double[] { qMax, qMax }),
                new InitialGuess(startPosition));

        function.value(res.getPoint());
        double slope = function.regression.getSlope();
        Amount<Dimensionless> c4 = getC4(function.regression);

        System.out.println("Final Result");
        String msg = StringUtils
                .join(new String[] { "   c4 ", c4.toString(), "   slope ", Double.toString(slope) }, " : ");
        System.out.println(msg);
        msg = StringUtils.join(new String[] { "Slice", ArrayUtils.toString(res.getPoint()), "R",
                Double.toString(function.regression.getR()) }, " : ");
        System.out.println(msg);
    } catch (MaxCountExceededException e) {
        return null;
    }
    return function.regression;
}