art.cc

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

//============================================================
// COOOL           version 1.1           ---     Nov,  1995
//   Center for Wave Phenomena, Colorado School of Mines
//============================================================
//
//   This code is part of a preliminary release of COOOL (CWP
// Object-Oriented Optimization Library) and associated class 
// libraries. 
//
// The COOOL library is a free software. You can do anything you want
// with it including make a fortune.  However, neither the authors,
// the Center for Wave Phenomena, nor anyone else you can think of
// makes any guarantees about anything in this package or any aspect
// of its functionality.
//
// Since you've got the source code you can also modify the
// library to suit your own purposes. We would appreciate it 
// if the headers that identify the authors are kept in the 
// source code.
//
//=============================
// author:  H. Lydia Deng, 08/08/95
// Implementation of the Algebraic Reconstruction Technique (ART)
//=============================

#include <ART.hh>

ART::~ART() 
{;}

ART::ART(int n, LinearForward* p, Vector<double>* data, 
			int it, double eps)
 : QuadraticOptima(n, p, data)
{
    	iterMax 	= 	it;
    	tol 		= 	eps;

	m		=	p->modelSize();
	n		=	p->dataSize();
}

ART::ART(int n, LinearForward* p, Vector<double>* data, 
			int it, double eps, int verb)
 : QuadraticOptima(n, p, data, verb)
{
    	iterMax 	= 	it;
    	tol 		= 	eps;

	m		=	p->modelSize();
	n		=	p->dataSize();
}

Model<double> ART::upDating(int irow, const Model<double>& model)
{
   Model<double> new_model(model);
   Vector<double> a(m), data(n);
   MisFitFcn* f = (MisFitFcn*) fp;
   
   data = f->getRefData();
   a = f->getOneOperator(irow);
   double scale = 1./a.norm2S();

   scale *= data[irow] - new_model.modParam() * a;
   new_model = new_model.update(1,scale,a);
   
   return new_model;
}


Model<long> ART::upDating(int i, const Model<long>& model)
{
   Model<double> new_model(model);
   Vector<double> a(m), data(n);
   MisFitFcn* f = (MisFitFcn*) fp;
   
   a = f->getOneOperator(i);
   data = f->getRefData();
   double scale = 1./a.norm2S();

   scale *= data[i] - new_model.modParam() * a;
   new_model = new_model.update(1,scale,a);
   
   return new_model;
}

Model<double> ART::optimizer(Model<double>& init_model)
{
   Model<double> model(init_model);
   double value = fp->performance(model);
   if (isVerbose) cerr << "the current residue: "<< value << endl;
   QuadraticOptima::appendResidue(value);

   int iter = 0;
   while(iter < iterMax)
   {
      int k = iter%m;
      model = upDating(k, model);
      iter ++;
      
      if (k == 0) 
      {
	 value = fp->performance(model);
	 if (isVerbose) cerr << "the current residue: "<< value << endl;
	 QuadraticOptima::appendResidue(value);
      }
      
    }
    return(model);
}

Model<long> ART::optimizer(Model<long>& model)
{
    Model<double> temp(model);
    temp = optimizer(temp);
    Model<long> optm(temp);
    return optm;
}