regress.h

Go to the documentation of this file.
00001 /* Gxsm - Gnome X Scanning Microscopy
00002  * universal STM/AFM/SARLS/SPALEED/... controlling and
00003  * data analysis software
00004  * 
00005  * Copyright (C) 1999,2000,2001,2002,2003 Percy Zahl
00006  *
00007  * Authors: Percy Zahl <zahl@users.sf.net>
00008  * additional features: Andreas Klust <klust@users.sf.net>
00009  * WWW Home: http://gxsm.sf.net
00010  *
00011  * This program is free software; you can redistribute it and/or modify
00012  * it under the terms of the GNU General Public License as published by
00013  * the Free Software Foundation; either version 2 of the License, or
00014  * (at your option) any later version.
00015  *
00016  * This program is distributed in the hope that it will be useful,
00017  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00018  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00019  * GNU General Public License for more details.
00020  *
00021  * You should have received a copy of the GNU General Public License
00022  * along with this program; if not, write to the Free Software
00023  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA.
00024  */
00025 
00026 /* -*- Mode: C++; indent-tabs-mode: nil; c-basic-offset: 8 c-style: "K&R" -*- */
00027 
00028 /********************************************************
00029  *                                                      *
00030  *              regress.h :                             *
00031  *              Regression                              *
00032  *              Für Polynome nten Grades                *
00033  *              by Percy Zahl, Version 1.0, 18.2.1994   *
00034  *
00035  *
00036  *
00037  ********************************************************/
00038 
00039 #ifndef __REGRESS__H
00040 #define __REGRESS__H
00041 
00042 #include <iostream>
00043 #include <fstream>
00044 
00045 #include <cstring>
00046 #include <cmath>
00047 #include <cstdlib>
00048 #include <cstdarg>
00049 #include <cstdio>
00050 
00051 #include "mem2d.h"
00052 
00053 
00054 #define         MAXGL   20
00055 #define         MAXVAR  50
00056 #define         Delta   1E-5
00057 #define         NEG             -1
00058 
00059 #define         KEINE_LOESUNG   0
00060 #define         R_IST_LOESUNG   1
00061 #define         SPEZ_LOESUNG    2
00062 
00063 #define         FREI    0
00064 #define         UNFREI  1
00065 #define         EINS    2
00066 
00067 // #define         TALK
00068 
00069 double  gpow(double X, int n);
00070 
00071 // Datenaufnahme und Aufbereitung für LGS
00072 class StatInp;
00073 
00074 // LGS - Behandlung
00075 
00076 // Zeile beinhaltet eine Gleichung
00077 
00078 class Zeile{
00079 private:
00080         int             n,m;
00081         double* gl;     // Pointer auf Koeffizienten und b : a1 x1 + .. + an xn = b
00082 public:
00083         double Wert(int i) { return(gl[i]); };
00084         int Test(int);
00085         void Mult(double);
00086         void Add(class Zeile*);
00087         void P_Zeile();
00088         Zeile(int, int, StatInp*);
00089         Zeile(int, int);
00090         ~Zeile() { delete [] gl; };
00091 };
00092 
00093 class LGSys{
00094 public:
00095         int             n,m;                    // n x m Koeffizienten-Matrix
00096         Zeile*  gln[MAXGL];             // Pointer auf Gleichungen
00097         int             st_fkt[MAXGL];  // Pointer auf Stufenfunktion
00098         int             r;                              // Rang
00099         int*    aus;                    // Auswahl Liste
00100         double* pn[MAXGL];              // Pointer auf Loesungen
00101 
00102         void Eingabe();
00103         void Eingabe(StatInp*);
00104         int SucheGl(int);
00105         int Umformen();
00106         int Untersuche(int);
00107         void Loese(int);
00108         void Berechne(int);
00109         void Ausgabe_Loesung(int);
00110         double CalcPoly(double, int);
00111         double GetPolyK(int);
00112         void Ausgabe_Glsys();
00113         void Tausche(int, int);
00114         int DoIt() { return Untersuche(Umformen()); };
00115         void Init(){
00116           aus=NULL;
00117           for(int i=0; i<MAXGL; ++i){
00118             gln[i]=NULL;
00119             pn[i]=NULL;
00120           }
00121         };
00122         LGSys() {  Init(); Eingabe(); Ausgabe_Loesung(Untersuche(Umformen())); };
00123         LGSys(int N) {  Init(); n = m = N>=MAXGL? MAXGL:N; };
00124         ~LGSys() { 
00125           for(int i=0; i<MAXGL; i++) {
00126             if(gln[i]) delete gln[i]; 
00127             if(pn[i]) delete [] pn[i]; 
00128           }
00129           if(aus) delete [] aus;
00130         };
00131 };
00132 
00133 class StatInp {
00134 public:
00135   StatInp(){
00136     gs=NULL;
00137     x=y=NULL;
00138     std::cout << "Ausgleichs-Polynom berechnen  V1.00 (C) 1994 by PZ\n";
00139     do{ 
00140       std::cout << "\nBis Grad ?   => "; std::cin >> PGrad; PGrad++;
00141     }while(PGrad > MAXGL || PGrad < 0);
00142     PGrad++;
00143     State=0;
00144     N=0;
00145     Read(NULL);
00146                         Optimize();
00147   };
00148   StatInp(char *Grad, char *FName, int yx){
00149     gs=NULL;
00150     x=y=NULL;
00151     State=0;
00152     std::cout << "Ausgleichs-Polynom berechnen  V1.00 (C) 1994 by PZ\n";
00153     PGrad=atoi(Grad)+1;
00154     if(PGrad >= MAXGL || PGrad <= 0)
00155       State=1;
00156     else{
00157       N=0;
00158       Read(FName, yx);
00159       Optimize();
00160     }
00161   };
00162   StatInp(int Grad, double *xp, double *yp, int n){
00163     gs=NULL;
00164     x=y=NULL;
00165     State=0;
00166     PGrad=Grad+1;
00167     if(PGrad >= MAXGL || PGrad <= 0)
00168       State=1;
00169     else{
00170       N=n;
00171       x=new double[N];
00172       y=new double[N];
00173       memcpy(x, xp, sizeof(double)*N);
00174       memcpy(y, yp, sizeof(double)*N);
00175       Optimize();
00176     }
00177   }
00178   StatInp(int Grad, Mem2d *ysrc, int line, 
00179           double (*transfkt)(double,double)=NULL, double a0=0.){
00180     gs=NULL;
00181     x=y=NULL;
00182     State=0;
00183     PGrad=Grad+1;
00184     if(PGrad >= MAXGL || PGrad <= 0)
00185       State=1;
00186     else{
00187       N=ysrc->GetNx();
00188       x=new double[N];
00189       y=new double[N];
00190       ysrc->data->SetPtr(0,line);
00191       if(transfkt)
00192         for(int i=0; i<N; ++i){
00193           x[i] = (double)i;
00194           y[i] = (*transfkt)(ysrc->data->GetNext(),a0);
00195         }
00196       else
00197         for(int i=0; i<N; ++i){
00198           x[i] = (double)i;
00199           y[i] = ysrc->data->GetNext();
00200         }
00201       Optimize();
00202     }
00203   }
00204   ~StatInp(){ if(x) delete [] x; if(y) delete [] y; if(gs) delete gs; };
00205   int Result() { return(State); };
00206   void Read(char *Fn, int yx=0){
00207     double      X,Y;
00208     int i;
00209     char        DatOutFile[100];
00210     double      Xfac, Xof,Yfac,Yof,Err;
00211     char FileName[100], l[1024], *tok;
00212     Xfac=Yfac=1.; Xof=Yof=0.;
00213     Err=0.;
00214     DatOutFile[0]=0;
00215     std::ifstream IDat;
00216     std::ofstream ODat;
00217     if(Fn)
00218       strcpy(FileName, Fn);
00219     else{
00220       std::cout << "Daten File ? => "; std::cin >> FileName;
00221     }
00222     IDat.open(FileName , std::ios::in);
00223     if(!IDat.good()){
00224       std::cout << "\n Kann Datei >>" << FileName << "<< nicht finden !\n";
00225       State=2;
00226       return;
00227     }
00228     while(!IDat.eof()){
00229       IDat.getline(l,1024);
00230       tok=strtok(l," ;,:|&!");
00231       if(*tok == '%' || *tok == '*' || *tok == '#') continue;
00232       if(*tok == '$') break;
00233       if(*tok == 'A') { Xof=atof(tok+2); continue; }
00234       if(*tok == 'B') { Xfac=atof(tok+2); continue; }
00235       if(*tok == 'C') { Yof=atof(tok+2); continue; }
00236       if(*tok == 'D') { Yfac=atof(tok+2); continue; }
00237       if(*tok == 'F') { strcpy(DatOutFile,tok+2); continue; }
00238       if(*tok == 'E') { Err=atof(tok+2); continue; }
00239       N++;      // Eintraege Zaehlen
00240     }
00241     std::cout << "A= " << Xof << " B= " << Xfac << " C= " << Yof << " D= " << Yfac << "\n";
00242     std::cout << "Anzahl Werte: " << N << "\n\n";
00243     IDat.close();
00244     x=new double[N];
00245     y=new double[N];
00246     i=0;
00247     IDat.open(FileName , std::ios::in);
00248     if(DatOutFile[0]) ODat.open(DatOutFile, std::ios::out);
00249     if(!IDat.good()){
00250       std::cout << "\n Kann Datei >>" << FileName << "<< nicht wieder finden !\n";
00251       State=2;
00252       return;
00253     }
00254     std::cout << "\\begin{tabular}{|r|r|}\\hline \\\\\n X   &   Y\\hline\\\\\n";
00255     if(DatOutFile[0])
00256       ODat << "# Regres Daten\n"
00257            << "#\\begin{tabular}{|r|r|}\\hline \\\\\n# X   &   Y  &  Err\\hline\\\\\n";
00258     while(IDat.good() && i < N){
00259       IDat.getline(l,1024);
00260       tok=strtok(l," ;,:&|!");
00261       if(*tok == '%' || *tok == '*' || *tok == '#') continue;
00262       if(*tok == 'A' || *tok == 'B' || *tok == 'C' || *tok== 'D') continue;
00263       if(*tok == 'F' || *tok == 'E' || *tok == 'G' || *tok== 'H') continue;
00264       if(*tok == '$') break;
00265       if(tok){
00266         X = atof(tok)*Xfac+Xof;
00267         tok=strtok(0," ;,&:|!");
00268         if(tok)
00269           Y = atof(tok)*Yfac+Yof;
00270         else{
00271           Y = 0.;
00272           std::cout << "Y-Fehler in Datenzeile: " << i << " !\n";
00273           State=3;
00274         }
00275         if(yx){
00276           y[i]=X; x[i]=Y;
00277         }else{
00278           x[i]=X; y[i]=Y;
00279         }
00280       }
00281       else{
00282         x[i]=y[i]=0.;
00283         std::cout << "X-Fehler in Datenzeile: " << i << " !\n";
00284         State=3;
00285       }
00286       std::cout << " " << x[i] << "  &   " << y[i] << " \\\\\n";
00287       if(DatOutFile[0]) 
00288         ODat << x[i] << " " << y[i] << " " << Err << " \\\\\n";
00289       
00290       ++i;
00291     }
00292     std::cout << " & \\hline\\end{tabular}\n";
00293     if(DatOutFile[0]){ 
00294       ODat << "# & \\hline\\end{tabular}\n";
00295       ODat.close();
00296     }
00297     IDat.close();
00298   }
00299   void Optimize();
00300   double CalcKoef(int k, int j){        // Br.: S.787
00301 /*
00302                            _N_
00303                            \
00304         [ f(x) ] :=         )   f(x )           // Gauß-Klammer
00305                            /       i
00306                            ~~~
00307                            i=1
00308         LGS:
00309 
00310 allgemein:
00311         f(x;a ,a ,...,a   ) = a f (x) + a f (x) + ... + a   f   (x)
00312              0  1      n-1     0 0       1 1             n-1 n-1
00313 
00314         n-1
00315         ___
00316         \
00317          )      a  [ f (x) f (x) ]  =  [ y f (x) ]              (j=0,1,.. n-1
00318         /        k    k     j               j
00319         ~~~
00320         k=0
00321 
00322 für Polynom:
00323                  k
00324         f (x) = x
00325          k
00326 
00327   PGrad-1
00328         ___
00329         \             k+j           j
00330          )      a  [ x   ]  =  [ y x  ]           (j=0,1,...PGrad-1)
00331         /        k
00332         ~~~
00333         k=0
00334 */
00335     int i;
00336     double ckj;
00337     if(k<PGrad) // Linke Seite: Koefizienten
00338       for(ckj=0., i=0; i<N; i++) ckj += gpow(x[i],k+j);
00339     else        // Rechte Seite: Konst.
00340       for(ckj=0., i=0; i<N; i++) ckj += y[i]*gpow(x[i],j);
00341     return ckj;
00342   };
00343 
00344   double GetPoly(double x){
00345     return gs->CalcPoly(x, -1);
00346   };
00347 
00348   double GetPolyKoef(int i){
00349     return gs->GetPolyK(i);
00350   };
00351 
00352 private:
00353   int           N,PGrad;
00354   int           State;
00355   double        *x, *y;
00356   double        a[MAXGL];
00357   LGSys *gs;
00358 };
00359 
00360 #endif

Generated on Sat Apr 1 09:04:03 2006 for GXSM by  doxygen 1.4.6