leastsq10.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 C[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_2DDirIni( const double X[], const double DX[], const int ROWLIST[], int LISTSIZE,
                               int NUMTHREAD, int NEWPT, int* NODRV, int NUMVAR, void* USRMEM )
{
int i, j, k;
 
/*  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         */
 
      k   = 0;
      for ( i=0; i<nobs; i++ ) {
         for ( j=0; j<dimx; j++ ) {
            C[k] = 2.0*B[k]*DX[j];
            k++;
         }
      }
   return 0;
}
 

int COI_CALLCONV Leastsq_2DDir( const double X[], const double DX[], double D2G[],
                            int ROWNO, const int JACNUM[], int* NODRV, int NUMVAR, int NUMJAC,
                            int THREAD, void* USRMEM )
{
int i, j, k;
 
   if ( ROWNO == nobs ) {
 
/*  Objective row: sum(i, res(i)**2 )                                 */
 
      for ( j=0; j<dimx; j++ ) D2G[j] = 0;
      for ( i=0; i<nobs; i++ ) D2G[dimx+i] = 2.0*DX[dimx+i];
      }
   else {
 
/*  Constraint row with b(i,j)*x(j)**2 as only nonlinear term         */
 
      k   = ROWNO*dimx;
      for ( j=0; j<dimx; j++ ) {
         D2G[j] = C[k];
         k++; }
      for ( i=0; i<nobs; i++ ) D2G[dimx+i] = 0;
   }
   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_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, -1 first */
   COI_Error +=  COIDEF_Debug2D   ( CntVect, 0             );       /* 0 means no debugging, 1 each iter, -1 first */
   COI_Error +=  COIDEF_StdOut    ( CntVect, 0             );       /* 1 means Allow output to StdOut          */
   COI_Error +=  COIDEF_Optfile   ( CntVect, "leastsq10.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_2DDir     ( CntVect, &Leastsq_2DDir );      /* Register the callback 2DDir       */
   COI_Error +=  COIDEF_2DDirIni  ( CntVect, &Leastsq_2DDirIni );   /* Register the callback 2DDirIni    */
 
#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("End of Least Squares Example 10. Return code =%d\n",COI_Error);
   if ( COI_Error ) {
      c_log( "Errors encountered during solve", COI_Error); }
   else if ( stacalls == 0 || solcalls == 0 ) {
      c_log( "Status or Solution routine was not called", -1); }
   else if ( sstat != 1 || mstat != 2 ) {
      c_log( "Solver or Model status was not as expected (1,2)", -1); }
   else if ( fabs( OBJ-19.4443 ) > 0.001 ) {
      c_log( "Incorrect objective returned from first solve", -1); };
 
   c_log( "Successful Solve", OK );
}