opt_quadr.cc

COOOL：CWP面向对象最优化库（CWP Object Oriented Optimization Library） COOOL是C++类的一个集合

#include <iostream.h>
#include <stdlib.h>
#include <GetOpt.h>
#include <C++.hh>
#include <ObjFcns.hh>
#include <Optim.hh>

static int verbose        = 0;
static AString SlaveP     = "";
static AString SlavePArgs = "";
static int modelsize      = 3;
static int itMaxLine 	  = 10;
static int itMaxCG 	  = 20;
static double tol	  = 0.00001;

static char *opt_string = "S:A:D:N:T:V:M:L:U:C:";
extern "C" { int getpid(void); }

int cstraints(Vector<double>& v)
{
   for (int i=1; i<v.size(); i++)
      if (fabs(v[i]) > 5) return -1;
   
   return 1;
}

int main(int argc, char *argv[])
{
	int c; 		// used in data input 
	int i;		// counters
	char *  instring;
	Vector<double> delta(1), upper(1), lower(1);

        GetOpt getopt(argc, argv, opt_string);
        while (( c = getopt()) != EOF)
        switch(c) {
	  case 'S':
	    SlaveP = getopt.optarg;
	    break;
	  case 'A':
	    SlavePArgs = getopt.optarg;
	    break;
	  case 'N':
	    modelsize = atoi(getopt.optarg);
	    SlavePArgs += " -N ";
	    SlavePArgs += toAS(modelsize);
	    break;
	  case 'M':
	    itMaxLine = atoi(getopt.optarg);
	    break;
	  case 'C':
	    itMaxCG = atoi(getopt.optarg);
	    break;
	  case 'T':
	    tol = atof(getopt.optarg);
	    break;
	    case 'D':
	    instring = getopt.optarg;
	    delta[0] = atof(instring);
	    while((instring = strpbrk(instring, ",")) != NULL)
	    {
	       instring++;
	       delta.addVal(atof(instring));
	    }
	    break;
	    case 'L':
	    instring = getopt.optarg;
	    lower[0] = atof(instring);
	    while((instring = strpbrk(instring, ",")) != NULL)
	    {
	       instring++;
	       lower.addVal(atof(instring));
	    }
	    break;
	    case 'U':
	    instring = getopt.optarg;
	    upper[0] = atof(instring);
	    while((instring = strpbrk(instring, ",")) != NULL)
	    {
	       instring++;
	       upper.addVal(atof(instring));
	    }
	    case 'V':
	    i = atoi(getopt.optarg);
	    if (i) verbose++;
	    SlavePArgs += " -V ";
	    break;
	  default:
	    cerr << argv[0] << ": bad argument \n";
	    exit(1);
	}

	if( argc == 1) {
	    cerr << "Usage:  " << argv[0] << endl;
	    cerr << " -S name of the slave file" << endl;
	    cerr << " -N number of unknowns" << endl;
	    cerr << " -A arguments for the slave file" << endl;
	    cerr << " -D fraction of step size for numerical gradient" << endl;
	    cerr << " -M maximum number of iterations in line search" << endl;
	    cerr << 
	    " -C maximum number of iterations in conjugate gradient" << endl;
	    cerr << " -L lower boundary of model parameters" << endl;
	    cerr << " -U upper boundary of model parameters" << endl;
	    cerr << " -T tolerance of the residue" << endl;
	    cerr << " -V verbose" << endl;
	    exit(-1);
	}


	// construct the model space and the initial guess, 
	Vector<double> LOWER(modelsize, lower.size(), lower.toPointer());
	Vector<double> UPPER(modelsize, upper.size(), upper.toPointer());
	Vector<double> DELTA(modelsize, delta.size(), delta.toPointer());
	cerr << "Upper Boundary: \n" << UPPER <<endl;
	cerr << "Lower Boundary: \n" << LOWER <<endl;
	cerr << "Interval: \n" << DELTA << endl;

	//***NEW: TO SPECIFY POINTER TO THE CONSTRAINTS FUNCTION
	int (*tp)(Vector<double>&);
	tp = cstraints;

	Model<double> m(modelsize, tp);

	//construct objective function, connecting with user executable
        SlaveForward* my_slave = 
		new SlaveForward(modelsize, 1, SlaveP, SlavePArgs, verbose);
	DirectObjFcn*  fp = new DirectObjFcn(modelsize, my_slave, 1);

	//constrcut optimization objects
	CubicLineSearch* ls = new CubicLineSearch(fp,itMaxLine,&DELTA);

	// starting model is random 
	int seed=getpid();
	srand(seed);
	for (i=0; i<modelsize; i++) 
	   m[i] = (UPPER[i]-LOWER[i])*rand()/RAND_MAX+LOWER[i];
	cout << m << endl;

	// construct conjugate gradient object
	ConjugateGradient* opt = new ConjugateGradient(ls,itMaxCG,tol,verbose);
//	BFGS* opt = new BFGS(ls,itMaxCG,tol,1);

	double initialSp = fp->performance(m);
	cerr <<"Initial Value of the OBJ: " << initialSp << endl;

	m = opt->optimizer(m);
	double finalSp = fp->performance(m); 
	cerr <<"Final Value of the OBJ: " << finalSp << endl;
	
	cerr <<"The number of iterations: "<< opt->numIterations() << endl;
	cerr <<"Normalized residue: " << opt->normResidue() <<endl;
}