irls.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, 03/06/94, 03/14/94
// Definition of the IRLS
//=============================

#include <Irls.hh>

void  IterativeReweightedLS::updateWeight(Vector<double>& r)
{
   int i, n = r.size();
   
   DiagMatrix<double> d(n);
   for (i=0; i<n; i++) 
      d[i] = pow(fabs(r[i]),(2-order));
   
   double temp = d.min();
   double eps = temp + taper*(d.max()-temp);
	
   for(i=0; i<n; i++) 
   {
      temp = d[i];
      if(temp < eps) temp = eps;
      d[i] = eps/temp;
   }
	
   weight[0] = d;
}

IterativeReweightedLS::IterativeReweightedLS(int n, LinearForward* p, Vector<double>* data, 
	int nOrder, int w_nit, int nit, double tape, double eps)
: QuadraticOptima(n, p, data)
{
    order		=	nOrder;
    iterMax 		= 	nit;
    if (order == 2) 
	niterWeight 	= 	1;
    else
	niterWeight	=	w_nit;
    taper		=	tape;
    pcgls		=	new PreconditionedCGLS(n, p, data, iterMax, eps);
    weight 		= 	new DiagMatrix<double>(p->dataSize());
    weight[0]		=	1.;
}

IterativeReweightedLS::IterativeReweightedLS(int n, LinearForward* p, Vector<double>* data, 
	int nOrder, int w_nit, int nit, double tape, double eps, int verb)
: QuadraticOptima(n, p, data, verb)
{
    order		=	nOrder;
    iterMax 		= 	nit;
    if (order == 2) 
       niterWeight 	= 	1;
    else
       niterWeight	=	w_nit;
    taper		=	tape;
    pcgls		=	new PreconditionedCGLS(n, p, data, iterMax, eps, verb);
    weight 		= 	new DiagMatrix<double>(p->dataSize());
    weight[0]		=	1.;
    numIter             =       0;
}

Model<double> IterativeReweightedLS::optimizer(Model<double>& model)
{
   Model<double> optm(model);
   
   pcgls->assignWeight(weight[0]);
   optm = pcgls->optimizer(model);
   Vector<double> error = pcgls->currentError();
   this->updateWeight(error);
   residue[0] = pcgls->allResidue();
   numIter = pcgls->numIterations();
   
   int iloop = 1;
   Model<double> temp(optm);
   while (iloop < niterWeight && residue->last() > tol)
   {
      pcgls->assignWeight(weight[0]);
      optm = pcgls->optimizer(temp);
      Vector<double> error = pcgls->currentError();
      this->updateWeight(error);
      residue[0] = pcgls->allResidue();
      numIter = pcgls->numIterations();
      temp = optm;

      iloop++;
   }

    return optm;
}


Model<long> IterativeReweightedLS::optimizer(Model<long>& model)
{
   Model<double> mtemp(model);
   mtemp = optimizer(mtemp);
   Model<long> mopt(mtemp);
   
   return mopt;
}