Planet
navi homePPSaboutscreenshotsdownloaddevelopmentforum

source: orxonox.OLD/orxonox/trunk/src/lib/newmat/newmatnl.cpp @ 4615

Last change on this file since 4615 was 4565, checked in by patrick, 19 years ago

orxonox/trunk: added the newmat library to the project. needs some translation in directory, temp under util/newmat. is needed by the collision detection engine to perform lin alg operations such as eigenvector decomposition. perhaps we will make our own library to do that later.

File size: 7.0 KB
Line 
1//$$ newmatnl.cpp         Non-linear optimisation
2
3// Copyright (C) 1993,4,5,6: R B Davies
4
5
6#define WANT_MATH
7#define WANT_STREAM
8
9#include "newmatap.h"
10#include "newmatnl.h"
11
12#ifdef use_namespace
13namespace NEWMAT {
14#endif
15
16
17
18void FindMaximum2::Fit(ColumnVector& Theta, int n_it)
19{
20   Tracer tr("FindMaximum2::Fit");
21   enum State {Start, Restart, Continue, Interpolate, Extrapolate,
22      Fail, Convergence};
23   State TheState = Start;
24   Real z,w,x,x2,g,l1,l2,l3,d1,d2=0,d3;
25   ColumnVector Theta1, Theta2, Theta3;
26   int np = Theta.Nrows();
27   ColumnVector H1(np), H3, HP(np), K, K1(np);
28   bool oorg, conv;
29   int counter = 0;
30   Theta1 = Theta; HP = 0.0; g = 0.0;
31
32   // This is really a set of gotos and labels, but they do not work
33   // correctly in AT&T C++ and Sun 4.01 C++.
34
35   for(;;)
36   {
37      switch (TheState)
38      {
39      case Start:
40         tr.ReName("FindMaximum2::Fit/Start");
41         Value(Theta1, true, l1, oorg);
42         if (oorg) Throw(ProgramException("invalid starting value\n"));
43
44      case Restart:
45         tr.ReName("FindMaximum2::Fit/ReStart");
46         conv = NextPoint(H1, d1);
47         if (conv) { TheState = Convergence; break; }
48         if (counter++ > n_it) { TheState = Fail; break; }
49
50         z = 1.0 / sqrt(d1);
51         H3 = H1 * z; K = (H3 - HP) * g; HP = H3;
52         g = 0.0;                     // de-activate to use curved projection
53         if (g==0.0) K1 = 0.0; else K1 = K * 0.2 + K1 * 0.6;
54         // (K - K1) * alpha + K1 * (1 - alpha)
55         //     = K * alpha + K1 * (1 - 2 * alpha)
56         K = K1 * d1; g = z;
57
58      case Continue:
59         tr.ReName("FindMaximum2::Fit/Continue");
60         Theta2 = Theta1 + H1 + K;
61         Value(Theta2, false, l2, oorg);
62         if (counter++ > n_it) { TheState = Fail; break; }
63         if (oorg)
64         {
65            H1 *= 0.5; K *= 0.25; d1 *= 0.5; g *= 2.0;
66            TheState =  Continue; break;
67         }
68         d2 = LastDerivative(H1 + K * 2.0);
69
70      case Interpolate:
71         tr.ReName("FindMaximum2::Fit/Interpolate");
72         z = d1 + d2 - 3.0 * (l2 - l1);
73         w = z * z - d1 * d2;
74         if (w < 0.0) { TheState = Extrapolate; break; }
75         w = z + sqrt(w);
76         if (1.5 * w + d1 < 0.0)
77            { TheState = Extrapolate; break; }
78         if (d2 > 0.0 && l2 > l1 && w > 0.0)
79            { TheState = Extrapolate; break; }
80         x = d1 / (w + d1); x2 = x * x; g /= x;
81         Theta3 = Theta1 + H1 * x + K * x2;
82         Value(Theta3, true, l3, oorg);
83         if (counter++ > n_it) { TheState = Fail; break; }
84         if (oorg)
85         {
86            if (x <= 1.0)
87               { x *= 0.5; x2 = x*x; g *= 2.0; d1 *= x; H1 *= x; K *= x2; }
88            else
89            {
90               x = 0.5 * (x-1.0); x2 = x*x; Theta1 = Theta2;
91               H1 = (H1 + K * 2.0) * x;
92               K *= x2; g = 0.0; d1 = x * d2; l1 = l2;
93            }
94            TheState = Continue; break;
95         }
96
97         if (l3 >= l1 && l3 >= l2)
98            { Theta1 = Theta3; l1 = l3; TheState =  Restart; break; }
99
100         d3 = LastDerivative(H1 + K * 2.0);
101         if (l1 > l2)
102            { H1 *= x; K *= x2; Theta2 = Theta3; d1 *= x; d2 = d3*x; }
103         else
104         {
105            Theta1 = Theta2; Theta2 = Theta3;
106            x -= 1.0; x2 = x*x; g = 0.0; H1 = (H1 + K * 2.0) * x;
107            K *= x2; l1 = l2; l2 = l3; d1 = x*d2; d2 = x*d3;
108            if (d1 <= 0.0) { TheState = Start; break; }
109         }
110         TheState =  Interpolate; break;
111
112      case Extrapolate:
113         tr.ReName("FindMaximum2::Fit/Extrapolate");
114         Theta1 = Theta2; g = 0.0; K *= 4.0; H1 = (H1 * 2.0 + K);
115         d1 = 2.0 * d2; l1 = l2;
116         TheState = Continue; break;
117
118      case Fail:
119         Throw(ConvergenceException(Theta));
120
121      case Convergence:
122         Theta = Theta1; return;
123      }
124   }
125}
126
127
128
129void NonLinearLeastSquares::Value
130   (const ColumnVector& Parameters, bool, Real& v, bool& oorg)
131{
132   Tracer tr("NonLinearLeastSquares::Value");
133   Y.ReSize(n_obs); X.ReSize(n_obs,n_param);
134   // put the fitted values in Y, the derivatives in X.
135   Pred.Set(Parameters);
136   if (!Pred.IsValid()) { oorg=true; return; }
137   for (int i=1; i<=n_obs; i++)
138   {
139      Y(i) = Pred(i);
140      X.Row(i) = Pred.Derivatives();
141   }
142   if (!Pred.IsValid()) { oorg=true; return; }  // check afterwards as well
143   Y = *DataPointer - Y; Real ssq = Y.SumSquare();
144   errorvar =  ssq / (n_obs - n_param);
145   cout << "\n" << setw(15) << setprecision(10) << " " << errorvar;
146   Derivs = Y.t() * X;          // get the derivative and stash it
147   oorg = false; v = -0.5 * ssq;
148}
149
150bool NonLinearLeastSquares::NextPoint(ColumnVector& Adj, Real& test)
151{
152   Tracer tr("NonLinearLeastSquares::NextPoint");
153   QRZ(X, U); QRZ(X, Y, M);     // do the QR decomposition
154   test = M.SumSquare();
155   cout << " " << setw(15) << setprecision(10)
156      << test << " " << Y.SumSquare() / (n_obs - n_param);
157   Adj = U.i() * M;
158   if (test < errorvar * criterion) return true;
159   else return false;
160}
161
162Real NonLinearLeastSquares::LastDerivative(const ColumnVector& H)
163{ return (Derivs * H).AsScalar(); }
164
165void NonLinearLeastSquares::Fit(const ColumnVector& Data,
166   ColumnVector& Parameters)
167{
168   Tracer tr("NonLinearLeastSquares::Fit");
169   n_param = Parameters.Nrows(); n_obs = Data.Nrows();
170   DataPointer = &Data;
171   FindMaximum2::Fit(Parameters, Lim);
172   cout << "\nConverged\n";
173}
174
175void NonLinearLeastSquares::MakeCovariance()
176{
177   if (Covariance.Nrows()==0)
178   {
179      UpperTriangularMatrix UI = U.i();
180      Covariance << UI * UI.t() * errorvar;
181      SE << Covariance;                 // get diagonals
182      for (int i = 1; i<=n_param; i++) SE(i) = sqrt(SE(i));
183   }
184}
185
186void NonLinearLeastSquares::GetStandardErrors(ColumnVector& SEX)
187   { MakeCovariance(); SEX = SE.AsColumn(); }
188
189void NonLinearLeastSquares::GetCorrelations(SymmetricMatrix& Corr)
190   { MakeCovariance(); Corr << SE.i() * Covariance * SE.i(); }
191
192void NonLinearLeastSquares::GetHatDiagonal(DiagonalMatrix& Hat) const
193{
194   Hat.ReSize(n_obs);
195   for (int i = 1; i<=n_obs; i++) Hat(i) = X.Row(i).SumSquare();
196}
197
198
199// the MLE_D_FI routines
200
201void MLE_D_FI::Value
202   (const ColumnVector& Parameters, bool wg, Real& v, bool& oorg)
203{
204   Tracer tr("MLE_D_FI::Value");
205   if (!LL.IsValid(Parameters,wg)) { oorg=true; return; }
206   v = LL.LogLikelihood();
207   if (!LL.IsValid()) { oorg=true; return; }     // check validity again
208   cout << "\n" << setw(20) << setprecision(10) << v;
209   oorg = false;
210   Derivs = LL.Derivatives();                    // Get derivatives
211}
212
213bool MLE_D_FI::NextPoint(ColumnVector& Adj, Real& test)
214{
215   Tracer tr("MLE_D_FI::NextPoint");
216   SymmetricMatrix FI = LL.FI();
217   LT = Cholesky(FI);
218   ColumnVector Adj1 = LT.i() * Derivs;
219   Adj = LT.t().i() * Adj1;
220   test = SumSquare(Adj1);
221   cout << "   " << setw(20) << setprecision(10) << test;
222   return (test < Criterion);
223}
224
225Real MLE_D_FI::LastDerivative(const ColumnVector& H)
226{ return (Derivs.t() * H).AsScalar(); }
227
228void MLE_D_FI::Fit(ColumnVector& Parameters)
229{
230   Tracer tr("MLE_D_FI::Fit");
231   FindMaximum2::Fit(Parameters,Lim);
232   cout << "\nConverged\n";
233}
234 
235void MLE_D_FI::MakeCovariance()
236{
237   if (Covariance.Nrows()==0)
238   {
239      LowerTriangularMatrix LTI = LT.i();
240      Covariance << LTI.t() * LTI;
241      SE << Covariance;                // get diagonal
242      int n = Covariance.Nrows();
243      for (int i=1; i <= n; i++) SE(i) = sqrt(SE(i));
244   }
245}
246
247void MLE_D_FI::GetStandardErrors(ColumnVector& SEX)
248{ MakeCovariance(); SEX = SE.AsColumn(); }
249   
250void MLE_D_FI::GetCorrelations(SymmetricMatrix& Corr)
251{ MakeCovariance(); Corr << SE.i() * Covariance * SE.i(); }
252
253
254
255#ifdef use_namespace
256}
257#endif
258
Note: See TracBrowser for help on using the repository browser.