marqfitp.h

二次最优化的C/C++源代码

// marqfitp.h
// written at Masaryk University, Brno, Czech Republic
// Freeware version. It is distributed without any warranty.
/**************************************************************
     Puvodni pascalskou unitu MarqFitP.pas prepsal do C++
     v kvetnu 1994 Petr Mikulik
     [e-mail], http://www.sci.muni.cz/~mikulik/

Toto je MarqFitP verze duben 1998.

Poznamky k verzi C++
* Petr Mikulik, kveten 1994
  - zadny z autoru neprotestuje, aby se tento program mohl sirit
    jako public domain, a to za predpokladu, ze jmena autoru
    zustanou uvedena. Pripadne zmeny v teto unite, ci chyby,
    hlaste autovi. Vzhledem k tomu, ze fitovaci rutina je objekt,
    muzes si napsat vlastni dedice s novymi vlastnostmi a tuto
    unitu pouzivat jako zdroj predchudcu.
  - do baliku techto unit (marqfitp.h + marqfitp.cpp)
    jsem pridal i program TESTMQ.CPP na otestovani tehle unity.
    (puvodni pascalska unita je jediny soubor MARQFITP.PAS,
    takze ten priklad je prilozen na jeji konec.)
  - vuci puvodni pascalske unite (viz dale) jsem do C++ verze
    pridal i moznost vypisovani chyb fitovanych parametru
    i behem vlastniho fitovani
  - na rozdil od pascalu se zde vyuziva specifickeho chovani
    Cecka vuci polim, ktera se deklaruji jako  parfloat *,
    datafloat * apod.
  - do budoucna: mozna ze pridelam i definici hranic,
    v nichz se mohou fitovane parametry menit.
* Radek Pavelka, leden 1995
  - hranice fitovanych parametru lze zadat tak, ze se krome poli
    parametru B vytvori dalsi dve pole - Bmin a Bmax. Ve volani Marq()
    se pak pridaji do seznamu argumentu za B. Pro vysledne parametry
    pak jsou Bmin[i] <= B[i] <= Bmax[i]. Ohraniceni lze vypnout
    pridanim parametru 0 ve volani Marq();
  - doplneni vraceni uplynuleho casu (v sekundach)
    GetElapsedTime() vraci dobu od predchoziho volani funkce,
    GetElapsedTime(1) vraci celkovy uplynuly cas.
* Petr Mikulik, leden 1995
  - preklad i pod ANSI C++
  - program PMFIT.CPP pro fitovani experimentalnich dat
    ctenych ze souboru pro lib. funkci: ulozen na ftp server
    a distrubuovan jako pmfit.lzh ci pmfit.zip.  Je podstatne
    rozsahlejsi nez demo-program TESTMQ.CPP, protoze se nejedna
    o demo, ale o plne funkcni program
* Petr Mikulik, leden 1997
  - par nevyznamnych zmen k lepsimu, aby to slo zkompilovat i Watcom Ceckem
* Petr Mikulik, duben 1998
  - opravy chyby v rutine Vystup, ktera chybne vypisovala pouze odmocniny 
    chyb (odmocninu z GetSigma)

Komentar k puvodni pascalske unite MARQFITP.PAS:

  Nelinearni regrese Marquartovou-Levenbergovou metodou
* proceduru Marq a pomocne napsal Jan Cely,
  katedra teoreticke fyziky, Prir. fak. MU, 611 37 Brno
  verse 2.0 / 01-III-90

Unitu dale upravovali (katedra fyziky pevne faze):
* Frantisek Vizda:
  - pole odchylek Vaha pro zapocet vah jednotlivych dat
    (mereni) pri vypoctu residua
  - dodana procedura na inverzi symetricke matice - pro hessian
* Petr Mikulik   (e-mail: [e-mail]):
  - zabudovani polozky Index pri volani modelove funkce
    (zjistuji se u ni jen fcni hodnoty v exp. bodech
  - prepsani puvodne procedur do objektu a patricne modifikace
    pouzitych procedur
  - zrychleni fitovani eliminovanim opakujicich se vypoctu se
    stejnym vektorem parametru B, zejmena u derivaci
  - statisticke funkce
  - doporuceni overene praxi: je vhodne fitovani provadet pouze
    v aritmetice  float  (rychlost, chyby ve vstupnich datech
    jsou vetsi nez presnost teoretickeho modelu s danymi vypocetnimi
    chybami). Tohle pravidlo fungovalo v pascalu, v C++ se stejne
    skoro vsechno pocita v double, takze...
  - jestlize navod k pouzivani MarqFitP neni zcela zrejmy,
    tak si ten testovaci programek vytahni do zvlastniho
    souboru, spust a podivej se, jak funguje
**************************************************************/


#ifndef MARQFIT_H
#define MARQFIT_H

#define ElaTime               // vraceni uplynuleho casu

#ifdef ElaTime
#include <time.h>
#endif

#include <fstream.h>


typedef float fitfloat;
typedef float parfloat;
typedef float datafloat;
typedef float longvecfloat;

/*
typedef double fitfloat;
typedef double parfloat;
typedef double datafloat;
typedef double longvecfloat;
*/

/* Jestli si myslis, ze tvoje spatne namerena data ma smysl
   fitovat scestnym teoretickym modelem v dvojnasobne
   presnosti, tak si sem dej double
Ovsem na druhou stranu, veskere vypocty v C jedou v aritmetice double,
takze... je treba vzdy otestovat, jakou aritmetiku ma smysl pouzivat
*/

extern const float FitMachEps;     // = 1.2e-7  FLT_EPSILON
extern const float FitSqrtMachEps; // = 3.5e-4

/*
extern const double FitMachEps;     // = 2.3e-16  FLT_EPSILON
extern const double FitSqrtMachEps; // = 1.52e-8
*/

// GLOBALNI POZNAMKA: vsechna data se berou od 1..n_dat, 1..n_par !!!
// neboli prvek 0 se nikde nepouziva ( x[0] of x[0..n+1] == x[-1] of x[1..n] )

// deklarace v C jsou dynamicke, pres pointery

class TMarq
{
  public:
    int citer;
    datafloat *D0;  // values calc. previously
  protected:
    int n_dat,n_par;
    int cresid,cderiv,cprint,cfitfun;
    int sizeofA, sizeofB;
    int ShowErrors, InvYes;
    datafloat *R;
    longvecfloat *SaveA;
    datafloat SumR1;   // na konci se sem da res. suma ctvercu
    int Eval1;     // zaznamena vysledek procedury Marq
  public: // 3.8.96
    int StopMarq; //Pro StopMarq=TRUE se ukonci iterace
  protected:
    int _Indx_Mr; //index prave pocitaneho rezidua
    datafloat Lambda;
    #ifdef ElaTime
    clock_t beginTime,endTime,lastTime;
    char elTime[15];
    #endif
    char TraceFile [80];
  public:
    datafloat SumY2; // sum of squares of Y data values
    TMarq ( const char *ATraceFile, int AShowErrors ); // constructor
    virtual ~TMarq ();                // destructor
    virtual datafloat FitFun ( datafloat X, parfloat *B,
		      int& Fres );
       //  Modelova funkce. Musi byt prepsana !
    virtual void DoJac ( int I, datafloat *X, parfloat *B,
	       parfloat *D, int& dres );
       // Procedura pro vypocet derivaci - jacobiho matice.
       // Muze se prepsat, pokud znas analytickou a nepouzivas numerickou derivaci
    virtual void FormMat ( int& Suc, longvecfloat *A,
	datafloat *X, datafloat *Vaha,
	parfloat *B, parfloat *V );
    virtual void Residue ( datafloat& Res, int& Suc,
      datafloat *X, datafloat *Y, datafloat *Vaha, parfloat *B );
    void Marq
      ( int _n_dat, datafloat *X, datafloat *Y, datafloat *Vaha,
	int _n_par, parfloat *B, datafloat& SumR, int& Eval,
	int Print );
    void Marq
      ( int _n_dat, datafloat *X, datafloat *Y, datafloat *Vaha,
	int _n_par, parfloat *B, parfloat *Bmin, parfloat *Bmax,
        datafloat& SumR, int& Eval, int Print, int b=1);
    virtual void Vystup ( datafloat SumR, parfloat *B, int Eval );
    datafloat GetDisperze();
    datafloat GetSigma ( int Par );
    datafloat GetKorelKoef ( int Par1, int Par2 );
    void PomVystup ( ostream& Out, datafloat SumR, parfloat *B );
    void PomInfoVystup ( ostream& Out );
    void InfoVystup ( ostream& Out );
    #ifdef ElaTime
    double GetElapsedTime(int i=0);
    #endif
};

class TMarqData: virtual public TMarq
{
  public:
    TMarqData ( const char *ATraceFile, int AShowErrors ):
      TMarq(ATraceFile,AShowErrors) {};
    virtual void FitFunData ( datafloat *X, parfloat *B,
      datafloat *D, int& Fres );
      // Modelova funkce. Musi byt prepsana !
    void DoJacData ( int Par, datafloat *X,
      parfloat *B, datafloat *Der, int& dres );
    virtual void Residue ( datafloat& Res, int& Suc,
      datafloat *X, datafloat *Y, datafloat *Vaha, parfloat *B );
    virtual void FormMat ( int& Suc, longvecfloat *A,
      datafloat *X, datafloat *Vaha,
      parfloat *B, parfloat *V );
};

/*-------------------------------------------------------------
Nasledujici help byl prevzat z pascalske unity, a to beze zmeny.
Do Cecka si to prelouskej sam.
Kdyztak si spust demo na fitovani.

POSTUP pri fitovani:
1. utvorime naslednika objektu s nami pozadovanou ucelovou
   funkci FitFun
2. umime-li spocitat analyticky derivace, prepiseme take DoJac
3. deklarace: var MujMarq = ^TMujMarq;
4. v programu: new(MujMarq,Init('trace.fit'));
   spusteni fitovani: MujMarq^.Marq ( _parametry_ );
   dale napr. writeln(MujMarq^.GetDisperze);
   konec: dispose(MujMarq,Done);
- nejlepe viz priklad, ktery je v baliku techto fitovacich programu

Soubor TraceFile: je-li to nenulovy retezec, pak se bude pri
kazdem vystupu postupu fitovani na obrazovku procedurou Vystup
provadet zapis vysledku i do tohoto souboru.

VSTUPni parametry v Marq:
 _n_dat = pocet vstupnich dvojic (x[i],y[i]),
      x = pole nezavisle promennych x[1],...,x[n_dat],
      y = pole zavisle promennych y[1],...,y[n_par],
    Vaha = pole vah v bodech x[i]
 FitFun = ucelova funkce (i),
  DoJac = jmeno procedury (ii), pocitajici radek Jacobiho matice
	  (standardne se pocita numericka derivace)
 _n_par = pocet hledanych parametru B[i] v modelu Fmarq,
      B = vektor, ktery ma v B[1],...,B[n_par] pocatecni odhad
	  hledanych parametru,
  print = urcuje vystup behem minimalizace:
	  print<=0 = bez vystupu, print>0 = krok vystupu.
  Vystup= procedura pro sledovani prubehu minimalizace
	  Je-li parametr print>0,vola se procedura po kazdem
	  iteracnim kroku delitelnem cislem "print".

VYSTUPni parametry v Marq:
    B = vektor, ktery ma v B[1],...,B[n_par] hodnoty parametru
	B pro nalezene minimum rezid. souctu ctvercu,
 sumr = rezidualni soucet ctvercu pro nalezene B,
 eval = [TRUE/FALSE] jestlize Marq ukoncil praci
	[normalne/nemohl vypocitat funkci nebo derivaci]


Fitujici procedura Marq v objektu vyzaduje deklaraci:

(a) Procedura pro MATEMATICKY MODEL (ucelova funkce)
    y = f(x,B)  funkci:
    function FitFun ( X: float; Index: integer; B: ParVec;
		      var Fres: boolean): float; virtual;
	  x = nezavisle promenna ( x==datax[index] )
	  B = vektor hledanych parametru
       fres = [TRUE/FALSE] byl-li vypocet funkce [mozny/nemozny]
    Tato funkce musi byt vzdy prepsana !

(b) Procedura pro VYPOCET PARCIALNICH DERIVACI rezidua r[i]
    podle parametru B[j]
    (tj. i-teho radku Jacobiho matice dr[i]/db[j],
     j=1,...,n_par a i=1,...,n_dat)
    procedure DoJac ( I: integer; var X: Data; var B: ParVec;
	       var D: ParVec; var dres: boolean ); virtual;
    kde navic k (a) je
	i = index rezidua (radku Jacobiho matice), i=1,...,n_dat,
    datax = vektor nezavisle promennych x[1],...,x[n_dat],
	d = vystupni vektor (i-ty radek) vypoctenych derivaci,
     dres = [TRUE/FALSE] byl-li vypocet [mozny/nemozny].

    Implicitne je naprogramovana procedura pro numerickou
    derivaci. Derivace se nahrazuji pravymi diferencemi s krokem
     h_j = (Abs(B[j] + _sqrteps)*_sqrteps
    Pokud pouzivas analytickou derivaci, muzes tuto proceduru
    prepsat.

---------------------------------------------------------------
POZNAMKA:
 Procedury Marq konci hledani minima az se zadny z hledanych
 parametru jiz nemeni v ramci presnosti stroje;
 Praci procedur lze z klavesnice prerusit stiskem <DEL>
 nebo nastavenim promenne StopMarq na hodnotu TRUE


  Statisticke funkce: smeji byt pouzity v pripade uspesne
dokonceneho fitovani procedurou Marq (v pripade neuspechu
budes mit chyby 1e38 nebo 1e19):
function GetDisperze: vyhodi disperzi
function GetSigma(Par): vyhodi stredni kvadratickou odchylku
  pro dany parametr   { sigma^2 = disperze veliciny }
function GetKorelKoef(Par1,Par2): vyhodi korelacni koeficient
  pro danou dvojici parametru
function GetElapsedTime vraci string s uplynulou dobou
procedura InfoVystup zapise konecne informace do souboru

  OBJEKT TMarqData: funguje stejne jako TMarq, ale je uzpusoben
pro prijimani vypoctenych hodnot zaroven pro vsechny x-ove
namerene hodnoty (vyhoda pro pripad konvoluci ci integralu).
FitFunData: pouziva se misto procedury FitFun a musi byt tedy
vzdy prepsana. Ve vektoru D se vraci spoctene hodnoty pro
vsechny x-ove souradnice.
DoJacData: pouzivaji se misto DoJac. Muze byt prepsana,
standardne je v DoJacData numericka derivace. Kvalitativne
se lisi zahlavi DoJac a DoJacData: TMarq se pocita derivace
podle vsech parametru v danem bode X, tady se chce derivace
podle jednoho parametru ve vsech bodech X.
-------------------------------------------------------------*/


char* TimeToString ( char* S, double time );

#endif

// uplny konec MarqFitP