leastsq10.java

Go to the documentation of this file.

 
import java.util.*;
import java.lang.Math;
import conopt.*;

public class leastsq10 {
   public static void main(String argv[]) {
      System.loadLibrary("conoptjni4");
 
      String name = "leastsq10";
 
      Conopt conopt = new Conopt(name);
      LeastSq10ModelData model = new LeastSq10ModelData(700, 500);
      MsgHandler msghdlr = new MsgHandler(name);
 
      model.buildModel();
 
      conopt.loadModel(model);
      conopt.setMessageHandler(msghdlr);
 
      // try to set the license using the environment variables
      try {
         int license_int_1 = Integer.parseInt(System.getenv("CONOPT_LICENSE_INT_1"));
         int license_int_2 = Integer.parseInt(System.getenv("CONOPT_LICENSE_INT_2"));
         int license_int_3 = Integer.parseInt(System.getenv("CONOPT_LICENSE_INT_3"));
         String license_text = System.getenv("CONOPT_LICENSE_TEXT");
 
         conopt.setLicense(license_int_1, license_int_2, license_int_3,
               license_text);
      } catch (Exception e) {
         System.out.println("Unable to set license: " + e.getMessage());
      }
 
      conopt.solve();
 
      conopt.printStatus();
 
      std s = new std();
      int retcode = s.checkSolve(name, conopt.modelStatus(), conopt.solutionStatus(),
            conopt.objectiveValue(), 19.4443, 0.001);
 
      msghdlr.close();
 
      System.exit(retcode);
   }
}
 
class LeastSq10ModelData extends ModelData {
   private int seed = 12359;
 
   // Model Parameters
   public double[] A;       // length = nobs * dimx
   public double[] B;       // length = nobs * dimx
   public double[] C;       // length = nobs * dimx
   public double[] Obs;     // length = nobs
   public int nobs;
   public int dimx;
 
   // Declaring the variable and constraint indices
   public int[] varx;           // length = dimx
   public int[] varres;         // length = nobs
   public int consobj;
   public int[] consresidual;   // length = nobs
 
 
   public LeastSq10ModelData(int numobs, int dimensionx) {
      super();
 
      this.nobs = numobs;
      this.dimx = dimensionx;
 
      this.A   = new double[nobs * dimx];
      this.B   = new double[nobs * dimx];
      this.C   = new double[nobs * dimx];
      this.Obs = new double[nobs];
 
      defineData();
   }
 
   // @brief Defines a pseudo random number between 0 and 1
   private double rndx() {
      seed = seed * 1027 + 25;
      int times = seed / 1_048_576;
      seed = seed - 1_048_576 * times;
      return (double) seed / 1_048_576.0;
   }
 
   // @brief Defines the data for the problem
   private void defineData() {
      final double Xtarg = -1.0;
      final double Noise = 1.0;
 
      int k = 0;
      for (int i = 0; i < nobs; ++i) {
         double O = 0.0;
         for (int j = 0; j < dimx; ++j) {
               A[k] = rndx();
               B[k] = rndx();
               O += A[k] * Xtarg + B[k] * (Xtarg * Xtarg);
               ++k;
         }
         Obs[i] = O + Noise * rndx();
      }
   }

   public void buildModel() {
      // adding the variables to the model
      List<Integer> varxList = new ArrayList<>(dimx);
      for (int i = 0; i < dimx; ++i) {
         int varidx = addVariable(-Conopt.Infinity, Conopt.Infinity, -0.8);
         varxList.add(varidx);
      }
 
      List<Integer> varresList = new ArrayList<>(nobs);
      for (int i = 0; i < nobs; ++i) {
         int varidx = addVariable(-Conopt.Infinity, Conopt.Infinity, 0.0);
         varresList.add(varidx);
      }
 
      this.varx   = varxList.stream().mapToInt(Integer::intValue).toArray();
      this.varres = varresList.stream().mapToInt(Integer::intValue).toArray();
 
 
 
      // adding the constraints to the model
      this.consresidual = new int[nobs];
      for (int i = 0; i < nobs; ++i) {
 
         int m = dimx + 1;                   // +1 for the residual var
         int[]    varidx = new int[m];
         double[] coeffs = new double[m];
         int[]    nlf    = new int[m];
 
         for (int j = 0; j < dimx; ++j) {
            varidx[j] = varx[j];
            coeffs[j] = 0.0;
            nlf[j]    = 1;
         }
 
         // Add residual var linearly with coefficient -1
         varidx[dimx] = varres[i];
         coeffs[dimx] = -1.0;
         nlf[dimx]    = 0;
 
         int considx = addConstraint(ConstraintType.Eq, Obs[i], varidx, coeffs, nlf);
         consresidual[i] = considx;
      }
 
      // Objective function
      double[] objCoeffs = new double[nobs];
      int[]    objNlf    = new int[nobs];
      for (int i = 0; i < nobs; ++i) {
         objCoeffs[i] = 1.0;
         objNlf[i]    = 1;    // residuals appear nonlinearly in the objective (r_i^2)
      }
 
      consobj = addConstraint(ConstraintType.Free, 0.0, varres, objCoeffs, objNlf);
 
      // setting the objective constraint
      setObjectiveElement(ObjectiveElement.Constraint, consobj);
 
      // setting the optimisation direction
      setOptimizationSense(Sense.Minimize);
 
      // setting the second derivative evaluation type
      setSDEvaluationType(SDEvaluationType.Constraint);
   }

   public double evaluateNonlinearTerm(double[] x, int rowno, boolean ignerr, int thread) {
      double g = 0;
      if (rowno == consobj) {
         double sum = 0.0;
         for (int i = 0; i < nobs; i++) {
            sum += Math.pow(x[varres[i]], 2);
         }
         g = sum;
      }
      else {
         int k = rowno * dimx;
         double sum = 0.0;
         for (int i = 0; i < dimx; i++) {
            sum += A[k] * x[varx[i]] + B[k] * Math.pow(x[varx[i]], 2);
            k++;
         }
         g = sum;
      }
 
      return g;
   }

   public void evaluateNonlinearJacobian(double[] x, double[] jac, int rowno, int[] jacnum, boolean ignerr, int thread) {
      assert x.length == jac.length;
 
      if (rowno == consobj) {
         for (int i = 0; i < nobs; i++) {
            jac[varres[i]] = 2 * x[varres[i]];
         }
      }
      else {
         int k = rowno * dimx;
         for (int i = 0; i < dimx; i++) {
            jac[varx[i]] = A[k] + 2 * B[k] * x[varx[i]];
            k++;
         }
      }
   }

   public void initDirectionalSDEval(double[] x, double[] dx, int[] rowlist, int numthread, boolean newpoint) {
      //   Constraint row with b(i,j)*x(j)**2 as only nonlinear term.
      //     Compute the directional term 2*b(i,j)*dx(j) and keep for later
 
      int k = 0;
      for (int i = 0; i < nobs; i++) {
         for (int j = 0; j < dimx; j++) {
            C[k] = 2.0 * B[k] * dx[varx[j]];
            k++;
         }
      }
   }
 

   public void evaluateDirectionalSD(double[] x, double[] dx, double[] d2g, int rowno, int[] jacnum, int thread) {
      assert x.length == d2g.length;
      assert dx.length == d2g.length;
 
      if (rowno == nobs) {
         /*  Objective row: sum(i, res(i)**2 )                                 */
 
         for (int i = 0; i < dimx; i++)
            d2g[varx[i]] = 0;
         for (int i = 0; i < nobs; i++)
            d2g[varres[i]] = 2.0 * dx[varres[i]];
      }
      else {
         /*  Constraint row with b(i,j)*x(j)**2 as only nonlinear term         */
 
         int k = rowno * dimx;
         for (int i = 0; i < dimx; i++) {
            d2g[varx[i]] = C[k];
            k++;
         }
         for (int i = 0; i < nobs; i++)
            d2g[varres[i]] = 0;
      }
   }
}