sirt.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/14/95
// Implementation of the Simulaneous Iterative Reconstruction Technique (SIRT)
//==========================================================================

#include <SIRT.hh>

SIRT::~SIRT()
{;}

SIRT::SIRT(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();
}

SIRT::SIRT(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> SIRT::upDating(int irow, const Model<double>& model)
{
   Model<double> new_model(model);
   Vector<double> a(m), data(n), b(m);
   MisFitFcn* f = (MisFitFcn*) fp;
   double scale, sigma=0;
   int i;
   
   data = f->getRefData();
   b = 0.;
   for (i=0; i<m; i++)
   {
      a = f->getOneOperator(i);
      scale = data[i] - model.modParam()*a;
      b += scale*a;
      sigma += a.norm2S();
   }
   sigma = 1./sigma;

   new_model = model.update(1,sigma,b);
   
   return new_model;
}


Model<long> SIRT::upDating(int irow, const Model<long>& model)
{
    Model<double> temp(model);
    temp = upDating(irow, temp);
    Model<long> new_model(temp);

    return new_model;
}

Model<double> SIRT::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> SIRT::optimizer(Model<long>& model)
{
    Model<double> temp(model);
    temp = optimizer(temp);
    Model<long> optm(temp);
    return optm;
}