leastsq10.cpp

Go to the documentation of this file.

 
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <iostream>
#include "conopt.hpp"
 
int seed = 12359; /* seed for random number generator */
 
class Leastsq10_ModelData : public ConoptModelData
{
public:
   std::vector<double> A;
   std::vector<double> B;
   std::vector<double> C;
   std::vector<double> Obs;
   int nobs;
   int dimx;
 
   std::vector<int> varx;
   std::vector<int> varres;
 
   int consobj;
   std::vector<int> consresidual;
 
   Leastsq10_ModelData(int numobs, int dimensionx) : nobs(numobs), dimx(dimensionx)
   {
      A.resize(nobs * dimx);
      B.resize(nobs * dimx);
      C.resize(nobs * dimx);
      Obs.resize(nobs);
      defineData();
   }

   double rndx()
   {
      int times;
      double result;
      seed = seed * 1027 + 25;
      times = seed / 1048576;
      seed = seed - 1048576 * times;
      result = static_cast<double>(seed);
      result = result / 1048576.0;
      return result;
   }

   void defineData()
   {
      double Xtarg = -1.0;
      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] * std::pow(Xtarg, 2);
            ++k;
         }
         Obs[i] = O + Noise * rndx();
      }
   }

   void buildModel()
   {
      /*                                                               */
      /*   Information about Variables:                                */
      /*   Default: Lower = -Inf, Curr = 0, and Upper = +inf.          */
      /*   Default: the status information in Vsta is not used.        */
      /*                                                               */
      for (int i = 0; i < dimx; ++i)
      {
         int varidx = addVariable(-Conopt::Infinity, Conopt::Infinity, -0.8);
         varx.push_back(varidx);
      }
      for (int i = 0; i < nobs; ++i)
      {
         int varidx = addVariable(-Conopt::Infinity, Conopt::Infinity, 0.0);
         varres.push_back(varidx);
      }
 
      /*   Information about Constraints:                              */
      /*   Default: Rhs = 0                                            */
      /*   Default: the status information in Esta and the function    */
      /*            value in FV are not used.                          */
      /*   Default: Type: There is no default.                         */
      /*            0 = Equality,                                      */
      /*            1 = Greater than or equal,                         */
      /*            2 = Less than or equal,                            */
      /*            3 = Non binding.                                   */
      /*                                                               */
      for (int i = 0; i < nobs; ++i)
      {
         std::vector<int> varidx(varx.begin(), varx.end());
         std::vector<double> coeffs(dimx, 0.0);
         std::vector<int> nlf(dimx, 1);
 
         // ADD THE OUTPUT RESIDUAL VARIABLE
         varidx.push_back(varres[i]); // residual variable index
         coeffs.push_back(-1.0);      // linear coefficient
         nlf.push_back(0);            // linear
 
         int considx = addConstraint(ConoptConstraintType::Eq, Obs[i], varidx, coeffs, nlf);
         consresidual.push_back(considx);
      }
 
      /*   Constraint nobs (Objective)                                 */
      /*   Rhs = 0.0 and type Non binding                              */
      /*                                                               */
      std::vector<int> objVarIdx(varres.begin(), varres.end());
      std::vector<double> objCoeffs(nobs, 1.0);
      std::vector<int> objNlf(nobs, 1);
 
      consobj = addConstraint(ConoptConstraintType::Free, 0.0, objVarIdx, objCoeffs, objNlf);
 
      setObjectiveElement(ConoptObjectiveElement::Constraint, consobj);
      setOptimizationSense(ConoptSense::Minimize);
 
      /* setting the second derivative evaluation type */
      setSDEvaluationType(ConoptSDEvaluationType::Constraint);
   }

   int FDEval(const double x[], double *g, double jac[], int rowno, const int jacnum[], int mode,
         int ignerr, int *errcnt, int numvar, int numjac, int thread) override
   {
      if (rowno == consobj)
      {
         /*                                                                   */
         /*   Mode = 1 or 3. Function value: g = P * Out                      */
         /*                                                                   */
         if (mode == 1 || mode == 3)
         {
            double sum = 0.0;
            for (int i = 0; i < nobs; i++)
            {
               sum += pow(x[varres[i]], 2);
            }
            *g = sum;
         }
         /*                                                                   */
         /*   Mode = 2 or 3: Derivative values:                               */
         /*                                                                   */
         if (mode == 2 || mode == 3)
         {
            for (int i = 0; i < nobs; i++)
            {
               jac[varres[i]] = 2 * x[varres[i]];
            }
         }
      }
      /*
      Row 0 to nobs-1: The observation definitions:                       */
      else
      {
         /*
         Mode = 1 or 3: Function value - x-part only                        */
 
         if (mode == 1 || mode == 3)
         {
            int k = rowno * dimx;
            double sum = 0.0;
            for (int i = 0; i < dimx; i++)
            {
               sum += A[k] * x[varx[i]] + B[k] * pow(x[varx[i]], 2);
               k++;
            }
            *g = sum;
         }
         /*
         Mode = 2 or 3: Derivatives                                         */
 
         if (mode == 2 || mode == 3)
         {
            int k = rowno * dimx;
            for (int i = 0; i < dimx; i++)
            {
               jac[varx[i]] = A[k] + 2 * B[k] * x[varx[i]];
               k++;
            }
         }
      }
 
      return 0;
   }

   int SDDirIni(const double x[], const double dx[], const int rowlist[], int listsize,
         int numthread, int newpt, int *nodrv, int numvar) override
   {
      /*  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++;
         }
      }
      return 0;
   }

   int SDDir(const double x[], const double dx[], double d2g[], int rowno, const int jacnum[],
         int *nodrv, int numvar, int numjac, int thread) override
   {
      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;
      }
 
      return 0;
   }
};
 
#include "std.cpp"

int main(int argc, char **argv)
{
   int COI_Error = 0;
 
   // getting the program name from the executable path
   std::string pname = getProgramName(argv[0]);
 
   // initialising the Conopt Object
   Conopt conopt(pname);
   Leastsq10_ModelData modeldata(700, 500);
   Tut_MessageHandler msghandler(pname);
 
   // adding the message handler to the conopt interface
   conopt.setMessageHandler(msghandler);
 
   // building the model
   modeldata.buildModel();
 
   // loading the model in the conopt object
   conopt.loadModel(modeldata);
 
#if defined(CONOPT_LICENSE_INT_1) && defined(CONOPT_LICENSE_INT_2) && defined(CONOPT_LICENSE_INT_3) && defined(CONOPT_LICENSE_TEXT)
   std::string license = CONOPT_LICENSE_TEXT;
   COI_Error += conopt.setLicense(CONOPT_LICENSE_INT_1, CONOPT_LICENSE_INT_2, CONOPT_LICENSE_INT_3, license);
#endif
 
   if (COI_Error)
      cpp_log(
            conopt, "Skipping COI_Solve due to license error. COI_Error = " + std::to_string(COI_Error), COI_Error);
 
   COI_Error = conopt.solve(); /* Optimize                      */
 
   // checking the status and objective value
   if (!(conopt.solutionStatus() == 1 && conopt.modelStatus() == 2))
   {
      cpp_log(conopt, "Incorrect Model or Solver Status", -1);
   }
   else if (fabs(conopt.objectiveValue() - 19.4443) > 0.001)
   {
      cpp_log(conopt, "Incorrect objective returned", -1);
   }
 
   // printing the final status of the optimisation
   conopt.printStatus();
 
   cpp_log(conopt, "Successful Solve", COI_Error);
}