my_nr.h

Go to the documentation of this file.

/*
 * last update 29 dec 98
 * all matrix and vector index ranges are from 0..n-1 ! 
 */

/* Changed 98/12/29 by AMB */
/*         98/01/16 by KG  */
/* $Id: my_nr.h,v 1.10.2.14 2022/11/03 17:28:11 garloff Exp $ */

#ifndef TBCI_MY_NR_H
#define TBCI_MY_NR_H  

#include "tbci/vector.h"
#include "tbci/matrix.h"
#include "tbci/constants.h"
#include "tbci/solver/gauss_jordan.h"

#ifdef HAVE_SIGNALS
# include <signal.h>
#endif
/*
#include <cstdlib>
#include <cmath>
#include <iostream>
#include <iomanip>
*/


NAMESPACE_TBCI

#ifdef PRAGMA_I
# pragma interface
#endif

INST (template <typename T> class dummy friend void do_nothing(T);)
template <typename T>
inline void do_nothing(T dummy) {}

// the kronecker delta
#define DELTA(i,j) ((i==j)?1:0)

/*
inline int delta (const int i, const int j)
{
   return ((i==j) ? 1 : 0);
}
 */

#if defined(SIGHANDLE) && defined(HAVE_SIGNALS)
volatile int sigint = 0;

void breakhandler (int signum)
{
   STD__ cerr << "\n SIGINT !\n";
   sigint++;
   if (sigint > 2) 
     {
       STD__ cerr << "\n ... aborted due to too many ignored SIGINTs !\n";
       exit (1);
     }
   else signal (signum, breakhandler);
}
#endif

#define DELTAX 1e-30

INST(template <typename T> class Vector friend int basis_trafo (long int, long int, Vector<T>&);)
template <typename T>
int basis_trafo (long int value, long int basis, Vector<T>& index)
/* transforms the integer 'value' into the basis 'basis'; 
   index contain the resulting digits with lowest significant first */
{
   int n=index.size();
   int i=0;

   index=(T)0;
   if (value==0) return 1; //special zero treatment
   if (n<(int)ceil(MATH__ sqrt(2*log((double)value)/log((double)basis))) || basis<=1) return 0;
   while( value>0 )
   {
      index(i)=value%basis;
      value/=basis;
      i++;
   }
   return 1;
}


#define CON2 2.0
#define CON SQRT2
#define BIG 1.0e30
#define NTAB 10
#define SAFE 2.0

INST(template <typename T> class Vector friend T partial_del (T (*func)(const Vector<T>&, const Vector<T>&), const Vector<T>&, const Vector<T>&, int, T, T&);)
template <typename T>
T partial_del (T (*func)(const Vector<T>&, const Vector<T>&),
               const Vector<T>& x, const Vector<T>& p, 
               int mu, T h, T& err)
  // calculate partial derivative of func with respect to mu
  // Acc. to Num. recipes, ch. 5.7 
{
   //static T maxarg1 ALIGN(MIN_ALIGN), maxarg2 ALIGN(MIN_ALIGN);
   T errt ALIGN(MIN_ALIGN), fac, hh, ans(0);
   int xlen = x.size();
   Matrix<T> a(NTAB,NTAB);
   Vector<T> buf_plus(xlen);
   Vector<T> buf_minus(xlen);
#ifdef FOURTHORDER
   Vector<T> buf_plus2(xlen);
   Vector<T> buf_minus2(xlen);
#endif
   
   if (mu >= xlen)
     STD__ cerr << "\n error: partial_del: derivate variable index out of range \n";
   hh = h ;
   buf_plus = x; buf_minus = x; buf_plus(mu) += hh; buf_minus(mu) -= hh;
#ifdef FOURTHORDER
   buf_plus2 = x; buf_minus2 = x; buf_plus2(mu) += hh*2; buf_minus2 -= hh*2;
   a(0,0) = (-(*func)(buf_plus2, p) + 8.0*(*func)(buf_plus, p) 
             - 8.0*(*func)(buf_minus, p) + (*func)(buf_minus2, p)) / (hh*12.0);
#else
   a(0,0) = ((*func)(buf_plus, p) - (*func)(buf_minus, p)) / (2.0*hh);
#endif
   err = BIG;
   for (int i = 1; i < NTAB; i++)
     {
        hh /= CON;
        buf_plus = x; buf_minus = x; buf_plus(mu) += hh; buf_minus(mu) -= hh;
#ifdef FOURTHORDER
        buf_plus2 = x; buf_minus2 = x; buf_plus2(mu) += hh*2; buf_minus2 -= hh*2;
        a(0,i) = (-(*func)(buf_plus2, p) + 8.0*(*func)(buf_plus, p) 
                  - 8.0*(*func)(buf_minus, p) + (*func)(buf_minus2, p)) / (hh*12.0);
#else
        a(0,i) = ((*func)(buf_plus, p) - (*func)(buf_minus, p)) / (2.0*hh);
#endif
        fac = CON2;
        for (REGISTER int j = 1; j <= i; j++)
          {
             a(j,i) = (a(j-1,i)*fac - a(j-1,i-1))/(fac-1.0);
             fac *= CON2;
             errt = MAX (MATH__ fabs(a(j,i)-a(j-1,i)), MATH__ fabs(a(j,i)-a(j-1,i-1)));
             if (errt <= err) {
                err = errt;
                ans = a(j,i);
             }
          }
        if (MATH__ fabs (a(i,i)-a(i-1,i-1)) >= MATH__ fabs (SAFE*err)) break;
     }
   return ans;
}

INST(template <typename T> class Vector friend T chisq (Vector<T>&, Vector<T>&, Vector<T>&);)
template <typename T>
T chisq (Vector<T>& data, Vector<T>& func, Vector<T>& sig)
   // calculates the chi-squared sum

{
   T chi ALIGN(MIN_ALIGN) = 0;
   int len=data.size();

   for(REGISTER int i=0;i<len;i++)  chi+=sqr((data(i)-func(i))/sig(i));
   return(chi);
}

INST(template <typename T> class Vector friend T chisquare (T (*func)(const Vector<T>&, const Vector<T>&), Vector<T>&, Vector<T>&, Vector<T>&, Vector<T>&);)
template <typename T>
T chisquare (T (*func)(const Vector<T>&, const Vector<T>&), 
            Vector<T>& param, Vector<T>& x, 
            Vector<T>& y, Vector<T>& sigma) 
  // calculates function values of func at position x with params param
  // and return chisq with given y values weighted with sigmas
{
  Vector<T> arg(1);
  
  T chi ALIGN(MIN_ALIGN) = 0, wfvaldiff;
  int anzpkt = x.size();
  
  for (REGISTER int i = 0; i < anzpkt; i++)
    {
      arg(0) = x(i);
      wfvaldiff = (y(i) - (*func)(param, arg)) / sigma(i);  
      chi += sqr (wfvaldiff);
    }
  return chi;
}

#if 0
//class nslist<minimum<T>,<T>> minima;

template <typename T>
class minimum
{
 public:
   Vector<T> koord;
   T value ALIGN(MIN_ALIGN);
   T getval () { return value; }
   minimum () { value = 1e38; }
};

INST(template <typename T> class Vector friend TVector<T> grid_mins(T (*func)(const Vector<T>&, const Vector<T>&), Vector<T>&, Vector<T>&, Vector<T>&, Vector<T>&, Vector<T>&, int, char, int, nsList<minimum<T>,T>& );)
template <typename T>
TVector<T> grid_mins(T (*func)(const Vector<T>&, const Vector<T>&),
                 Vector<T>& x, Vector<T>& y,
                 Vector<T>& sigma, Vector<T>& pmin, Vector<T>& pmax,
                 int res, char v, int anzmin,
                 nsList<minimum<T>,T>& minima)
                 //return absmin;
  // find minima of func on grid with resolution res
{
  int anzpar = pmin.size();
  minimum<T>* min; 
  minimum<T> min2; min2.koord.resize (anzpar);
  Vector<T> actpar(anzpar);
  TVector<T> absmin(anzpar);
  //absmin.resize(anzpar);
  Vector<T> factor(anzpar);
  Vector<T> plim(anzpar);
  long int len = (long int) MATH__ pow ((double)res, anzpar);
  bool overflow; int k;
  
  actpar = pmin;
  T actval ALIGN(MIN_ALIGN) = chisquare (func, actpar, x, y, sigma);
  min = minima.setlast();
  STD__ cerr << "Parameter space discretisation in process ...\n";
  //factors
  for (REGISTER int j = 0; j < anzpar; j++)
    {
      factor(j) = (pmax(j) - pmin(j)) / (res - 1);
      plim(j) = pmax(j) + 0.5 * factor(j);
    }
  //discretisation
  for (long int i = 0; i < len; i++)
    {
      actval = chisquare (func, actpar, x, y, sigma);
      if (actval < min->value)
        {
           if (v=='y' || v=='Y') 
                STD__ cerr << "Found new min " << actval << " at " << actpar << '\n';
#ifdef USEQSORT
           min->koord = actpar;
           min->value = actval;
           if (anzmin > 1) minima.qsort ();
           min = minima.setlast ();
#else
           min2.koord = actpar; min2.value = actval;
           min = minima.sortin_r ( min2 );
#endif
        }
      if (!(i % 1000)) STD__ cerr << i << '\r';
      //avoid basis_trafo
      overflow = true; k = 0;
      while (overflow && (k < anzpar))
        {
          actpar(k) += factor(k);
          if (actpar(k) > plim(k)) actpar(k) = pmin(k);
          else overflow = false;
          k++;
        }
#ifdef SIGHANDLE
       if (sigint) { sigint--; return 0; }
#endif
    }
  //STD__ cerr << "\nLast par: " << actpar << "\n";
  STD__ cerr << "\n ... done.\n\n";
  absmin = minima.getfirst()->koord;
  return absmin;
}
#endif /* 0 */

INST(template <typename T> class Vector friend TVector<T> grid_min(T (*func)(const Vector<T>&, const Vector<T>&), Vector<T>&, Vector<T>&, Vector<T>&, Vector<T>&, Vector<T>&, Vector<T>&, long int, int, char);)
template <typename T>
TVector<T> grid_min(T (*func)(const Vector<T>&, const Vector<T>&), 
                    Vector<T>& x, Vector<T>& y,
                    Vector<T>& sig, Vector<T>& pmin, 
                    Vector<T>& pmax, Vector<T>& lambda,
                    long int res, int max_it, char v)
   // find minimun of func on grid with resolution res
{

   int plen=pmin.size();
   int xlen=x.size();
   long int len=0;
   Vector<T> p(pmin);
   TVector<T> min_pos(pmin);
   Vector<T> factor(plen);
   Vector<T> arg(1);
   Vector<T> fval(x.size());
   Vector<T> index(plen);
   Vector<T> b(plen,0);
   Vector<T> a(plen,0);
   T min ALIGN(MIN_ALIGN) = 0, chi=0, err=0.0, fac=1e-8;

   len=(long int) pow((double)res, plen);
   for (int s=0; s<xlen; s++) {  arg(0)=x(s); fval(s)=(*func)(p, arg);  }
   min=chisq(y, fval, sig);
   STD__ cerr << "\n parameter space discretisation in progress ... \n ";
   for (REGISTER int j = 0; j < plen; j++)
     factor(j) = (pmax(j) - pmin(j)) / ((T)res - (T)1);
   for (long int i=0;i<len;i++)
   {
      // linearize recursive for structure by basis trafo
      if ( !basis_trafo(i, res, index) )
         STD__ cerr << "\n error: basis trafo failed\n";
      for (REGISTER int j=0;j<plen;j++)
         p(j) = factor(j) * index(j) + pmin(j); //p-space discr.
      //do max_it steepest gradient steps
      for (int l=0;l<max_it;l++)
      {
         // setup Vector
         err=0.0;
         for (int k=0;k<plen;k++)
         {
            for (int j=0;j<xlen;j++)
            {
               arg(0)=x(j);
               b(k)+=(y(j)-(*func)(p, arg))/(sig(j)*sig(j))*
                  partial_del(func, p, arg, k, x(j)*fac, err);
            }
         }
         // setup diagonal Matrix for steepest gradient const
         err = 0.0;
         for (int t=0; t<plen; t++)
         {
            for (REGISTER int u=0; u<xlen; u++)
            {
               arg(0) = x(u);
               a(t) += sqr(partial_del(func, p, arg, t, x(u)*fac, err))/
                  (sig(u)*sig(u));
            }
         }
         for (int q=0; q<plen; q++)
           p(q) += b(q)*lambda(q)/(a(q)+DELTAX);
      }
      // calc func values with new parameter
      for (int q=0; q<xlen; q++) {  arg(0)=x(q); fval(q)=(*func)(p, arg);  }
      // calc chi2 sqared
      chi=chisq(y, fval, sig);
      if ( chi<min ) {  min=chi; min_pos=p;  } //choose new min
      if ( i%1000==0 ) STD__ cerr << i << "\r";
#ifdef SIGHANDLE
      if (sigint) { sigint--; return 0; }
#endif
   }
   STD__ cerr << " \n done! \n\n\n";
   return min_pos;
}

INST(template <typename T> class Vector friend T lev_mar(Vector<T>&, Vector<T>&, Vector<T>&, Vector<T>&, Vector<T>&, Vector<T>&, T (*func)(const Vector<T>&, const Vector<T>&), Vector<T>, Vector<T>&, T, int, char);)
template <typename T>
T lev_mar(Vector<T>& x, Vector<T>& y, Vector<T>& sig, Vector<T>& param,
          Vector<T>& plow, Vector<T>& phigh,
          T (*func)(const Vector<T>&, const Vector<T>&),
          Vector<T> lambda, Vector<T>& lscale,
          T tol, int iter, char ver)
/* fit method: levenberg-marquardt
   leads to an linear equ
   derivatives are calculated numerically by num_derivatives 
   parameter space is restricted by plow and phigh */

{
   T chi ALIGN(MIN_ALIGN), chi_new, chi_min, chi_dif, err=0.0, fac=1.0e-8;
   int len=x.size();
   int plen=param.size();
   Vector<T> param2(plen); // Vector<T> for new parameter set
   Vector<T> param_best(plen); // Vector<T> for best parameter set
   Vector<T> fval(len); // Vector for fit-func values
   Vector<T> arg(1); // Vector for arguments of fit-func (now only 1d)

   // comp chi-squared
   chi=chisquare(func,param,x,y,sig);
   STD__ cerr << " \n chi-squared with start parameters is: " << chi;
   chi_min=chi;
   // ini matrices
   Matrix<T> a(plen,plen);
   Matrix<T> b(plen,1);
   STD__ cerr << "\n initializing matrices and starting loop... \n";
   int j=0; // iteration counter
   param_best=param;
   do
   {
      // setup Matrix
      for (int k=0;k<plen;k++)
         for (int l=0;l<plen;l++)
         {
            for (int i=0;i<len;i++)
            {
               arg(0)=x(i);
               a(k,l)+=partial_del(func, param, (Vector<T>&)arg, k, x(i)*fac, err)*
                 partial_del(func, param, (Vector<T>&)arg, l, x(i)*fac, err) /(sig(i)*sig(i));
            }
            a(k,l)*=((T)1+lambda(k)*(T)DELTA(k,l));
         }
      // setup Vector<T>
      for (int q=0; q<plen; q++)
      {
         for (int i=0;i<len;i++)
         {
            arg(0)=x(i);
            b(q,0)+=(y(i)-(*func)(param, (Vector<T>&)arg))/(sig(i)*sig(i))*
              partial_del(func, param, (Vector<T>&)arg, q, x(i)*fac, err);
         }
      }
      // solve equ a.x=b
      char exceed = gaussj (a, b);
      // calc new parameter
      for (REGISTER int t=0; t<plen; t++) param2(t)=param(t)+b(t,0);
      // check parameter range
      exceed = 0;
      for (REGISTER int i=0; i<plen; i++)
      {
         if ((param2(i) > phigh(i)) && (phigh(i)!=plow(i)) ) exceed = 1;
         if ((param2(i) < plow(i))  && (phigh(i)!=plow(i)) ) exceed = 1;
      }
      // calc chi2 sqared
      chi_new = chisquare(func, param2, x, y, sig);
      chi_dif = MATH__ fabs (chi_new - chi);
      if ((chi_new >= chi) || (exceed == 1)) lambda *= lscale(0);
      //increase lambda cause downhill is exhausted or limits are exceeded
      else {
              // decrease lambda and update trial solution
              lambda*=(T)1/lscale(1);
              param=param2;
              // update best est. param if new min is reached
              if (chi_new<chi_min) {
                                      chi_min=chi_new;
                                      param_best=param;
                                   }
           }
      chi=chi_new; // update chi
      if (ver=='y' || ver=='Y')
         STD__ cerr << "\n iteration " << j << "   " << param <<
         "   " << param2 << "   " << chi_new;
#ifdef SIGHANDLE
      if (sigint) { sigint--; param = param_best; return 1e38; }
#endif
      j++;
      if (j%10 == 0)  STD__ cerr << ".";
      // end loop if tolerance is reached or number of iterations exceeded
   } while( (chi_dif >= tol) && (j < iter) );
   param=param_best; // return best parameter set found
   STD__ cerr << "\n loop finished after " << j << " iterations with last delta chi: "
      << chi_dif << "\n and chi-squared: " << chi_new << "\n";
   return (chi_new);
}

INST(template <typename T> class Matrix friend void hpsort (Matrix<T>&, int);)
template <typename T>
void hpsort (Matrix<T>& x, int pos)
/* sorts rows of input Matrix x into ascending numerical order of values 
   in colomn pos;
   uses heap sort algorithm, rf. nr p 336 */

{
   unsigned long i, ir, j, l, n;
   Vector<T> rra(x.columns());

   n=x.rows();
   if ( n<2 ) return;
   l=(n >> 1)+1;
   ir=n;
   for (;;)
   {
      if ( l>1 ) {
                    rra=x(--l-1);
                 }
      else {
              rra=x(ir-1);
              x(ir-1)=x(0);
              if (--ir==1) {
                              x(0)=rra;
                              break;
                           }
           }
      i=l;
      j=l+l;
      while ( j<=ir ) {
                         if ( j<ir && x(j-1,pos) < x(j,pos) ) j++;
                         if ( rra(pos)<x(j-1,pos) ) {
                                                       x(i-1)=x(j-1);
                                                       i=j;
                                                       j <<= 1;
                                                    }
                         else j=ir+1;
                      }
      x(i-1)=rra;
   }
}

NAMESPACE_END

#endif /* TBCI_MY_NR_H */