leastsq2.py

Go to the documentation of this file.

 
import sys
import os
import math
 
import conopt as co
 
sys.path.append('../common/')
import std
 
 
class LeastSqModelData(co.ModelData):
   def __init__(self, numobs: int, dimensionx: int):
      self.seed = 12359  # seed for random number generator
 
      self.nobs = numobs
      self.dimx = dimensionx
 
      self.A = [0] * (self.nobs * self.dimx)
      self.B = [0] * (self.nobs * self.dimx)
      self.Obs = [0] * self.nobs
 
      self.varx = []
      self.varres = []
      self.consresidual = []
 
      self.define_data()
 
      super().__init__()
 
   def rndx(self):
      """
      Defines a pseudo random number between 0 and 1
 
      NOTE: it would be possible to use random module to generate the random
      numbers. We have written our own random number generator to be consistent
      with the original Fortran example.
      """
      self.seed = self.seed * 1027 + 25
      times = int(self.seed / 1048576)
      self.seed = self.seed - 1048576 * times
      return self.seed / 1048576.0
 
   def define_data(self):
      """
      Defines the data for the problem
      """
      Xtarg = -1.0
      Noise = 1.0
      k = 0
 
      for i in range(self.nobs):
         O = 0
         for j in range(self.dimx):
            self.A[k] = self.rndx()
            self.B[k] = self.rndx()
            O += self.A[k] * Xtarg + self.B[k] * math.pow(Xtarg, 2)
            k = k + 1
         self.Obs[i] = O + Noise * self.rndx()
 
   def buildModel(self):
      """
      adding the variables and constraints to the model
      @ingroup PYTHON1THREAD_LEASTSQ2
      """
      # adding the variables to the model
      for i in range(self.dimx):
         varidx = self.addVariable(
            -co.Conopt.Infinity, co.Conopt.Infinity, -0.8
         )
         self.varx.append(varidx)
 
      for i in range(self.nobs):
         varidx = self.addVariable(-co.Conopt.Infinity, co.Conopt.Infinity, 0.0)
         self.varres.append(varidx)
 
      # adding the constraints to the model
      for i in range(self.nobs):
         varidx = self.varx.copy()
         coeffs = [0.0] * self.dimx
         nlf = [1] * self.dimx
 
         # ADD THE OUTPUT RESIDUAL VARIABLE
         varidx.append(self.varres[i])  # residual variable index
         coeffs.append(-1.0)  # linear coefficient
         nlf.append(0)  # linear
 
         considx = self.addConstraint(
            co.ConstraintType_Eq, self.Obs[i], varidx, coeffs, nlf
         )
         self.consresidual.append(considx)
 
      # Constraint nobs (Objective)
      objVarIdx = self.varres.copy()
      objCoeffs = [1.0] * self.nobs
      objNlf = [1] * self.nobs
 
      self.consobj = self.addConstraint(
         co.ConstraintType_Free, 0.0, objVarIdx, objCoeffs, objNlf
      )
 
      # setting the objective constraint
      self.setObjectiveElement(co.ObjectiveElement_Constraint, self.consobj)
 
      # setting the optimisation direction
      self.setOptimizationSense(co.Sense_Minimize)
 
      # We assume that this is a Large Residual model, where the most
      # important 2nd derivatives are those that appear directly in the
      # objective.
      # We will therefore only define the 2nd derivatives corresponding to
      # the objective constraint where the 2nd derivatives is 2*I.
      # CONOPT will by default complain that some nonlinear variables do
      # not appear in the Hessian. We will therefore define one dummy element
      # on the diagonal of the Hessian for each of the nonlinear X-variables
      # and give it the numerical value 0.
 
      # setting the structure of the hessian
      Qrow = [0] * (self.dimx + self.nobs)
      Qcol = [0] * (self.dimx + self.nobs)
 
      for i in range(self.dimx + self.nobs):
         Qrow[i] = i
         Qcol[i] = i
 
      self.setSDLagrangianStructure(Qrow, Qcol)
 
   def evaluateNonlinearTerm(self, x, rowno, ignerr, thread):
      """
      @copydoc conopt.ModelData.evaluateNonlinearTerm
      @ingroup PYTHON1THREAD_LEASTSQ2
      """
      g = 0
      if rowno == self.consobj:
         sum = 0.0
         for i in range(self.nobs):
            sum += pow(x[self.varres[i]], 2)
         g = sum
 
      else:
         k = rowno * self.dimx
         sum = 0.0
         for i in range(self.dimx):
            sum += self.A[k] * x[self.varx[i]] + self.B[k] * pow(
               x[self.varx[i]], 2
            )
            k += 1
         g = sum
 
      return g
 
   def evaluateNonlinearJacobian(self, x, rowno, jacnum, ignerr, thread):
      """
      @copydoc conopt.ModelData.evaluateNonlinearJacobian
      @ingroup PYTHON1THREAD_LEASTSQ2
      """
 
      jac = []
      if rowno == self.consobj:
         for i in range(self.nobs):
            jac.append(2 * x[self.varres[i]])
 
      else:
         k = rowno * self.dimx
         for i in range(self.dimx):
            jac.append(self.A[k] + 2 * self.B[k] * x[self.varx[i]])
            k += 1
 
      return jac
 
   def evaluateSDLagrangian(self, x, u, hessianrow, hessiancol):
      """
      @copydoc conopt.ModelData.evaluateSDLagrangian
      @ingroup PYTHON1THREAD_LEASTSQ2
      """
      numhessian = len(hessianrow)
      hessian = [0 for i in range(numhessian)]
      for i in range(self.dimx):
         hessian[self.varx[i]] = 0.0
      for i in range(self.nobs):
         hessian[self.varres[i]] = 2.0 * u[self.nobs]
 
      return hessian
 
 
if __name__ == '__main__':
   name = os.path.basename(__file__)[:-3]
 
   conopt = co.Conopt(name)
   model = LeastSqModelData(700, 500)
   msghdlr = std.TutMessageHandler(name)
 
   model.buildModel()
 
   conopt.loadModel(model)
   conopt.setMessageHandler(msghdlr)
 
   # getting the license variables
   license_int_1 = os.environ.get('CONOPT_LICENSE_INT_1', None)
   license_int_2 = os.environ.get('CONOPT_LICENSE_INT_2', None)
   license_int_3 = os.environ.get('CONOPT_LICENSE_INT_3', None)
   license_text = os.environ.get('CONOPT_LICENSE_TEXT', None)
   if (
      license_int_1 is not None
      and license_int_2 is not None
      and license_int_3 is not None
      and license_text is not None
   ):
      conopt.setLicense(
         int(license_int_1),
         int(license_int_2),
         int(license_int_3),
         license_text,
      )
 
   coi_error = conopt.solve()
 
   retcode = std.checkSolve(conopt, 19.4443, coi_error, 0.001)
 
   sys.exit(retcode)