Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/ALGLIB/lbfgs.cs @ 2575

Last change on this file since 2575 was 2563, checked in by gkronber, 15 years ago

Updated ALGLIB to latest version. #751 (Plugin for for data-modeling with ANN (integrated into CEDMA))

File size: 47.7 KB
Line 
1/*************************************************************************
2Copyright (c) 2007-2008, 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 lbfgs
26    {
27        public struct lbfgsstate
28        {
29            public int n;
30            public int m;
31            public double epsg;
32            public double epsf;
33            public double epsx;
34            public int maxits;
35            public int flags;
36            public int nfev;
37            public int mcstage;
38            public int k;
39            public int q;
40            public int p;
41            public double[] rho;
42            public double[,] y;
43            public double[,] s;
44            public double[] theta;
45            public double[] d;
46            public double stp;
47            public double[] work;
48            public double fold;
49            public double gammak;
50            public double[] x;
51            public double f;
52            public double[] g;
53            public bool xupdated;
54            public AP.rcommstate rstate;
55            public int repiterationscount;
56            public int repnfev;
57            public int repterminationtype;
58            public bool brackt;
59            public bool stage1;
60            public int infoc;
61            public double dg;
62            public double dgm;
63            public double dginit;
64            public double dgtest;
65            public double dgx;
66            public double dgxm;
67            public double dgy;
68            public double dgym;
69            public double finit;
70            public double ftest1;
71            public double fm;
72            public double fx;
73            public double fxm;
74            public double fy;
75            public double fym;
76            public double stx;
77            public double sty;
78            public double stmin;
79            public double stmax;
80            public double width;
81            public double width1;
82            public double xtrapf;
83        };
84
85
86        public struct lbfgsreport
87        {
88            public int iterationscount;
89            public int nfev;
90            public int terminationtype;
91        };
92
93
94
95
96        public const double ftol = 0.0001;
97        public const double xtol = 100*AP.Math.MachineEpsilon;
98        public const double gtol = 0.9;
99        public const int maxfev = 20;
100        public const double stpmin = 1.0E-20;
101        public const double stpmax = 1.0E20;
102
103
104        /*************************************************************************
105                LIMITED MEMORY BFGS METHOD FOR LARGE SCALE OPTIMIZATION
106
107        The subroutine minimizes function F(x) of N arguments by  using  a  quasi-
108        Newton method (LBFGS scheme) which is optimized to use  a  minimum  amount
109        of memory.
110
111        The subroutine generates the approximation of an inverse Hessian matrix by
112        using information about the last M steps of the algorithm  (instead of N).
113        It lessens a required amount of memory from a value  of  order  N^2  to  a
114        value of order 2*N*M.
115
116        Input parameters:
117            N   -   problem dimension. N>0
118            M   -   number of corrections in the BFGS scheme of Hessian
119                    approximation update. Recommended value:  3<=M<=7. The smaller
120                    value causes worse convergence, the bigger will  not  cause  a
121                    considerably better convergence, but will cause a fall in  the
122                    performance. M<=N.
123            X   -   initial solution approximation, array[0..N-1].
124            EpsG -  positive number which  defines  a  precision  of  search.  The
125                    subroutine finishes its work if the condition ||G|| < EpsG  is
126                    satisfied, where ||.|| means Euclidian norm, G - gradient, X -
127                    current approximation.
128            EpsF -  positive number which  defines  a  precision  of  search.  The
129                    subroutine finishes its work if on iteration  number  k+1  the
130                    condition |F(k+1)-F(k)| <= EpsF*max{|F(k)|, |F(k+1)|, 1}    is
131                    satisfied.
132            EpsX -  positive number which  defines  a  precision  of  search.  The
133                    subroutine finishes its work if on iteration number k+1    the
134                    condition |X(k+1)-X(k)| <= EpsX is fulfilled.
135            MaxIts- maximum number of iterations. If MaxIts=0, the number of
136                    iterations is unlimited.
137            Flags - additional settings:
138                    * Flags = 0     means no additional settings
139                    * Flags = 1     "do not allocate memory". used when solving
140                                    a many subsequent tasks with  same N/M  values.
141                                    First  call MUST  be without this flag bit set,
142                                    subsequent calls of MinLBFGS with same LBFGSState
143                                    structure can set Flags to 1.
144
145        Output parameters:
146            State - structure used for reverse communication.
147           
148        See also MinLBFGSIteration, MinLBFGSResults
149
150          -- ALGLIB --
151             Copyright 14.11.2007 by Bochkanov Sergey
152        *************************************************************************/
153        public static void minlbfgs(int n,
154            int m,
155            ref double[] x,
156            double epsg,
157            double epsf,
158            double epsx,
159            int maxits,
160            int flags,
161            ref lbfgsstate state)
162        {
163            bool allocatemem = new bool();
164            int i_ = 0;
165
166            System.Diagnostics.Debug.Assert(n>=1, "MinLBFGS: N too small!");
167            System.Diagnostics.Debug.Assert(m>=1, "MinLBFGS: M too small!");
168            System.Diagnostics.Debug.Assert(m<=n, "MinLBFGS: M too large!");
169            System.Diagnostics.Debug.Assert((double)(epsg)>=(double)(0), "MinLBFGS: negative EpsG!");
170            System.Diagnostics.Debug.Assert((double)(epsf)>=(double)(0), "MinLBFGS: negative EpsF!");
171            System.Diagnostics.Debug.Assert((double)(epsx)>=(double)(0), "MinLBFGS: negative EpsX!");
172            System.Diagnostics.Debug.Assert(maxits>=0, "MinLBFGS: negative MaxIts!");
173           
174            //
175            // Initialize
176            //
177            state.n = n;
178            state.m = m;
179            state.epsg = epsg;
180            state.epsf = epsf;
181            state.epsx = epsx;
182            state.maxits = maxits;
183            state.flags = flags;
184            allocatemem = flags%2==0;
185            flags = flags/2;
186            if( allocatemem )
187            {
188                state.rho = new double[m-1+1];
189                state.theta = new double[m-1+1];
190                state.y = new double[m-1+1, n-1+1];
191                state.s = new double[m-1+1, n-1+1];
192                state.d = new double[n-1+1];
193                state.x = new double[n-1+1];
194                state.g = new double[n-1+1];
195                state.work = new double[n-1+1];
196            }
197           
198            //
199            // Initialize Rep structure
200            //
201            state.xupdated = false;
202           
203            //
204            // Prepare first run
205            //
206            state.k = 0;
207            for(i_=0; i_<=n-1;i_++)
208            {
209                state.x[i_] = x[i_];
210            }
211            state.rstate.ia = new int[6+1];
212            state.rstate.ra = new double[4+1];
213            state.rstate.stage = -1;
214        }
215
216
217        /*************************************************************************
218        One L-BFGS iteration
219
220        Called after initialization with MinLBFGS.
221        See HTML documentation for examples.
222
223        Input parameters:
224            State   -   structure which stores algorithm state between calls and
225                        which is used for reverse communication. Must be initialized
226                        with MinLBFGS.
227
228        If suborutine returned False, iterative proces has converged.
229
230        If subroutine returned True, caller should calculate function value
231        State.F an gradient State.G[0..N-1] at State.X[0..N-1] and call
232        MinLBFGSIteration again.
233
234          -- ALGLIB --
235             Copyright 20.04.2009 by Bochkanov Sergey
236        *************************************************************************/
237        public static bool minlbfgsiteration(ref lbfgsstate state)
238        {
239            bool result = new bool();
240            int n = 0;
241            int m = 0;
242            int maxits = 0;
243            double epsf = 0;
244            double epsg = 0;
245            double epsx = 0;
246            int i = 0;
247            int j = 0;
248            int ic = 0;
249            int mcinfo = 0;
250            double v = 0;
251            double vv = 0;
252            int i_ = 0;
253
254           
255            //
256            // Reverse communication preparations
257            // I know it looks ugly, but it works the same way
258            // anywhere from C++ to Python.
259            //
260            // This code initializes locals by:
261            // * random values determined during code
262            //   generation - on first subroutine call
263            // * values from previous call - on subsequent calls
264            //
265            if( state.rstate.stage>=0 )
266            {
267                n = state.rstate.ia[0];
268                m = state.rstate.ia[1];
269                maxits = state.rstate.ia[2];
270                i = state.rstate.ia[3];
271                j = state.rstate.ia[4];
272                ic = state.rstate.ia[5];
273                mcinfo = state.rstate.ia[6];
274                epsf = state.rstate.ra[0];
275                epsg = state.rstate.ra[1];
276                epsx = state.rstate.ra[2];
277                v = state.rstate.ra[3];
278                vv = state.rstate.ra[4];
279            }
280            else
281            {
282                n = -983;
283                m = -989;
284                maxits = -834;
285                i = 900;
286                j = -287;
287                ic = 364;
288                mcinfo = 214;
289                epsf = -338;
290                epsg = -686;
291                epsx = 912;
292                v = 585;
293                vv = 497;
294            }
295            if( state.rstate.stage==0 )
296            {
297                goto lbl_0;
298            }
299            if( state.rstate.stage==1 )
300            {
301                goto lbl_1;
302            }
303           
304            //
305            // Routine body
306            //
307           
308            //
309            // Unload frequently used variables from State structure
310            // (just for typing convinience)
311            //
312            n = state.n;
313            m = state.m;
314            epsg = state.epsg;
315            epsf = state.epsf;
316            epsx = state.epsx;
317            maxits = state.maxits;
318            state.repterminationtype = 0;
319            state.repiterationscount = 0;
320            state.repnfev = 0;
321           
322            //
323            // Update info
324            //
325            state.xupdated = false;
326           
327            //
328            // Calculate F/G
329            //
330            state.rstate.stage = 0;
331            goto lbl_rcomm;
332        lbl_0:
333            state.repnfev = 1;
334           
335            //
336            // Preparations
337            //
338            state.fold = state.f;
339            v = 0.0;
340            for(i_=0; i_<=n-1;i_++)
341            {
342                v += state.g[i_]*state.g[i_];
343            }
344            v = Math.Sqrt(v);
345            if( (double)(v)==(double)(0) )
346            {
347                state.repterminationtype = 4;
348                result = false;
349                return result;
350            }
351            state.stp = 1.0/v;
352            for(i_=0; i_<=n-1;i_++)
353            {
354                state.d[i_] = -state.g[i_];
355            }
356           
357            //
358            // Main cycle
359            //
360        lbl_2:
361            if( false )
362            {
363                goto lbl_3;
364            }
365           
366            //
367            // Main cycle: prepare to 1-D line search
368            //
369            state.p = state.k%m;
370            state.q = Math.Min(state.k, m-1);
371           
372            //
373            // Store X[k], G[k]
374            //
375            for(i_=0; i_<=n-1;i_++)
376            {
377                state.s[state.p,i_] = -state.x[i_];
378            }
379            for(i_=0; i_<=n-1;i_++)
380            {
381                state.y[state.p,i_] = -state.g[i_];
382            }
383           
384            //
385            // Minimize F(x+alpha*d)
386            //
387            state.mcstage = 0;
388            if( state.k!=0 )
389            {
390                state.stp = 1.0;
391            }
392            mcsrch(n, ref state.x, ref state.f, ref state.g, ref state.d, ref state.stp, ref mcinfo, ref state.nfev, ref state.work, ref state, ref state.mcstage);
393        lbl_4:
394            if( state.mcstage==0 )
395            {
396                goto lbl_5;
397            }
398            state.rstate.stage = 1;
399            goto lbl_rcomm;
400        lbl_1:
401            mcsrch(n, ref state.x, ref state.f, ref state.g, ref state.d, ref state.stp, ref mcinfo, ref state.nfev, ref state.work, ref state, ref state.mcstage);
402            goto lbl_4;
403        lbl_5:
404           
405            //
406            // Main cycle: update information and Hessian.
407            // Check stopping conditions.
408            //
409            state.repnfev = state.repnfev+state.nfev;
410            state.repiterationscount = state.repiterationscount+1;
411           
412            //
413            // Calculate S[k], Y[k], Rho[k], GammaK
414            //
415            for(i_=0; i_<=n-1;i_++)
416            {
417                state.s[state.p,i_] = state.s[state.p,i_] + state.x[i_];
418            }
419            for(i_=0; i_<=n-1;i_++)
420            {
421                state.y[state.p,i_] = state.y[state.p,i_] + state.g[i_];
422            }
423           
424            //
425            // Stopping conditions
426            //
427            if( state.repiterationscount>=maxits & maxits>0 )
428            {
429               
430                //
431                // Too many iterations
432                //
433                state.repterminationtype = 5;
434                result = false;
435                return result;
436            }
437            v = 0.0;
438            for(i_=0; i_<=n-1;i_++)
439            {
440                v += state.g[i_]*state.g[i_];
441            }
442            if( (double)(Math.Sqrt(v))<=(double)(epsg) )
443            {
444               
445                //
446                // Gradient is small enough
447                //
448                state.repterminationtype = 4;
449                result = false;
450                return result;
451            }
452            if( (double)(state.fold-state.f)<=(double)(epsf*Math.Max(Math.Abs(state.fold), Math.Max(Math.Abs(state.f), 1.0))) )
453            {
454               
455                //
456                // F(k+1)-F(k) is small enough
457                //
458                state.repterminationtype = 1;
459                result = false;
460                return result;
461            }
462            v = 0.0;
463            for(i_=0; i_<=n-1;i_++)
464            {
465                v += state.s[state.p,i_]*state.s[state.p,i_];
466            }
467            if( (double)(Math.Sqrt(v))<=(double)(epsx) )
468            {
469               
470                //
471                // X(k+1)-X(k) is small enough
472                //
473                state.repterminationtype = 2;
474                result = false;
475                return result;
476            }
477           
478            //
479            // Calculate Rho[k], GammaK
480            //
481            v = 0.0;
482            for(i_=0; i_<=n-1;i_++)
483            {
484                v += state.y[state.p,i_]*state.s[state.p,i_];
485            }
486            vv = 0.0;
487            for(i_=0; i_<=n-1;i_++)
488            {
489                vv += state.y[state.p,i_]*state.y[state.p,i_];
490            }
491            if( (double)(v)==(double)(0) | (double)(vv)==(double)(0) )
492            {
493               
494                //
495                // Rounding errors make further iterations impossible.
496                //
497                state.repterminationtype = -2;
498                result = false;
499                return result;
500            }
501            state.rho[state.p] = 1/v;
502            state.gammak = v/vv;
503           
504            //
505            //  Calculate d(k+1) = H(k+1)*g(k+1)
506            //
507            //  for I:=K downto K-Q do
508            //      V = s(i)^T * work(iteration:I)
509            //      theta(i) = V
510            //      work(iteration:I+1) = work(iteration:I) - V*Rho(i)*y(i)
511            //  work(last iteration) = H0*work(last iteration)
512            //  for I:=K-Q to K do
513            //      V = y(i)^T*work(iteration:I)
514            //      work(iteration:I+1) = work(iteration:I) +(-V+theta(i))*Rho(i)*s(i)
515            //
516            //  NOW WORK CONTAINS d(k+1)
517            //
518            for(i_=0; i_<=n-1;i_++)
519            {
520                state.work[i_] = state.g[i_];
521            }
522            for(i=state.k; i>=state.k-state.q; i--)
523            {
524                ic = i%m;
525                v = 0.0;
526                for(i_=0; i_<=n-1;i_++)
527                {
528                    v += state.s[ic,i_]*state.work[i_];
529                }
530                state.theta[ic] = v;
531                vv = v*state.rho[ic];
532                for(i_=0; i_<=n-1;i_++)
533                {
534                    state.work[i_] = state.work[i_] - vv*state.y[ic,i_];
535                }
536            }
537            v = state.gammak;
538            for(i_=0; i_<=n-1;i_++)
539            {
540                state.work[i_] = v*state.work[i_];
541            }
542            for(i=state.k-state.q; i<=state.k; i++)
543            {
544                ic = i%m;
545                v = 0.0;
546                for(i_=0; i_<=n-1;i_++)
547                {
548                    v += state.y[ic,i_]*state.work[i_];
549                }
550                vv = state.rho[ic]*(-v+state.theta[ic]);
551                for(i_=0; i_<=n-1;i_++)
552                {
553                    state.work[i_] = state.work[i_] + vv*state.s[ic,i_];
554                }
555            }
556            for(i_=0; i_<=n-1;i_++)
557            {
558                state.d[i_] = -state.work[i_];
559            }
560           
561            //
562            // Next step
563            //
564            state.fold = state.f;
565            state.k = state.k+1;
566            state.xupdated = true;
567            goto lbl_2;
568        lbl_3:
569            result = false;
570            return result;
571           
572            //
573            // Saving state
574            //
575        lbl_rcomm:
576            result = true;
577            state.rstate.ia[0] = n;
578            state.rstate.ia[1] = m;
579            state.rstate.ia[2] = maxits;
580            state.rstate.ia[3] = i;
581            state.rstate.ia[4] = j;
582            state.rstate.ia[5] = ic;
583            state.rstate.ia[6] = mcinfo;
584            state.rstate.ra[0] = epsf;
585            state.rstate.ra[1] = epsg;
586            state.rstate.ra[2] = epsx;
587            state.rstate.ra[3] = v;
588            state.rstate.ra[4] = vv;
589            return result;
590        }
591
592
593        /*************************************************************************
594        L-BFGS algorithm results
595
596        Called after MinLBFGSIteration returned False.
597
598        Input parameters:
599            State   -   algorithm state (used by MinLBFGSIteration).
600
601        Output parameters:
602            X       -   array[0..N-1], solution
603            Rep     -   optimization report:
604                        * Rep.TerminationType completetion code:
605                            * -2    rounding errors prevent further improvement.
606                                    X contains best point found.
607                            * -1    incorrect parameters were specified
608                            *  1    relative function improvement is no more than
609                                    EpsF.
610                            *  2    relative step is no more than EpsX.
611                            *  4    gradient norm is no more than EpsG
612                            *  5    MaxIts steps was taken
613                        * Rep.IterationsCount contains iterations count
614                        * NFEV countains number of function calculations
615
616          -- ALGLIB --
617             Copyright 14.11.2007 by Bochkanov Sergey
618        *************************************************************************/
619        public static void minlbfgsresults(ref lbfgsstate state,
620            ref double[] x,
621            ref lbfgsreport rep)
622        {
623            int i_ = 0;
624
625            x = new double[state.n-1+1];
626            for(i_=0; i_<=state.n-1;i_++)
627            {
628                x[i_] = state.x[i_];
629            }
630            rep.iterationscount = state.repiterationscount;
631            rep.nfev = state.repnfev;
632            rep.terminationtype = state.repterminationtype;
633        }
634
635
636        /*************************************************************************
637        THE  PURPOSE  OF  MCSRCH  IS  TO  FIND A STEP WHICH SATISFIES A SUFFICIENT
638        DECREASE CONDITION AND A CURVATURE CONDITION.
639
640        AT EACH STAGE THE SUBROUTINE  UPDATES  AN  INTERVAL  OF  UNCERTAINTY  WITH
641        ENDPOINTS  STX  AND  STY.  THE INTERVAL OF UNCERTAINTY IS INITIALLY CHOSEN
642        SO THAT IT CONTAINS A MINIMIZER OF THE MODIFIED FUNCTION
643
644            F(X+STP*S) - F(X) - FTOL*STP*(GRADF(X)'S).
645
646        IF  A STEP  IS OBTAINED FOR  WHICH THE MODIFIED FUNCTION HAS A NONPOSITIVE
647        FUNCTION  VALUE  AND  NONNEGATIVE  DERIVATIVE,   THEN   THE   INTERVAL  OF
648        UNCERTAINTY IS CHOSEN SO THAT IT CONTAINS A MINIMIZER OF F(X+STP*S).
649
650        THE  ALGORITHM  IS  DESIGNED TO FIND A STEP WHICH SATISFIES THE SUFFICIENT
651        DECREASE CONDITION
652
653            F(X+STP*S) .LE. F(X) + FTOL*STP*(GRADF(X)'S),
654
655        AND THE CURVATURE CONDITION
656
657            ABS(GRADF(X+STP*S)'S)) .LE. GTOL*ABS(GRADF(X)'S).
658
659        IF  FTOL  IS  LESS  THAN GTOL AND IF, FOR EXAMPLE, THE FUNCTION IS BOUNDED
660        BELOW,  THEN  THERE  IS  ALWAYS  A  STEP  WHICH SATISFIES BOTH CONDITIONS.
661        IF  NO  STEP  CAN BE FOUND  WHICH  SATISFIES  BOTH  CONDITIONS,  THEN  THE
662        ALGORITHM  USUALLY STOPS  WHEN  ROUNDING ERRORS  PREVENT FURTHER PROGRESS.
663        IN THIS CASE STP ONLY SATISFIES THE SUFFICIENT DECREASE CONDITION.
664
665        PARAMETERS DESCRIPRION
666
667        N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER OF VARIABLES.
668
669        X IS  AN  ARRAY  OF  LENGTH N. ON INPUT IT MUST CONTAIN THE BASE POINT FOR
670        THE LINE SEARCH. ON OUTPUT IT CONTAINS X+STP*S.
671
672        F IS  A  VARIABLE. ON INPUT IT MUST CONTAIN THE VALUE OF F AT X. ON OUTPUT
673        IT CONTAINS THE VALUE OF F AT X + STP*S.
674
675        G IS AN ARRAY OF LENGTH N. ON INPUT IT MUST CONTAIN THE GRADIENT OF F AT X.
676        ON OUTPUT IT CONTAINS THE GRADIENT OF F AT X + STP*S.
677
678        S IS AN INPUT ARRAY OF LENGTH N WHICH SPECIFIES THE SEARCH DIRECTION.
679
680        STP  IS  A NONNEGATIVE VARIABLE. ON INPUT STP CONTAINS AN INITIAL ESTIMATE
681        OF A SATISFACTORY STEP. ON OUTPUT STP CONTAINS THE FINAL ESTIMATE.
682
683        FTOL AND GTOL ARE NONNEGATIVE INPUT VARIABLES. TERMINATION OCCURS WHEN THE
684        SUFFICIENT DECREASE CONDITION AND THE DIRECTIONAL DERIVATIVE CONDITION ARE
685        SATISFIED.
686
687        XTOL IS A NONNEGATIVE INPUT VARIABLE. TERMINATION OCCURS WHEN THE RELATIVE
688        WIDTH OF THE INTERVAL OF UNCERTAINTY IS AT MOST XTOL.
689
690        STPMIN AND STPMAX ARE NONNEGATIVE INPUT VARIABLES WHICH SPECIFY LOWER  AND
691        UPPER BOUNDS FOR THE STEP.
692
693        MAXFEV IS A POSITIVE INTEGER INPUT VARIABLE. TERMINATION OCCURS WHEN THE
694        NUMBER OF CALLS TO FCN IS AT LEAST MAXFEV BY THE END OF AN ITERATION.
695
696        INFO IS AN INTEGER OUTPUT VARIABLE SET AS FOLLOWS:
697            INFO = 0  IMPROPER INPUT PARAMETERS.
698
699            INFO = 1  THE SUFFICIENT DECREASE CONDITION AND THE
700                      DIRECTIONAL DERIVATIVE CONDITION HOLD.
701
702            INFO = 2  RELATIVE WIDTH OF THE INTERVAL OF UNCERTAINTY
703                      IS AT MOST XTOL.
704
705            INFO = 3  NUMBER OF CALLS TO FCN HAS REACHED MAXFEV.
706
707            INFO = 4  THE STEP IS AT THE LOWER BOUND STPMIN.
708
709            INFO = 5  THE STEP IS AT THE UPPER BOUND STPMAX.
710
711            INFO = 6  ROUNDING ERRORS PREVENT FURTHER PROGRESS.
712                      THERE MAY NOT BE A STEP WHICH SATISFIES THE
713                      SUFFICIENT DECREASE AND CURVATURE CONDITIONS.
714                      TOLERANCES MAY BE TOO SMALL.
715
716        NFEV IS AN INTEGER OUTPUT VARIABLE SET TO THE NUMBER OF CALLS TO FCN.
717
718        WA IS A WORK ARRAY OF LENGTH N.
719
720        ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JUNE 1983
721        JORGE J. MORE', DAVID J. THUENTE
722        *************************************************************************/
723        private static void mcsrch(int n,
724            ref double[] x,
725            ref double f,
726            ref double[] g,
727            ref double[] s,
728            ref double stp,
729            ref int info,
730            ref int nfev,
731            ref double[] wa,
732            ref lbfgsstate state,
733            ref int stage)
734        {
735            double v = 0;
736            double p5 = 0;
737            double p66 = 0;
738            double zero = 0;
739            int i_ = 0;
740
741           
742            //
743            // init
744            //
745            p5 = 0.5;
746            p66 = 0.66;
747            state.xtrapf = 4.0;
748            zero = 0;
749           
750            //
751            // Main cycle
752            //
753            while( true )
754            {
755                if( stage==0 )
756                {
757                   
758                    //
759                    // NEXT
760                    //
761                    stage = 2;
762                    continue;
763                }
764                if( stage==2 )
765                {
766                    state.infoc = 1;
767                    info = 0;
768                   
769                    //
770                    //     CHECK THE INPUT PARAMETERS FOR ERRORS.
771                    //
772                    if( n<=0 | (double)(stp)<=(double)(0) | (double)(ftol)<(double)(0) | (double)(gtol)<(double)(zero) | (double)(xtol)<(double)(zero) | (double)(stpmin)<(double)(zero) | (double)(stpmax)<(double)(stpmin) | maxfev<=0 )
773                    {
774                        stage = 0;
775                        return;
776                    }
777                   
778                    //
779                    //     COMPUTE THE INITIAL GRADIENT IN THE SEARCH DIRECTION
780                    //     AND CHECK THAT S IS A DESCENT DIRECTION.
781                    //
782                    v = 0.0;
783                    for(i_=0; i_<=n-1;i_++)
784                    {
785                        v += g[i_]*s[i_];
786                    }
787                    state.dginit = v;
788                    if( (double)(state.dginit)>=(double)(0) )
789                    {
790                        stage = 0;
791                        return;
792                    }
793                   
794                    //
795                    //     INITIALIZE LOCAL VARIABLES.
796                    //
797                    state.brackt = false;
798                    state.stage1 = true;
799                    nfev = 0;
800                    state.finit = f;
801                    state.dgtest = ftol*state.dginit;
802                    state.width = stpmax-stpmin;
803                    state.width1 = state.width/p5;
804                    for(i_=0; i_<=n-1;i_++)
805                    {
806                        wa[i_] = x[i_];
807                    }
808                   
809                    //
810                    //     THE VARIABLES STX, FX, DGX CONTAIN THE VALUES OF THE STEP,
811                    //     FUNCTION, AND DIRECTIONAL DERIVATIVE AT THE BEST STEP.
812                    //     THE VARIABLES STY, FY, DGY CONTAIN THE VALUE OF THE STEP,
813                    //     FUNCTION, AND DERIVATIVE AT THE OTHER ENDPOINT OF
814                    //     THE INTERVAL OF UNCERTAINTY.
815                    //     THE VARIABLES STP, F, DG CONTAIN THE VALUES OF THE STEP,
816                    //     FUNCTION, AND DERIVATIVE AT THE CURRENT STEP.
817                    //
818                    state.stx = 0;
819                    state.fx = state.finit;
820                    state.dgx = state.dginit;
821                    state.sty = 0;
822                    state.fy = state.finit;
823                    state.dgy = state.dginit;
824                   
825                    //
826                    // NEXT
827                    //
828                    stage = 3;
829                    continue;
830                }
831                if( stage==3 )
832                {
833                   
834                    //
835                    //     START OF ITERATION.
836                    //
837                    //     SET THE MINIMUM AND MAXIMUM STEPS TO CORRESPOND
838                    //     TO THE PRESENT INTERVAL OF UNCERTAINTY.
839                    //
840                    if( state.brackt )
841                    {
842                        if( (double)(state.stx)<(double)(state.sty) )
843                        {
844                            state.stmin = state.stx;
845                            state.stmax = state.sty;
846                        }
847                        else
848                        {
849                            state.stmin = state.sty;
850                            state.stmax = state.stx;
851                        }
852                    }
853                    else
854                    {
855                        state.stmin = state.stx;
856                        state.stmax = stp+state.xtrapf*(stp-state.stx);
857                    }
858                   
859                    //
860                    //        FORCE THE STEP TO BE WITHIN THE BOUNDS STPMAX AND STPMIN.
861                    //
862                    if( (double)(stp)>(double)(stpmax) )
863                    {
864                        stp = stpmax;
865                    }
866                    if( (double)(stp)<(double)(stpmin) )
867                    {
868                        stp = stpmin;
869                    }
870                   
871                    //
872                    //        IF AN UNUSUAL TERMINATION IS TO OCCUR THEN LET
873                    //        STP BE THE LOWEST POINT OBTAINED SO FAR.
874                    //
875                    if( state.brackt & ((double)(stp)<=(double)(state.stmin) | (double)(stp)>=(double)(state.stmax)) | nfev>=maxfev-1 | state.infoc==0 | state.brackt & (double)(state.stmax-state.stmin)<=(double)(xtol*state.stmax) )
876                    {
877                        stp = state.stx;
878                    }
879                   
880                    //
881                    //        EVALUATE THE FUNCTION AND GRADIENT AT STP
882                    //        AND COMPUTE THE DIRECTIONAL DERIVATIVE.
883                    //
884                    for(i_=0; i_<=n-1;i_++)
885                    {
886                        x[i_] = wa[i_];
887                    }
888                    for(i_=0; i_<=n-1;i_++)
889                    {
890                        x[i_] = x[i_] + stp*s[i_];
891                    }
892                   
893                    //
894                    // NEXT
895                    //
896                    stage = 4;
897                    return;
898                }
899                if( stage==4 )
900                {
901                    info = 0;
902                    nfev = nfev+1;
903                    v = 0.0;
904                    for(i_=0; i_<=n-1;i_++)
905                    {
906                        v += g[i_]*s[i_];
907                    }
908                    state.dg = v;
909                    state.ftest1 = state.finit+stp*state.dgtest;
910                   
911                    //
912                    //        TEST FOR CONVERGENCE.
913                    //
914                    if( state.brackt & ((double)(stp)<=(double)(state.stmin) | (double)(stp)>=(double)(state.stmax)) | state.infoc==0 )
915                    {
916                        info = 6;
917                    }
918                    if( (double)(stp)==(double)(stpmax) & (double)(f)<=(double)(state.ftest1) & (double)(state.dg)<=(double)(state.dgtest) )
919                    {
920                        info = 5;
921                    }
922                    if( (double)(stp)==(double)(stpmin) & ((double)(f)>(double)(state.ftest1) | (double)(state.dg)>=(double)(state.dgtest)) )
923                    {
924                        info = 4;
925                    }
926                    if( nfev>=maxfev )
927                    {
928                        info = 3;
929                    }
930                    if( state.brackt & (double)(state.stmax-state.stmin)<=(double)(xtol*state.stmax) )
931                    {
932                        info = 2;
933                    }
934                    if( (double)(f)<=(double)(state.ftest1) & (double)(Math.Abs(state.dg))<=(double)(-(gtol*state.dginit)) )
935                    {
936                        info = 1;
937                    }
938                   
939                    //
940                    //        CHECK FOR TERMINATION.
941                    //
942                    if( info!=0 )
943                    {
944                        stage = 0;
945                        return;
946                    }
947                   
948                    //
949                    //        IN THE FIRST STAGE WE SEEK A STEP FOR WHICH THE MODIFIED
950                    //        FUNCTION HAS A NONPOSITIVE VALUE AND NONNEGATIVE DERIVATIVE.
951                    //
952                    if( state.stage1 & (double)(f)<=(double)(state.ftest1) & (double)(state.dg)>=(double)(Math.Min(ftol, gtol)*state.dginit) )
953                    {
954                        state.stage1 = false;
955                    }
956                   
957                    //
958                    //        A MODIFIED FUNCTION IS USED TO PREDICT THE STEP ONLY IF
959                    //        WE HAVE NOT OBTAINED A STEP FOR WHICH THE MODIFIED
960                    //        FUNCTION HAS A NONPOSITIVE FUNCTION VALUE AND NONNEGATIVE
961                    //        DERIVATIVE, AND IF A LOWER FUNCTION VALUE HAS BEEN
962                    //        OBTAINED BUT THE DECREASE IS NOT SUFFICIENT.
963                    //
964                    if( state.stage1 & (double)(f)<=(double)(state.fx) & (double)(f)>(double)(state.ftest1) )
965                    {
966                       
967                        //
968                        //           DEFINE THE MODIFIED FUNCTION AND DERIVATIVE VALUES.
969                        //
970                        state.fm = f-stp*state.dgtest;
971                        state.fxm = state.fx-state.stx*state.dgtest;
972                        state.fym = state.fy-state.sty*state.dgtest;
973                        state.dgm = state.dg-state.dgtest;
974                        state.dgxm = state.dgx-state.dgtest;
975                        state.dgym = state.dgy-state.dgtest;
976                       
977                        //
978                        //           CALL CSTEP TO UPDATE THE INTERVAL OF UNCERTAINTY
979                        //           AND TO COMPUTE THE NEW STEP.
980                        //
981                        mcstep(ref state.stx, ref state.fxm, ref state.dgxm, ref state.sty, ref state.fym, ref state.dgym, ref stp, state.fm, state.dgm, ref state.brackt, state.stmin, state.stmax, ref state.infoc);
982                       
983                        //
984                        //           RESET THE FUNCTION AND GRADIENT VALUES FOR F.
985                        //
986                        state.fx = state.fxm+state.stx*state.dgtest;
987                        state.fy = state.fym+state.sty*state.dgtest;
988                        state.dgx = state.dgxm+state.dgtest;
989                        state.dgy = state.dgym+state.dgtest;
990                    }
991                    else
992                    {
993                       
994                        //
995                        //           CALL MCSTEP TO UPDATE THE INTERVAL OF UNCERTAINTY
996                        //           AND TO COMPUTE THE NEW STEP.
997                        //
998                        mcstep(ref state.stx, ref state.fx, ref state.dgx, ref state.sty, ref state.fy, ref state.dgy, ref stp, f, state.dg, ref state.brackt, state.stmin, state.stmax, ref state.infoc);
999                    }
1000                   
1001                    //
1002                    //        FORCE A SUFFICIENT DECREASE IN THE SIZE OF THE
1003                    //        INTERVAL OF UNCERTAINTY.
1004                    //
1005                    if( state.brackt )
1006                    {
1007                        if( (double)(Math.Abs(state.sty-state.stx))>=(double)(p66*state.width1) )
1008                        {
1009                            stp = state.stx+p5*(state.sty-state.stx);
1010                        }
1011                        state.width1 = state.width;
1012                        state.width = Math.Abs(state.sty-state.stx);
1013                    }
1014                   
1015                    //
1016                    //  NEXT.
1017                    //
1018                    stage = 3;
1019                    continue;
1020                }
1021            }
1022        }
1023
1024
1025        private static void mcstep(ref double stx,
1026            ref double fx,
1027            ref double dx,
1028            ref double sty,
1029            ref double fy,
1030            ref double dy,
1031            ref double stp,
1032            double fp,
1033            double dp,
1034            ref bool brackt,
1035            double stmin,
1036            double stmax,
1037            ref int info)
1038        {
1039            bool bound = new bool();
1040            double gamma = 0;
1041            double p = 0;
1042            double q = 0;
1043            double r = 0;
1044            double s = 0;
1045            double sgnd = 0;
1046            double stpc = 0;
1047            double stpf = 0;
1048            double stpq = 0;
1049            double theta = 0;
1050
1051            info = 0;
1052           
1053            //
1054            //     CHECK THE INPUT PARAMETERS FOR ERRORS.
1055            //
1056            if( brackt & ((double)(stp)<=(double)(Math.Min(stx, sty)) | (double)(stp)>=(double)(Math.Max(stx, sty))) | (double)(dx*(stp-stx))>=(double)(0) | (double)(stmax)<(double)(stmin) )
1057            {
1058                return;
1059            }
1060           
1061            //
1062            //     DETERMINE IF THE DERIVATIVES HAVE OPPOSITE SIGN.
1063            //
1064            sgnd = dp*(dx/Math.Abs(dx));
1065           
1066            //
1067            //     FIRST CASE. A HIGHER FUNCTION VALUE.
1068            //     THE MINIMUM IS BRACKETED. IF THE CUBIC STEP IS CLOSER
1069            //     TO STX THAN THE QUADRATIC STEP, THE CUBIC STEP IS TAKEN,
1070            //     ELSE THE AVERAGE OF THE CUBIC AND QUADRATIC STEPS IS TAKEN.
1071            //
1072            if( (double)(fp)>(double)(fx) )
1073            {
1074                info = 1;
1075                bound = true;
1076                theta = 3*(fx-fp)/(stp-stx)+dx+dp;
1077                s = Math.Max(Math.Abs(theta), Math.Max(Math.Abs(dx), Math.Abs(dp)));
1078                gamma = s*Math.Sqrt(AP.Math.Sqr(theta/s)-dx/s*(dp/s));
1079                if( (double)(stp)<(double)(stx) )
1080                {
1081                    gamma = -gamma;
1082                }
1083                p = gamma-dx+theta;
1084                q = gamma-dx+gamma+dp;
1085                r = p/q;
1086                stpc = stx+r*(stp-stx);
1087                stpq = stx+dx/((fx-fp)/(stp-stx)+dx)/2*(stp-stx);
1088                if( (double)(Math.Abs(stpc-stx))<(double)(Math.Abs(stpq-stx)) )
1089                {
1090                    stpf = stpc;
1091                }
1092                else
1093                {
1094                    stpf = stpc+(stpq-stpc)/2;
1095                }
1096                brackt = true;
1097            }
1098            else
1099            {
1100                if( (double)(sgnd)<(double)(0) )
1101                {
1102                   
1103                    //
1104                    //     SECOND CASE. A LOWER FUNCTION VALUE AND DERIVATIVES OF
1105                    //     OPPOSITE SIGN. THE MINIMUM IS BRACKETED. IF THE CUBIC
1106                    //     STEP IS CLOSER TO STX THAN THE QUADRATIC (SECANT) STEP,
1107                    //     THE CUBIC STEP IS TAKEN, ELSE THE QUADRATIC STEP IS TAKEN.
1108                    //
1109                    info = 2;
1110                    bound = false;
1111                    theta = 3*(fx-fp)/(stp-stx)+dx+dp;
1112                    s = Math.Max(Math.Abs(theta), Math.Max(Math.Abs(dx), Math.Abs(dp)));
1113                    gamma = s*Math.Sqrt(AP.Math.Sqr(theta/s)-dx/s*(dp/s));
1114                    if( (double)(stp)>(double)(stx) )
1115                    {
1116                        gamma = -gamma;
1117                    }
1118                    p = gamma-dp+theta;
1119                    q = gamma-dp+gamma+dx;
1120                    r = p/q;
1121                    stpc = stp+r*(stx-stp);
1122                    stpq = stp+dp/(dp-dx)*(stx-stp);
1123                    if( (double)(Math.Abs(stpc-stp))>(double)(Math.Abs(stpq-stp)) )
1124                    {
1125                        stpf = stpc;
1126                    }
1127                    else
1128                    {
1129                        stpf = stpq;
1130                    }
1131                    brackt = true;
1132                }
1133                else
1134                {
1135                    if( (double)(Math.Abs(dp))<(double)(Math.Abs(dx)) )
1136                    {
1137                       
1138                        //
1139                        //     THIRD CASE. A LOWER FUNCTION VALUE, DERIVATIVES OF THE
1140                        //     SAME SIGN, AND THE MAGNITUDE OF THE DERIVATIVE DECREASES.
1141                        //     THE CUBIC STEP IS ONLY USED IF THE CUBIC TENDS TO INFINITY
1142                        //     IN THE DIRECTION OF THE STEP OR IF THE MINIMUM OF THE CUBIC
1143                        //     IS BEYOND STP. OTHERWISE THE CUBIC STEP IS DEFINED TO BE
1144                        //     EITHER STPMIN OR STPMAX. THE QUADRATIC (SECANT) STEP IS ALSO
1145                        //     COMPUTED AND IF THE MINIMUM IS BRACKETED THEN THE THE STEP
1146                        //     CLOSEST TO STX IS TAKEN, ELSE THE STEP FARTHEST AWAY IS TAKEN.
1147                        //
1148                        info = 3;
1149                        bound = true;
1150                        theta = 3*(fx-fp)/(stp-stx)+dx+dp;
1151                        s = Math.Max(Math.Abs(theta), Math.Max(Math.Abs(dx), Math.Abs(dp)));
1152                       
1153                        //
1154                        //        THE CASE GAMMA = 0 ONLY ARISES IF THE CUBIC DOES NOT TEND
1155                        //        TO INFINITY IN THE DIRECTION OF THE STEP.
1156                        //
1157                        gamma = s*Math.Sqrt(Math.Max(0, AP.Math.Sqr(theta/s)-dx/s*(dp/s)));
1158                        if( (double)(stp)>(double)(stx) )
1159                        {
1160                            gamma = -gamma;
1161                        }
1162                        p = gamma-dp+theta;
1163                        q = gamma+(dx-dp)+gamma;
1164                        r = p/q;
1165                        if( (double)(r)<(double)(0) & (double)(gamma)!=(double)(0) )
1166                        {
1167                            stpc = stp+r*(stx-stp);
1168                        }
1169                        else
1170                        {
1171                            if( (double)(stp)>(double)(stx) )
1172                            {
1173                                stpc = stmax;
1174                            }
1175                            else
1176                            {
1177                                stpc = stmin;
1178                            }
1179                        }
1180                        stpq = stp+dp/(dp-dx)*(stx-stp);
1181                        if( brackt )
1182                        {
1183                            if( (double)(Math.Abs(stp-stpc))<(double)(Math.Abs(stp-stpq)) )
1184                            {
1185                                stpf = stpc;
1186                            }
1187                            else
1188                            {
1189                                stpf = stpq;
1190                            }
1191                        }
1192                        else
1193                        {
1194                            if( (double)(Math.Abs(stp-stpc))>(double)(Math.Abs(stp-stpq)) )
1195                            {
1196                                stpf = stpc;
1197                            }
1198                            else
1199                            {
1200                                stpf = stpq;
1201                            }
1202                        }
1203                    }
1204                    else
1205                    {
1206                       
1207                        //
1208                        //     FOURTH CASE. A LOWER FUNCTION VALUE, DERIVATIVES OF THE
1209                        //     SAME SIGN, AND THE MAGNITUDE OF THE DERIVATIVE DOES
1210                        //     NOT DECREASE. IF THE MINIMUM IS NOT BRACKETED, THE STEP
1211                        //     IS EITHER STPMIN OR STPMAX, ELSE THE CUBIC STEP IS TAKEN.
1212                        //
1213                        info = 4;
1214                        bound = false;
1215                        if( brackt )
1216                        {
1217                            theta = 3*(fp-fy)/(sty-stp)+dy+dp;
1218                            s = Math.Max(Math.Abs(theta), Math.Max(Math.Abs(dy), Math.Abs(dp)));
1219                            gamma = s*Math.Sqrt(AP.Math.Sqr(theta/s)-dy/s*(dp/s));
1220                            if( (double)(stp)>(double)(sty) )
1221                            {
1222                                gamma = -gamma;
1223                            }
1224                            p = gamma-dp+theta;
1225                            q = gamma-dp+gamma+dy;
1226                            r = p/q;
1227                            stpc = stp+r*(sty-stp);
1228                            stpf = stpc;
1229                        }
1230                        else
1231                        {
1232                            if( (double)(stp)>(double)(stx) )
1233                            {
1234                                stpf = stmax;
1235                            }
1236                            else
1237                            {
1238                                stpf = stmin;
1239                            }
1240                        }
1241                    }
1242                }
1243            }
1244           
1245            //
1246            //     UPDATE THE INTERVAL OF UNCERTAINTY. THIS UPDATE DOES NOT
1247            //     DEPEND ON THE NEW STEP OR THE CASE ANALYSIS ABOVE.
1248            //
1249            if( (double)(fp)>(double)(fx) )
1250            {
1251                sty = stp;
1252                fy = fp;
1253                dy = dp;
1254            }
1255            else
1256            {
1257                if( (double)(sgnd)<(double)(0.0) )
1258                {
1259                    sty = stx;
1260                    fy = fx;
1261                    dy = dx;
1262                }
1263                stx = stp;
1264                fx = fp;
1265                dx = dp;
1266            }
1267           
1268            //
1269            //     COMPUTE THE NEW STEP AND SAFEGUARD IT.
1270            //
1271            stpf = Math.Min(stmax, stpf);
1272            stpf = Math.Max(stmin, stpf);
1273            stp = stpf;
1274            if( brackt & bound )
1275            {
1276                if( (double)(sty)>(double)(stx) )
1277                {
1278                    stp = Math.Min(stx+0.66*(sty-stx), stp);
1279                }
1280                else
1281                {
1282                    stp = Math.Max(stx+0.66*(sty-stx), stp);
1283                }
1284            }
1285        }
1286    }
1287}
Note: See TracBrowser for help on using the repository browser.