Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.3.0/ALGLIB-2.3.0/minlm.cs @ 3021

Last change on this file since 3021 was 2806, checked in by gkronber, 15 years ago

Added plugin for new version of ALGLIB. #875 (Update ALGLIB sources)

File size: 35.2 KB
Line 
1/*************************************************************************
2Copyright (c) 2009, Sergey Bochkanov (ALGLIB project).
3
4>>> SOURCE LICENSE >>>
5This program is free software; you can redistribute it and/or modify
6it under the terms of the GNU General Public License as published by
7the Free Software Foundation (www.fsf.org); either version 2 of the
8License, or (at your option) any later version.
9
10This program is distributed in the hope that it will be useful,
11but WITHOUT ANY WARRANTY; without even the implied warranty of
12MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13GNU General Public License for more details.
14
15A copy of the GNU General Public License is available at
16http://www.fsf.org/licensing/licenses
17
18>>> END OF LICENSE >>>
19*************************************************************************/
20
21using System;
22
23namespace alglib
24{
25    public class minlm
26    {
27        public struct lmstate
28        {
29            public bool wrongparams;
30            public int n;
31            public int m;
32            public double epsf;
33            public double epsx;
34            public int maxits;
35            public int flags;
36            public int usermode;
37            public double[] x;
38            public double f;
39            public double[] fi;
40            public double[,] j;
41            public double[,] h;
42            public double[] g;
43            public bool needf;
44            public bool needfg;
45            public bool needfgh;
46            public bool needfij;
47            public bool xupdated;
48            public lbfgs.lbfgsstate internalstate;
49            public lbfgs.lbfgsreport internalrep;
50            public double[] xprec;
51            public double[] xbase;
52            public double[] xdir;
53            public double[] gbase;
54            public double[,] rawmodel;
55            public double[,] model;
56            public double[] work;
57            public AP.rcommstate rstate;
58            public int repiterationscount;
59            public int repterminationtype;
60            public int repnfunc;
61            public int repnjac;
62            public int repngrad;
63            public int repnhess;
64            public int repncholesky;
65            public int solverinfo;
66            public densesolver.densesolverreport solverrep;
67        };
68
69
70        public struct lmreport
71        {
72            public int iterationscount;
73            public int terminationtype;
74            public int nfunc;
75            public int njac;
76            public int ngrad;
77            public int nhess;
78            public int ncholesky;
79        };
80
81
82
83
84        public const int lmmodefj = 0;
85        public const int lmmodefgj = 1;
86        public const int lmmodefgh = 2;
87        public const int lmflagnoprelbfgs = 1;
88        public const int lmflagnointlbfgs = 2;
89        public const int lmprelbfgsm = 5;
90        public const int lmintlbfgsits = 5;
91        public const int lbfgsnorealloc = 1;
92
93
94        /*************************************************************************
95            LEVENBERG-MARQUARDT-LIKE METHOD FOR NON-LINEAR OPTIMIZATION
96
97        Optimization using function gradient and Hessian.  Algorithm -  Levenberg-
98        Marquardt   modification   with   L-BFGS   pre-optimization  and  internal
99        pre-conditioned L-BFGS optimization after each Levenberg-Marquardt step.
100
101        Function F has general form (not "sum-of-squares"):
102
103            F = F(x[0], ..., x[n-1])
104
105        EXAMPLE
106
107        See HTML-documentation.
108
109        INPUT PARAMETERS:
110            N       -   dimension, N>1
111            X       -   initial solution, array[0..N-1]
112            EpsF    -   stopping criterion. Algorithm stops if
113                        |F(k+1)-F(k)| <= EpsF*max{|F(k)|, |F(k+1)|, 1}
114            EpsX    -   stopping criterion. Algorithm stops if
115                        |X(k+1)-X(k)| <= EpsX*(1+|X(k)|)
116            MaxIts  -   stopping criterion. Algorithm stops after MaxIts iterations.
117                        MaxIts=0 means no stopping criterion.
118
119        OUTPUT PARAMETERS:
120            State   -   structure which stores algorithm state between subsequent
121                        calls of MinLMIteration. Used for reverse communication.
122                        This structure should be passed to MinLMIteration subroutine.
123
124        See also MinLMIteration, MinLMResults.
125
126        NOTE
127
128        Passing EpsF=0, EpsX=0 and MaxIts=0 (simultaneously) will lead to automatic
129        stopping criterion selection (small EpsX).
130
131          -- ALGLIB --
132             Copyright 30.03.2009 by Bochkanov Sergey
133        *************************************************************************/
134        public static void minlmfgh(int n,
135            ref double[] x,
136            double epsf,
137            double epsx,
138            int maxits,
139            ref lmstate state)
140        {
141            int i_ = 0;
142
143           
144            //
145            // Prepare RComm
146            //
147            state.rstate.ia = new int[3+1];
148            state.rstate.ba = new bool[0+1];
149            state.rstate.ra = new double[8+1];
150            state.rstate.stage = -1;
151           
152            //
153            // prepare internal structures
154            //
155            lmprepare(n, 0, true, ref state);
156           
157            //
158            // initialize, check parameters
159            //
160            state.xupdated = false;
161            state.n = n;
162            state.m = 0;
163            state.epsf = epsf;
164            state.epsx = epsx;
165            state.maxits = maxits;
166            state.flags = 0;
167            if( (double)(state.epsf)==(double)(0) & (double)(state.epsx)==(double)(0) & state.maxits==0 )
168            {
169                state.epsx = 1.0E-6;
170            }
171            state.usermode = lmmodefgh;
172            state.wrongparams = false;
173            if( n<1 | (double)(epsf)<(double)(0) | (double)(epsx)<(double)(0) | maxits<0 )
174            {
175                state.wrongparams = true;
176                return;
177            }
178            for(i_=0; i_<=n-1;i_++)
179            {
180                state.x[i_] = x[i_];
181            }
182        }
183
184
185        /*************************************************************************
186            LEVENBERG-MARQUARDT-LIKE METHOD FOR NON-LINEAR OPTIMIZATION
187
188        Optimization using function gradient and Jacobian.  Algorithm -  Levenberg-
189        Marquardt   modification   with   L-BFGS   pre-optimization  and  internal
190        pre-conditioned L-BFGS optimization after each Levenberg-Marquardt step.
191
192        Function F is represented as sum of squares:
193
194            F = f[0]^2(x[0],...,x[n-1]) + ... + f[m-1]^2(x[0],...,x[n-1])
195
196        EXAMPLE
197
198        See HTML-documentation.
199
200        INPUT PARAMETERS:
201            N       -   dimension, N>1
202            M       -   number of functions f[i]
203            X       -   initial solution, array[0..N-1]
204            EpsF    -   stopping criterion. Algorithm stops if
205                        |F(k+1)-F(k)| <= EpsF*max{|F(k)|, |F(k+1)|, 1}
206            EpsX    -   stopping criterion. Algorithm stops if
207                        |X(k+1)-X(k)| <= EpsX*(1+|X(k)|)
208            MaxIts  -   stopping criterion. Algorithm stops after MaxIts iterations.
209                        MaxIts=0 means no stopping criterion.
210
211        OUTPUT PARAMETERS:
212            State   -   structure which stores algorithm state between subsequent
213                        calls of MinLMIteration. Used for reverse communication.
214                        This structure should be passed to MinLMIteration subroutine.
215
216        See also MinLMIteration, MinLMResults.
217
218        NOTE
219
220        Passing EpsF=0, EpsX=0 and MaxIts=0 (simultaneously) will lead to automatic
221        stopping criterion selection (small EpsX).
222
223          -- ALGLIB --
224             Copyright 30.03.2009 by Bochkanov Sergey
225        *************************************************************************/
226        public static void minlmfgj(int n,
227            int m,
228            ref double[] x,
229            double epsf,
230            double epsx,
231            int maxits,
232            ref lmstate state)
233        {
234            int i_ = 0;
235
236           
237            //
238            // Prepare RComm
239            //
240            state.rstate.ia = new int[3+1];
241            state.rstate.ba = new bool[0+1];
242            state.rstate.ra = new double[8+1];
243            state.rstate.stage = -1;
244           
245            //
246            // prepare internal structures
247            //
248            lmprepare(n, m, true, ref state);
249           
250            //
251            // initialize, check parameters
252            //
253            state.xupdated = false;
254            state.n = n;
255            state.m = m;
256            state.epsf = epsf;
257            state.epsx = epsx;
258            state.maxits = maxits;
259            state.flags = 0;
260            if( (double)(state.epsf)==(double)(0) & (double)(state.epsx)==(double)(0) & state.maxits==0 )
261            {
262                state.epsx = 1.0E-6;
263            }
264            state.usermode = lmmodefgj;
265            state.wrongparams = false;
266            if( n<1 | m<1 | (double)(epsf)<(double)(0) | (double)(epsx)<(double)(0) | maxits<0 )
267            {
268                state.wrongparams = true;
269                return;
270            }
271            for(i_=0; i_<=n-1;i_++)
272            {
273                state.x[i_] = x[i_];
274            }
275        }
276
277
278        /*************************************************************************
279            CLASSIC LEVENBERG-MARQUARDT METHOD FOR NON-LINEAR OPTIMIZATION
280
281        Optimization using Jacobi matrix. Algorithm  -  classic Levenberg-Marquardt
282        method.
283
284        Function F is represented as sum of squares:
285
286            F = f[0]^2(x[0],...,x[n-1]) + ... + f[m-1]^2(x[0],...,x[n-1])
287
288        EXAMPLE
289
290        See HTML-documentation.
291
292        INPUT PARAMETERS:
293            N       -   dimension, N>1
294            M       -   number of functions f[i]
295            X       -   initial solution, array[0..N-1]
296            EpsF    -   stopping criterion. Algorithm stops if
297                        |F(k+1)-F(k)| <= EpsF*max{|F(k)|, |F(k+1)|, 1}
298            EpsX    -   stopping criterion. Algorithm stops if
299                        |X(k+1)-X(k)| <= EpsX*(1+|X(k)|)
300            MaxIts  -   stopping criterion. Algorithm stops after MaxIts iterations.
301                        MaxIts=0 means no stopping criterion.
302
303        OUTPUT PARAMETERS:
304            State   -   structure which stores algorithm state between subsequent
305                        calls of MinLMIteration. Used for reverse communication.
306                        This structure should be passed to MinLMIteration subroutine.
307
308        See also MinLMIteration, MinLMResults.
309
310        NOTE
311
312        Passing EpsF=0, EpsX=0 and MaxIts=0 (simultaneously) will lead to automatic
313        stopping criterion selection (small EpsX).
314
315          -- ALGLIB --
316             Copyright 30.03.2009 by Bochkanov Sergey
317        *************************************************************************/
318        public static void minlmfj(int n,
319            int m,
320            ref double[] x,
321            double epsf,
322            double epsx,
323            int maxits,
324            ref lmstate state)
325        {
326            int i_ = 0;
327
328           
329            //
330            // Prepare RComm
331            //
332            state.rstate.ia = new int[3+1];
333            state.rstate.ba = new bool[0+1];
334            state.rstate.ra = new double[8+1];
335            state.rstate.stage = -1;
336           
337            //
338            // prepare internal structures
339            //
340            lmprepare(n, m, true, ref state);
341           
342            //
343            // initialize, check parameters
344            //
345            state.xupdated = false;
346            state.n = n;
347            state.m = m;
348            state.epsf = epsf;
349            state.epsx = epsx;
350            state.maxits = maxits;
351            state.flags = 0;
352            if( (double)(state.epsf)==(double)(0) & (double)(state.epsx)==(double)(0) & state.maxits==0 )
353            {
354                state.epsx = 1.0E-6;
355            }
356            state.usermode = lmmodefj;
357            state.wrongparams = false;
358            if( n<1 | m<1 | (double)(epsf)<(double)(0) | (double)(epsx)<(double)(0) | maxits<0 )
359            {
360                state.wrongparams = true;
361                return;
362            }
363            for(i_=0; i_<=n-1;i_++)
364            {
365                state.x[i_] = x[i_];
366            }
367        }
368
369
370        /*************************************************************************
371        One Levenberg-Marquardt iteration.
372
373        Called after inialization of State structure with MinLMXXX subroutine.
374        See HTML docs for examples.
375
376        Input parameters:
377            State   -   structure which stores algorithm state between subsequent
378                        calls and which is used for reverse communication. Must be
379                        initialized with MinLMXXX call first.
380
381        If subroutine returned False, iterative algorithm has converged.
382
383        If subroutine returned True, then:
384        * if State.NeedF=True,      -   function value F at State.X[0..N-1]
385                                        is required
386        * if State.NeedFG=True      -   function value F and gradient G
387                                        are required
388        * if State.NeedFiJ=True     -   function vector f[i] and Jacobi matrix J
389                                        are required
390        * if State.NeedFGH=True     -   function value F, gradient G and Hesian H
391                                        are required
392
393        One and only one of this fields can be set at time.
394
395        Results are stored:
396        * function value            -   in LMState.F
397        * gradient                  -   in LMState.G[0..N-1]
398        * Jacobi matrix             -   in LMState.J[0..M-1,0..N-1]
399        * Hessian                   -   in LMState.H[0..N-1,0..N-1]
400
401          -- ALGLIB --
402             Copyright 10.03.2009 by Bochkanov Sergey
403        *************************************************************************/
404        public static bool minlmiteration(ref lmstate state)
405        {
406            bool result = new bool();
407            int n = 0;
408            int m = 0;
409            int i = 0;
410            double xnorm = 0;
411            double stepnorm = 0;
412            bool spd = new bool();
413            double fbase = 0;
414            double fnew = 0;
415            double lambda = 0;
416            double nu = 0;
417            double lambdaup = 0;
418            double lambdadown = 0;
419            int lbfgsflags = 0;
420            double v = 0;
421            int i_ = 0;
422
423           
424            //
425            // Reverse communication preparations
426            // I know it looks ugly, but it works the same way
427            // anywhere from C++ to Python.
428            //
429            // This code initializes locals by:
430            // * random values determined during code
431            //   generation - on first subroutine call
432            // * values from previous call - on subsequent calls
433            //
434            if( state.rstate.stage>=0 )
435            {
436                n = state.rstate.ia[0];
437                m = state.rstate.ia[1];
438                i = state.rstate.ia[2];
439                lbfgsflags = state.rstate.ia[3];
440                spd = state.rstate.ba[0];
441                xnorm = state.rstate.ra[0];
442                stepnorm = state.rstate.ra[1];
443                fbase = state.rstate.ra[2];
444                fnew = state.rstate.ra[3];
445                lambda = state.rstate.ra[4];
446                nu = state.rstate.ra[5];
447                lambdaup = state.rstate.ra[6];
448                lambdadown = state.rstate.ra[7];
449                v = state.rstate.ra[8];
450            }
451            else
452            {
453                n = -983;
454                m = -989;
455                i = -834;
456                lbfgsflags = 900;
457                spd = true;
458                xnorm = 364;
459                stepnorm = 214;
460                fbase = -338;
461                fnew = -686;
462                lambda = 912;
463                nu = 585;
464                lambdaup = 497;
465                lambdadown = -271;
466                v = -581;
467            }
468            if( state.rstate.stage==0 )
469            {
470                goto lbl_0;
471            }
472            if( state.rstate.stage==1 )
473            {
474                goto lbl_1;
475            }
476            if( state.rstate.stage==2 )
477            {
478                goto lbl_2;
479            }
480            if( state.rstate.stage==3 )
481            {
482                goto lbl_3;
483            }
484            if( state.rstate.stage==4 )
485            {
486                goto lbl_4;
487            }
488            if( state.rstate.stage==5 )
489            {
490                goto lbl_5;
491            }
492            if( state.rstate.stage==6 )
493            {
494                goto lbl_6;
495            }
496           
497            //
498            // Routine body
499            //
500            System.Diagnostics.Debug.Assert(state.usermode==lmmodefj | state.usermode==lmmodefgj | state.usermode==lmmodefgh, "LM: internal error");
501            if( state.wrongparams )
502            {
503                state.repterminationtype = -1;
504                result = false;
505                return result;
506            }
507           
508            //
509            // prepare params
510            //
511            n = state.n;
512            m = state.m;
513            lambdaup = 10;
514            lambdadown = 0.3;
515            nu = 2;
516            lbfgsflags = 0;
517           
518            //
519            // if we have F and G
520            //
521            if( ! ((state.usermode==lmmodefgj | state.usermode==lmmodefgh) & state.flags/lmflagnoprelbfgs%2==0) )
522            {
523                goto lbl_7;
524            }
525           
526            //
527            // First stage of the hybrid algorithm: LBFGS
528            //
529            lbfgs.minlbfgs(n, Math.Min(n, lmprelbfgsm), ref state.x, 0.0, 0.0, 0.0, Math.Max(5, n), 0, ref state.internalstate);
530        lbl_9:
531            if( ! lbfgs.minlbfgsiteration(ref state.internalstate) )
532            {
533                goto lbl_10;
534            }
535           
536            //
537            // RComm
538            //
539            for(i_=0; i_<=n-1;i_++)
540            {
541                state.x[i_] = state.internalstate.x[i_];
542            }
543            lmclearrequestfields(ref state);
544            state.needfg = true;
545            state.rstate.stage = 0;
546            goto lbl_rcomm;
547        lbl_0:
548            state.repnfunc = state.repnfunc+1;
549            state.repngrad = state.repngrad+1;
550           
551            //
552            // Call LBFGS
553            //
554            state.internalstate.f = state.f;
555            for(i_=0; i_<=n-1;i_++)
556            {
557                state.internalstate.g[i_] = state.g[i_];
558            }
559            goto lbl_9;
560        lbl_10:
561            lbfgs.minlbfgsresults(ref state.internalstate, ref state.x, ref state.internalrep);
562        lbl_7:
563           
564            //
565            // Second stage of the hybrid algorithm: LM
566            // Initialize quadratic model.
567            //
568            if( state.usermode!=lmmodefgh )
569            {
570                goto lbl_11;
571            }
572           
573            //
574            // RComm
575            //
576            lmclearrequestfields(ref state);
577            state.needfgh = true;
578            state.rstate.stage = 1;
579            goto lbl_rcomm;
580        lbl_1:
581            state.repnfunc = state.repnfunc+1;
582            state.repngrad = state.repngrad+1;
583            state.repnhess = state.repnhess+1;
584           
585            //
586            // generate raw quadratic model
587            //
588            for(i=0; i<=n-1; i++)
589            {
590                for(i_=0; i_<=n-1;i_++)
591                {
592                    state.rawmodel[i,i_] = state.h[i,i_];
593                }
594            }
595            for(i_=0; i_<=n-1;i_++)
596            {
597                state.gbase[i_] = state.g[i_];
598            }
599            fbase = state.f;
600        lbl_11:
601            if( ! (state.usermode==lmmodefgj | state.usermode==lmmodefj) )
602            {
603                goto lbl_13;
604            }
605           
606            //
607            // RComm
608            //
609            lmclearrequestfields(ref state);
610            state.needfij = true;
611            state.rstate.stage = 2;
612            goto lbl_rcomm;
613        lbl_2:
614            state.repnfunc = state.repnfunc+1;
615            state.repnjac = state.repnjac+1;
616           
617            //
618            // generate raw quadratic model
619            //
620            blas.matrixmatrixmultiply(ref state.j, 0, m-1, 0, n-1, true, ref state.j, 0, m-1, 0, n-1, false, 1.0, ref state.rawmodel, 0, n-1, 0, n-1, 0.0, ref state.work);
621            blas.matrixvectormultiply(ref state.j, 0, m-1, 0, n-1, true, ref state.fi, 0, m-1, 1.0, ref state.gbase, 0, n-1, 0.0);
622            fbase = 0.0;
623            for(i_=0; i_<=m-1;i_++)
624            {
625                fbase += state.fi[i_]*state.fi[i_];
626            }
627        lbl_13:
628            lambda = 0.001;
629        lbl_15:
630            if( false )
631            {
632                goto lbl_16;
633            }
634           
635            //
636            // 1. Model = RawModel+lambda*I
637            // 2. Try to solve (RawModel+Lambda*I)*dx = -g.
638            //    Increase lambda if left part is not positive definite.
639            //
640            for(i=0; i<=n-1; i++)
641            {
642                for(i_=0; i_<=n-1;i_++)
643                {
644                    state.model[i,i_] = state.rawmodel[i,i_];
645                }
646                state.model[i,i] = state.model[i,i]+lambda;
647            }
648            spd = trfac.spdmatrixcholesky(ref state.model, n, true);
649            state.repncholesky = state.repncholesky+1;
650            if( !spd )
651            {
652                lambda = lambda*lambdaup*nu;
653                nu = nu*2;
654                goto lbl_15;
655            }
656            densesolver.spdmatrixcholeskysolve(ref state.model, n, true, ref state.gbase, ref state.solverinfo, ref state.solverrep, ref state.xdir);
657            if( state.solverinfo<0 )
658            {
659                lambda = lambda*lambdaup*nu;
660                nu = nu*2;
661                goto lbl_15;
662            }
663            for(i_=0; i_<=n-1;i_++)
664            {
665                state.xdir[i_] = -1*state.xdir[i_];
666            }
667           
668            //
669            // Candidate lambda found.
670            // 1. Save old w in WBase
671            // 1. Test some stopping criterions
672            // 2. If error(w+wdir)>error(w), increase lambda
673            //
674            for(i_=0; i_<=n-1;i_++)
675            {
676                state.xbase[i_] = state.x[i_];
677            }
678            for(i_=0; i_<=n-1;i_++)
679            {
680                state.x[i_] = state.x[i_] + state.xdir[i_];
681            }
682            xnorm = 0.0;
683            for(i_=0; i_<=n-1;i_++)
684            {
685                xnorm += state.xbase[i_]*state.xbase[i_];
686            }
687            stepnorm = 0.0;
688            for(i_=0; i_<=n-1;i_++)
689            {
690                stepnorm += state.xdir[i_]*state.xdir[i_];
691            }
692            xnorm = Math.Sqrt(xnorm);
693            stepnorm = Math.Sqrt(stepnorm);
694            if( (double)(stepnorm)<=(double)(state.epsx*(1+xnorm)) )
695            {
696               
697                //
698                // step size if small enough
699                //
700                state.repterminationtype = 2;
701                goto lbl_16;
702            }
703            lmclearrequestfields(ref state);
704            state.needf = true;
705            state.rstate.stage = 3;
706            goto lbl_rcomm;
707        lbl_3:
708            state.repnfunc = state.repnfunc+1;
709            fnew = state.f;
710            if( (double)(Math.Abs(fnew-fbase))<=(double)(state.epsf*Math.Max(1, Math.Max(Math.Abs(fbase), Math.Abs(fnew)))) )
711            {
712               
713                //
714                // function change is small enough
715                //
716                state.repterminationtype = 1;
717                goto lbl_16;
718            }
719            if( (double)(fnew)>(double)(fbase) )
720            {
721               
722                //
723                // restore state and continue out search for lambda
724                //
725                for(i_=0; i_<=n-1;i_++)
726                {
727                    state.x[i_] = state.xbase[i_];
728                }
729                lambda = lambda*lambdaup*nu;
730                nu = nu*2;
731                goto lbl_15;
732            }
733            if( ! ((state.usermode==lmmodefgj | state.usermode==lmmodefgh) & state.flags/lmflagnointlbfgs%2==0) )
734            {
735                goto lbl_17;
736            }
737           
738            //
739            // Optimize using inv(cholesky(H)) as preconditioner
740            //
741            if( ! trinverse.rmatrixtrinverse(ref state.model, n, true, false) )
742            {
743                goto lbl_19;
744            }
745           
746            //
747            // if matrix can be inverted use it.
748            // just silently move to next iteration otherwise.
749            // (will be very, very rare, mostly for specially
750            // designed near-degenerate tasks)
751            //
752            for(i_=0; i_<=n-1;i_++)
753            {
754                state.xbase[i_] = state.x[i_];
755            }
756            for(i=0; i<=n-1; i++)
757            {
758                state.xprec[i] = 0;
759            }
760            lbfgs.minlbfgs(n, Math.Min(n, lmintlbfgsits), ref state.xprec, 0.0, 0.0, 0.0, lmintlbfgsits, lbfgsflags, ref state.internalstate);
761        lbl_21:
762            if( ! lbfgs.minlbfgsiteration(ref state.internalstate) )
763            {
764                goto lbl_22;
765            }
766           
767            //
768            // convert XPrec to unpreconditioned form, then call RComm.
769            //
770            for(i=0; i<=n-1; i++)
771            {
772                v = 0.0;
773                for(i_=i; i_<=n-1;i_++)
774                {
775                    v += state.internalstate.x[i_]*state.model[i,i_];
776                }
777                state.x[i] = state.xbase[i]+v;
778            }
779            lmclearrequestfields(ref state);
780            state.needfg = true;
781            state.rstate.stage = 4;
782            goto lbl_rcomm;
783        lbl_4:
784            state.repnfunc = state.repnfunc+1;
785            state.repngrad = state.repngrad+1;
786           
787            //
788            // 1. pass State.F to State.InternalState.F
789            // 2. convert gradient back to preconditioned form
790            //
791            state.internalstate.f = state.f;
792            for(i=0; i<=n-1; i++)
793            {
794                state.internalstate.g[i] = 0;
795            }
796            for(i=0; i<=n-1; i++)
797            {
798                v = state.g[i];
799                for(i_=i; i_<=n-1;i_++)
800                {
801                    state.internalstate.g[i_] = state.internalstate.g[i_] + v*state.model[i,i_];
802                }
803            }
804           
805            //
806            // next iteration
807            //
808            goto lbl_21;
809        lbl_22:
810           
811            //
812            // change LBFGS flags to NoRealloc.
813            // L-BFGS subroutine will use memory allocated from previous run.
814            // it is possible since all subsequent calls will be with same N/M.
815            //
816            lbfgsflags = lbfgsnorealloc;
817           
818            //
819            // back to unpreconditioned X
820            //
821            lbfgs.minlbfgsresults(ref state.internalstate, ref state.xprec, ref state.internalrep);
822            for(i=0; i<=n-1; i++)
823            {
824                v = 0.0;
825                for(i_=i; i_<=n-1;i_++)
826                {
827                    v += state.xprec[i_]*state.model[i,i_];
828                }
829                state.x[i] = state.xbase[i]+v;
830            }
831        lbl_19:
832        lbl_17:
833           
834            //
835            // Accept new position.
836            // Calculate Hessian
837            //
838            if( state.usermode!=lmmodefgh )
839            {
840                goto lbl_23;
841            }
842           
843            //
844            // RComm
845            //
846            lmclearrequestfields(ref state);
847            state.needfgh = true;
848            state.rstate.stage = 5;
849            goto lbl_rcomm;
850        lbl_5:
851            state.repnfunc = state.repnfunc+1;
852            state.repngrad = state.repngrad+1;
853            state.repnhess = state.repnhess+1;
854           
855            //
856            // Update raw quadratic model
857            //
858            for(i=0; i<=n-1; i++)
859            {
860                for(i_=0; i_<=n-1;i_++)
861                {
862                    state.rawmodel[i,i_] = state.h[i,i_];
863                }
864            }
865            for(i_=0; i_<=n-1;i_++)
866            {
867                state.gbase[i_] = state.g[i_];
868            }
869            fbase = state.f;
870        lbl_23:
871            if( ! (state.usermode==lmmodefgj | state.usermode==lmmodefj) )
872            {
873                goto lbl_25;
874            }
875           
876            //
877            // RComm
878            //
879            lmclearrequestfields(ref state);
880            state.needfij = true;
881            state.rstate.stage = 6;
882            goto lbl_rcomm;
883        lbl_6:
884            state.repnfunc = state.repnfunc+1;
885            state.repnjac = state.repnjac+1;
886           
887            //
888            // generate raw quadratic model
889            //
890            blas.matrixmatrixmultiply(ref state.j, 0, m-1, 0, n-1, true, ref state.j, 0, m-1, 0, n-1, false, 1.0, ref state.rawmodel, 0, n-1, 0, n-1, 0.0, ref state.work);
891            blas.matrixvectormultiply(ref state.j, 0, m-1, 0, n-1, true, ref state.fi, 0, m-1, 1.0, ref state.gbase, 0, n-1, 0.0);
892            fbase = 0.0;
893            for(i_=0; i_<=m-1;i_++)
894            {
895                fbase += state.fi[i_]*state.fi[i_];
896            }
897        lbl_25:
898            state.repiterationscount = state.repiterationscount+1;
899            if( state.repiterationscount>=state.maxits & state.maxits>0 )
900            {
901                state.repterminationtype = 5;
902                goto lbl_16;
903            }
904           
905            //
906            // Update lambda
907            //
908            lambda = lambda*lambdadown;
909            nu = 2;
910            goto lbl_15;
911        lbl_16:
912            result = false;
913            return result;
914           
915            //
916            // Saving state
917            //
918        lbl_rcomm:
919            result = true;
920            state.rstate.ia[0] = n;
921            state.rstate.ia[1] = m;
922            state.rstate.ia[2] = i;
923            state.rstate.ia[3] = lbfgsflags;
924            state.rstate.ba[0] = spd;
925            state.rstate.ra[0] = xnorm;
926            state.rstate.ra[1] = stepnorm;
927            state.rstate.ra[2] = fbase;
928            state.rstate.ra[3] = fnew;
929            state.rstate.ra[4] = lambda;
930            state.rstate.ra[5] = nu;
931            state.rstate.ra[6] = lambdaup;
932            state.rstate.ra[7] = lambdadown;
933            state.rstate.ra[8] = v;
934            return result;
935        }
936
937
938        /*************************************************************************
939        Levenberg-Marquardt algorithm results
940
941        Called after MinLMIteration returned False.
942
943        Input parameters:
944            State   -   algorithm state (used by MinLMIteration).
945
946        Output parameters:
947            X       -   array[0..N-1], solution
948            Rep     -   optimization report:
949                        * Rep.TerminationType completetion code:
950                            * -1    incorrect parameters were specified
951                            *  1    relative function improvement is no more than
952                                    EpsF.
953                            *  2    relative step is no more than EpsX.
954                            *  5    MaxIts steps was taken
955                        * Rep.IterationsCount contains iterations count
956                        * Rep.NFunc     - number of function calculations
957                        * Rep.NJac      - number of Jacobi matrix calculations
958                        * Rep.NGrad     - number of gradient calculations
959                        * Rep.NHess     - number of Hessian calculations
960                        * Rep.NCholesky - number of Cholesky decomposition calculations
961
962          -- ALGLIB --
963             Copyright 10.03.2009 by Bochkanov Sergey
964        *************************************************************************/
965        public static void minlmresults(ref lmstate state,
966            ref double[] x,
967            ref lmreport rep)
968        {
969            int i_ = 0;
970
971            x = new double[state.n-1+1];
972            for(i_=0; i_<=state.n-1;i_++)
973            {
974                x[i_] = state.x[i_];
975            }
976            rep.iterationscount = state.repiterationscount;
977            rep.terminationtype = state.repterminationtype;
978            rep.nfunc = state.repnfunc;
979            rep.njac = state.repnjac;
980            rep.ngrad = state.repngrad;
981            rep.nhess = state.repnhess;
982            rep.ncholesky = state.repncholesky;
983        }
984
985
986        /*************************************************************************
987        Prepare internal structures (except for RComm).
988
989        Note: M must be zero for FGH mode, non-zero for FJ/FGJ mode.
990        *************************************************************************/
991        private static void lmprepare(int n,
992            int m,
993            bool havegrad,
994            ref lmstate state)
995        {
996            state.repiterationscount = 0;
997            state.repterminationtype = 0;
998            state.repnfunc = 0;
999            state.repnjac = 0;
1000            state.repngrad = 0;
1001            state.repnhess = 0;
1002            state.repncholesky = 0;
1003            if( n<0 | m<0 )
1004            {
1005                return;
1006            }
1007            if( havegrad )
1008            {
1009                state.g = new double[n-1+1];
1010            }
1011            if( m!=0 )
1012            {
1013                state.j = new double[m-1+1, n-1+1];
1014                state.fi = new double[m-1+1];
1015                state.h = new double[0+1, 0+1];
1016            }
1017            else
1018            {
1019                state.j = new double[0+1, 0+1];
1020                state.fi = new double[0+1];
1021                state.h = new double[n-1+1, n-1+1];
1022            }
1023            state.x = new double[n-1+1];
1024            state.rawmodel = new double[n-1+1, n-1+1];
1025            state.model = new double[n-1+1, n-1+1];
1026            state.xbase = new double[n-1+1];
1027            state.xprec = new double[n-1+1];
1028            state.gbase = new double[n-1+1];
1029            state.xdir = new double[n-1+1];
1030            state.work = new double[Math.Max(n, m)+1];
1031        }
1032
1033
1034        /*************************************************************************
1035        Clears request fileds (to be sure that we don't forgot to clear something)
1036        *************************************************************************/
1037        private static void lmclearrequestfields(ref lmstate state)
1038        {
1039            state.needf = false;
1040            state.needfg = false;
1041            state.needfgh = false;
1042            state.needfij = false;
1043        }
1044    }
1045}
Note: See TracBrowser for help on using the repository browser.