mp_leastsq2.c

Go to the documentation of this file.

 
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include "coiheader.h"
#include "comdecl.h"
 
int seed = 12359;   /* seed for random number generator */
#define nobs 700
#define dimx 500
double A[nobs*dimx];
double B[nobs*dimx];
double Obs[nobs];
 
float rndx()
{
/*
   Defines a pseudo random number between 0 and 1    */
   int   times;
   float result;
   seed   = seed*1027+25;
   times  = seed / 1048576;
   seed   = seed - 1048576*times;
   result = (float)seed;
   result = result / 1048576;
/*   fprintf(stdout,"Seed=%d, result=%f8 \n",seed, result); */
   return result;
}
 
void defdata()
{
   int i,j,k;
   double Xtarg = -1.0;
   double Noise = 1.0;
   double O;
 
   k = 0;
   for ( i = 0; i < nobs; i++ ) {
      O = 0.0;
      for ( j = 0; j < dimx; j++ ) {
         A[k] = rndx();
         B[k] = rndx();
         O    = O + A[k] * Xtarg + B[k] * pow(Xtarg,2);
         k    = k + 1;
      }
      Obs[i] = O + Noise * rndx();
   }
}
 

int COI_CALLCONV Leastsq_ReadMatrix( double LOWER[], double CURR[], double UPPER[], int VSTA[], int TYPE[], double RHS[],
                            int ESTA[], int COLSTA[], int ROWNO[], double VALUE[],
                            int NLFLAG[], int NUMVAR, int NUMCON, int NUMNZ, void* USRMEM )
{
  int i, j, k;
 /*                                                               */
 /*   Information about Variables:                                */
 /*   Default: Lower = -Inf, Curr = 0, and Upper = +inf.          */
 /*   Default: the status information in Vsta is not used.        */
 /*                                                               */
   for ( i=0; i<dimx; i++ ) {
      CURR[i] = -0.8 ;
      }
 /*   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 ( i=0; i<nobs; i++ ) {
      RHS[i]  = Obs[i] ;
      TYPE[i] = 0;
      }
 
 /*   Constraint nobs (Objective)                                 */
 /*   Rhs = 0.0 and type Non binding                              */
 /*                                                               */
   TYPE[nobs] = 3    ;
 /*                                                               */
 /*   Information about the Jacobian. We use the standard method with */
 /*   Rowno, Value, Nlflag and Colsta and we do not use Colno.    */
 /*                                                               */
 /*   Colsta = Start of column indices (No Defaults):             */
 /*   Rowno  = Row indices                                        */
 /*   Value  = Value of derivative (by default only linear        */
 /*                                 derivatives are used)         */
 /*   Nlflag = 0 for linear and 1 for nonlinear derivative        */
 /*            (not needed for completely linear models)          */
 /*                                                               */
 /*  Indices                                                      */
 /*       x(i)  res(j)                                            */
 /*  j:    NL   L=1                                               */
 /*  obj:       NL                                                */
 
   k = 0;  /* counter for current nonzero   */
   for ( j=0; j<dimx; j++ ) {
      COLSTA[j] = k;
      for (i=0; i<nobs; i++ ) {
         ROWNO[k]  = i;
         NLFLAG[k] = 1;
         k         = k + 1;
      }
   }
   for (i=0; i<nobs; i++ ) {
      COLSTA[dimx+i] = k;
      ROWNO[k]  = i;
      NLFLAG[k] = 0;
      VALUE[k]  = 1.0;
      k         = k + 1;
      ROWNO[k]  = nobs;
      NLFLAG[k] = 1;
      k         = k + 1;
   }
   COLSTA[dimx+nobs] = k;
   return 0;
}
 

int COI_CALLCONV Leastsq_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, void* USRMEM )
{
   int i, j, k;
   double sum;
 /*                                                                   */
 /*   Row nobs: the objective function is nonlinear                      */
 /*                                                                   */
 
   if ( ROWNO == nobs ) {
 /*                                                                   */
 /*   Mode = 1 or 3. Function value: G = P * Out                      */
 /*                                                                   */
      if ( MODE == 1 || MODE == 3 ) {
         sum = 0.0;
         for (i=0; i<nobs; i++ ) {
            sum = sum + pow(X[dimx+i],2);
            }
         *G = sum;
         }
 /*                                                                   */
 /*   Mode = 2 or 3: Derivative values:                               */
 /*                                                                   */
      if ( MODE == 2 || MODE == 3 ) {
         for (i=0; i<nobs; i++ ) {
            JAC[dimx+i] = 2*X[dimx+i];
            }
         }
      }
/*
  Row 0 to nobs-1: The observation definitions:                       */
   else {
/*
   Mode = 1 or 3: Function value - x-part only                        */
 
      if ( MODE == 1 || MODE == 3 ) {
         k   = ROWNO*dimx;
         sum = 0.0;
         for ( j=0; j<dimx; j++ ) {
            sum = sum + A[k]*X[j] + B[k]*pow(X[j],2);
            k = k + 1;
            }
         *G  = sum;
      }
/*
   Mode = 2 or 3: Derivatives                                         */
 
      if ( MODE == 2 || MODE == 3 ) {
         k   = ROWNO*dimx;
         for ( j=0; j<dimx; j++ ) {
            JAC[j] = A[k] + 2*B[k]*X[j];
            k = k + 1;
            }
      }
   }
   return 0;
}
 

int COI_CALLCONV Leastsq_2DLagrStr( int HSRW[], int HSCL[], int* NODRV, int NUMVAR,
                                    int NUMCON, int NHESS, void* USRMEM )
{
   int i;
/*
  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.
                                                                      */
/*  Define structure of Hessian                                       */
 
   for (i=0; i<dimx+nobs; i++ ) {
      HSRW[i] = i;
      HSCL[i] = i;
   }
 
   return 0;
}
 

int COI_CALLCONV Leastsq_2DLagrVal( const double X[], const double U[], const int HSRW[], const int HSCL[], double HSVL[],
                                    int* NODRV, int NUMVAR, int NUMCON, int NHESS, void* USRMEM )
{
   int i;
/*
  See assumptions above
                                                                      */
/*  Evaluation mode                                                   */
 
   for (i=0; i<dimx; i++ ) HSVL[i] = 0.0;
   for (i=0; i<nobs; i++ ) HSVL[dimx+i] = 2.0*U[nobs];
 
   return 0;
}
 
#include "std.c"

int main(int argc, char** argv)
{
   c_log( "Starting to execute", START );
 
   defdata();
 /*
   Tell CONOPT about the sizes in the model
 */
   COI_Error +=  COIDEF_NumVar    ( CntVect, nobs+dimx    );       /* nobs+dimx variables                     */
   COI_Error +=  COIDEF_NumCon    ( CntVect, nobs+1       );       /* nobs+1 constraints                      */
   COI_Error +=  COIDEF_NumNz     ( CntVect, nobs*(dimx+2));       /* nobs*(dimx+2) nonzeros in the Jacobian  */
   COI_Error +=  COIDEF_NumNlNz   ( CntVect, nobs*(dimx+1));       /* nobs*(dimx+1) of which are nonlinear    */
   COI_Error +=  COIDEF_NumHess   ( CntVect, nobs+dimx );           /* nobs+dimx nonzeros in the Hessian       */
   COI_Error +=  COIDEF_OptDir    ( CntVect, -1           );       /* Minimize                                */
   COI_Error +=  COIDEF_ObjCon    ( CntVect, nobs         );       /* Objective is constraint nobs            */
   COI_Error +=  COIDEF_DebugFV   ( CntVect, 0            );       /* 0 means no debugging, 1 each iter       */
   COI_Error +=  COIDEF_StdOut    ( CntVect, 0            );       /* 1 means Allow output to StdOut          */
   COI_Error +=  COIDEF_Optfile   ( CntVect, "leastsq2.opt");      /* Register the options file               */
 /*
   Register the necessary callback routines with CONOPT
 */
   COI_Error +=  COIDEF_Message   ( CntVect, &Std_Message  );  /* Register the callback Message     */
   COI_Error +=  COIDEF_ErrMsg    ( CntVect, &Std_ErrMsg   );  /* Register the callback ErrMsg      */
   COI_Error +=  COIDEF_Status    ( CntVect, &Std_Status   );  /* Register the callback Status      */
   COI_Error +=  COIDEF_Solution  ( CntVect, &Std_Solution );  /* Register the callback Solution    */
   COI_Error +=  COIDEF_ReadMatrix( CntVect, &Leastsq_ReadMatrix);  /* Register the callback ReadMatrix  */
   COI_Error +=  COIDEF_FDEval    ( CntVect, &Leastsq_FDEval);      /* Register the callback FDEval      */
   COI_Error +=  COIDEF_2DLagrStr ( CntVect, &Leastsq_2DLagrStr);   /* Register the callback 2DLagrStr   */
   COI_Error +=  COIDEF_2DLagrVal ( CntVect, &Leastsq_2DLagrVal);   /* Register the callback 2DLagrVal   */
 
#if defined(LICENSE_INT_1) && defined(LICENSE_INT_2) && defined(LICENSE_INT_3) && defined(LICENSE_TEXT)
   COI_Error +=  COIDEF_License   ( CntVect, LICENSE_INT_1, LICENSE_INT_2, LICENSE_INT_3, LICENSE_TEXT);
#endif
 
   if ( COI_Error ) {
      printf("Skipping COI_Solve due to setup errors. COI_Error = %d\n",COI_Error);
      c_log( "Skipping Solve due to setup errors", COI_Error); }
   COI_Error =  COI_Solve        ( CntVect );           /* Optimize                      */
   printf("1 Thread: End of Least Squares Example 2. Return code =%d\n",COI_Error);
   if ( COI_Error ) {
      c_log( "1 Thread: Errors encountered during first solve", COI_Error); }
   else if ( stacalls == 0 || solcalls == 0 ) {
      c_log( "1 Thread: Status or Solution routine was not during first solve called", -1); }
   else if ( sstat != 1 || mstat != 2 ) {
      c_log( "1 Thread: Solver or Model status was not as expected (1,2)", -1); }
   else if ( fabs( OBJ-19.4443 ) > 0.001 ) {
      c_log( "1 Thread: Incorrect objective returned from first solve", -1); };
 
/* Register free use of threads and run again                                            */
 
   COI_Error +=  COIDEF_ThreadS  ( CntVect, 0 ); /* 0 means take what you can get        */
   COI_Error =  COI_Solve        ( CntVect );           /* Optimize                      */
   printf("Multi Thread: End of Least Squares Example 2. Return code =%d\n",COI_Error);
   if ( COI_Error ) {
      c_log( "Multi Thread: Errors encountered during first solve", COI_Error); }
   else if ( stacalls == 0 || solcalls == 0 ) {
      c_log( "Multi Thread: Status or Solution routine was not during first solve called", -1); }
   else if ( !( sstat == 1 && mstat == 2 ) ) {
      c_log( "Multi Thread: Solver or Model status was not as expected (1,2)", -1); }
   else if ( fabs( OBJ-19.4443 ) > 0.001 ) {
      c_log( "Multi Thread: Incorrect objective returned from first solve", -1); };
 
   c_log( "Successful Solve", OK );
}