Assembla home | Assembla project page
 

Changeset 204

Show
Ignore:
Timestamp:
06/07/08 18:11:00 (4 months ago)
Author:
snovik
Message:

New: Tests of Optimization
Fix: Optimzations

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

  • trunk/QLNet/QLNet.vsmdi

    r203 r204  
    1010  <TestList name="Swap" id="1fb9da58-c37b-40da-9772-c8899ed4f8c6" parentListId="8c43106b-9dc1-4907-a29f-aa66a61bf5b6"> 
    1111    <TestLinks> 
    12       <TestLink id="bfdd9348-cf70-0cc9-aefe-30e8f235d572" name="testSpreadDependency" storage="test2008\bin\debug\test2008.dll" type="Microsoft.VisualStudio.TestTools.TestTypes.Unit.UnitTestElement, Microsoft.VisualStudio.QualityTools.Tips.UnitTest.ObjectModel,   PublicKeyToken=b03f5f7f11d50a3a" /> 
     12      <TestLink id="019427b9-dea4-5598-dd43-a00afdcf7984" name="testCachedValue" storage="test2008\bin\debug\test2008.dll" type="Microsoft.VisualStudio.TestTools.TestTypes.Unit.UnitTestElement, Microsoft.VisualStudio.QualityTools.Tips.UnitTest.ObjectModel,   PublicKeyToken=b03f5f7f11d50a3a" /> 
    1313      <TestLink id="cd832602-6407-a10e-5585-88a96542e91e" name="testInArrears" storage="test2008\bin\debug\test2008.dll" type="Microsoft.VisualStudio.TestTools.TestTypes.Unit.UnitTestElement, Microsoft.VisualStudio.QualityTools.Tips.UnitTest.ObjectModel,   PublicKeyToken=b03f5f7f11d50a3a" /> 
    1414      <TestLink id="4a7f54ba-aa40-17eb-4f2c-4e956babf8b3" name="testRateDependency" storage="test2008\bin\debug\test2008.dll" type="Microsoft.VisualStudio.TestTools.TestTypes.Unit.UnitTestElement, Microsoft.VisualStudio.QualityTools.Tips.UnitTest.ObjectModel,   PublicKeyToken=b03f5f7f11d50a3a" /> 
    1515      <TestLink id="6d8e3375-fefb-38b0-2f10-3bb7dcdc7cc6" name="testFairRate" storage="test2008\bin\debug\test2008.dll" type="Microsoft.VisualStudio.TestTools.TestTypes.Unit.UnitTestElement, Microsoft.VisualStudio.QualityTools.Tips.UnitTest.ObjectModel,   PublicKeyToken=b03f5f7f11d50a3a" /> 
    16       <TestLink id="019427b9-dea4-5598-dd43-a00afdcf7984" name="testCachedValue" storage="test2008\bin\debug\test2008.dll" type="Microsoft.VisualStudio.TestTools.TestTypes.Unit.UnitTestElement, Microsoft.VisualStudio.QualityTools.Tips.UnitTest.ObjectModel,   PublicKeyToken=b03f5f7f11d50a3a" /> 
     16      <TestLink id="bfdd9348-cf70-0cc9-aefe-30e8f235d572" name="testSpreadDependency" storage="test2008\bin\debug\test2008.dll" type="Microsoft.VisualStudio.TestTools.TestTypes.Unit.UnitTestElement, Microsoft.VisualStudio.QualityTools.Tips.UnitTest.ObjectModel,   PublicKeyToken=b03f5f7f11d50a3a" /> 
    1717      <TestLink id="f909c484-4167-4782-e0bc-cd5f18e918ba" name="testFairSpread" storage="test2008\bin\debug\test2008.dll" type="Microsoft.VisualStudio.TestTools.TestTypes.Unit.UnitTestElement, Microsoft.VisualStudio.QualityTools.Tips.UnitTest.ObjectModel,   PublicKeyToken=b03f5f7f11d50a3a" /> 
    1818    </TestLinks> 
    …… 
    7676    </TestLinks> 
    7777  </TestList> 
     78  <TestList name="Optimizations" id="a9871493-94b3-4f38-87f8-daa5d63d8c21" parentListId="8c43106b-9dc1-4907-a29f-aa66a61bf5b6"> 
     79    <TestLinks> 
     80      <TestLink id="ae1cd6c7-4ff4-b62a-9381-2a05f7f7cbf0" name="nestedOptimizationTest" storage="test2008\bin\debug\test2008.dll" type="Microsoft.VisualStudio.TestTools.TestTypes.Unit.UnitTestElement, Microsoft.VisualStudio.QualityTools.Tips.UnitTest.ObjectModel,   PublicKeyToken=b03f5f7f11d50a3a" /> 
     81      <TestLink id="d534c5b9-8146-acd9-ee15-708139fcfc23" name="OptimizersTest" storage="test2008\bin\debug\test2008.dll" type="Microsoft.VisualStudio.TestTools.TestTypes.Unit.UnitTestElement, Microsoft.VisualStudio.QualityTools.Tips.UnitTest.ObjectModel,   PublicKeyToken=b03f5f7f11d50a3a" /> 
     82    </TestLinks> 
     83  </TestList> 
    7884  <TestList name="American Option" id="da9a808d-acf7-4baf-9dfb-e4b091d30fef" parentListId="8c43106b-9dc1-4907-a29f-aa66a61bf5b6"> 
    7985    <TestLinks> 
    80       <TestLink id="2bbf531b-68cb-84c0-a06f-6bedfb89f20f" name="testBaroneAdesiWhaleyValues" storage="test2008\bin\debug\test2008.dll" type="Microsoft.VisualStudio.TestTools.TestTypes.Unit.UnitTestElement, Microsoft.VisualStudio.QualityTools.Tips.UnitTest.ObjectModel,   PublicKeyToken=b03f5f7f11d50a3a" /> 
     86      <TestLink id="33c4d801-48ac-44c4-329e-d6fe84f5e034" name="testJuValues" storage="test2008\bin\debug\test2008.dll" type="Microsoft.VisualStudio.TestTools.TestTypes.Unit.UnitTestElement, Microsoft.VisualStudio.QualityTools.Tips.UnitTest.ObjectModel,   PublicKeyToken=b03f5f7f11d50a3a" /> 
    8187      <TestLink id="3be50179-04ee-902d-73d8-3c94e9305a67" name="testFdShoutGreeks" storage="test2008\bin\debug\test2008.dll" type="Microsoft.VisualStudio.TestTools.TestTypes.Unit.UnitTestElement, Microsoft.VisualStudio.QualityTools.Tips.UnitTest.ObjectModel,   PublicKeyToken=b03f5f7f11d50a3a" /> 
    8288      <TestLink id="783f96db-cce8-bf61-b280-f5010f7d7eb3" name="testFdValues" storage="test2008\bin\debug\test2008.dll" type="Microsoft.VisualStudio.TestTools.TestTypes.Unit.UnitTestElement, Microsoft.VisualStudio.QualityTools.Tips.UnitTest.ObjectModel,   PublicKeyToken=b03f5f7f11d50a3a" /> 
    8389      <TestLink id="1e42c90a-8e82-1369-9ed9-52b036b1c07e" name="testBjerksundStenslandValues" storage="test2008\bin\debug\test2008.dll" type="Microsoft.VisualStudio.TestTools.TestTypes.Unit.UnitTestElement, Microsoft.VisualStudio.QualityTools.Tips.UnitTest.ObjectModel,   PublicKeyToken=b03f5f7f11d50a3a" /> 
    8490      <TestLink id="9f4ff6f8-9938-8ecb-5608-7dbc928cbda3" name="testFdAmericanGreeks" storage="test2008\bin\debug\test2008.dll" type="Microsoft.VisualStudio.TestTools.TestTypes.Unit.UnitTestElement, Microsoft.VisualStudio.QualityTools.Tips.UnitTest.ObjectModel,   PublicKeyToken=b03f5f7f11d50a3a" /> 
    85       <TestLink id="33c4d801-48ac-44c4-329e-d6fe84f5e034" name="testJuValues" storage="test2008\bin\debug\test2008.dll" type="Microsoft.VisualStudio.TestTools.TestTypes.Unit.UnitTestElement, Microsoft.VisualStudio.QualityTools.Tips.UnitTest.ObjectModel,   PublicKeyToken=b03f5f7f11d50a3a" /> 
     91      <TestLink id="2bbf531b-68cb-84c0-a06f-6bedfb89f20f" name="testBaroneAdesiWhaleyValues" storage="test2008\bin\debug\test2008.dll" type="Microsoft.VisualStudio.TestTools.TestTypes.Unit.UnitTestElement, Microsoft.VisualStudio.QualityTools.Tips.UnitTest.ObjectModel,   PublicKeyToken=b03f5f7f11d50a3a" /> 
    8692    </TestLinks> 
    8793  </TestList> 

  • trunk/QLNet/QLNet/Math/CostFunction.cs

    r117 r204  
    3030        public abstract Vector values(Vector x); 
    3131 
    32         //! method to overload to compute the cost function values in x 
    33         //virtual Disposable<Array> values(const Array& x) const =0; 
    34  
    3532        //! method to overload to compute grad_f, the first derivative of 
    3633        //  the cost function with respect to x 
    3734        public virtual void gradient(Vector grad, Vector x) { 
    3835            double eps = finiteDifferenceEpsilon(), fp, fm; 
    39             Vector xx = (Vector)x.Clone(); 
     36            Vector xx = new Vector(x); 
    4037            for (int i=0; i<x.Count; i++) { 
    4138                xx[i] += eps; 

  • trunk/QLNet/QLNet/Math/Optimization/ArmijoLineSearch.cs

    r145 r204  
    11/* 
    22 Copyright (C) 2008 Toyin Akin (toyin_akin@hotmail.com) 
    3  *  
     3 Copyright (C) 2008 Siarhei Novik (snovik@gmail.com) 
     4  
    45 This file is part of QLNet Project http://www.qlnet.org 
    56 
    …… 
    2223using System.Text; 
    2324 
    24 namespace QLNet 
    25 { 
     25namespace QLNet { 
    2626 
    27         //! Armijo line search. 
    28 //    ! Let \f$ \alpha \f$ and \f$ \beta \f$ be 2 scalars in \f$ [0,1] 
    29 //        \f$.  Let \f$ x \f$ be the current value of the unknown, \f$ d 
    30 //        \f$ the search direction and \f$ t \f$ the step. Let \f$ f \f$ 
    31 //        be the function to minimize.  The line search stops when \f$ t 
    32 //        \f$ verifies 
    33 //        \f[ f(x + t \cdot d) - f(x) \leq -\alpha t f'(x+t \cdot d) \f] 
    34 //        and 
    35 //        \f[ f(x+\frac{t}{\beta} \cdot d) - f(x) > -\frac{\alpha}{\beta} 
    36 //            t f'(x+t \cdot d) \f] 
    37 // 
    38 //        (see Polak, Algorithms and consistent approximations, Optimization, 
    39 //        volume 124 of Applied Mathematical Sciences, Springer-Verlag, NY, 
    40 //        1997) 
    41 //     
    42         public class ArmijoLineSearch : LineSearch 
    43         { 
    44                 //! Default constructor 
    45                 public ArmijoLineSearch(double eps, double alpha) : this(eps, alpha, 0.65) 
    46                 { 
    47                 } 
    48                 public ArmijoLineSearch(double eps) : this(eps, 0.05, 0.65) 
    49                 { 
    50                 } 
    51                 public ArmijoLineSearch() : this(1e-8, 0.05, 0.65) 
    52                 { 
    53                 } 
    54                 public ArmijoLineSearch(double eps, double alpha, double beta) : base(eps) 
    55                 { 
    56                         alpha_ = alpha; 
    57                         beta_ = beta; 
    58                 } 
     27    //! Armijo line search. 
     28    //    ! Let \f$ \alpha \f$ and \f$ \beta \f$ be 2 scalars in \f$ [0,1] 
     29    //        \f$.  Let \f$ x \f$ be the current value of the unknown, \f$ d 
     30    //        \f$ the search direction and \f$ t \f$ the step. Let \f$ f \f$ 
     31    //        be the function to minimize.  The line search stops when \f$ t 
     32    //        \f$ verifies 
     33    //        \f[ f(x + t \cdot d) - f(x) \leq -\alpha t f'(x+t \cdot d) \f] 
     34    //        and 
     35    //        \f[ f(x+\frac{t}{\beta} \cdot d) - f(x) > -\frac{\alpha}{\beta} 
     36    //            t f'(x+t \cdot d) \f] 
     37    // 
     38    //        (see Polak, Algorithms and consistent approximations, Optimization, 
     39    //        volume 124 of Applied Mathematical Sciences, Springer-Verlag, NY, 
     40    //        1997) 
     41    //     
     42    public class ArmijoLineSearch : LineSearch { 
     43        //! Default constructor 
     44        public ArmijoLineSearch(double eps, double alpha) : this(eps, alpha, 0.65) { } 
     45        public ArmijoLineSearch(double eps) : this(eps, 0.05, 0.65) { } 
     46        public ArmijoLineSearch() : this(1e-8, 0.05, 0.65) { } 
     47        public ArmijoLineSearch(double eps, double alpha, double beta) 
     48            : base(eps) { 
     49            alpha_ = alpha; 
     50            beta_ = beta; 
     51        } 
    5952 
    60                 //! Perform line search 
    61                 public override double value(Problem P, ref EndCriteria.Type ecType, EndCriteria endCriteria, double t_ini) 
    62                 { 
    63                         //OptimizationMethod& method = P.method(); 
    64                         Constraint constraint = P.constraint(); 
    65                         succeed_=true; 
    66                         bool maxIter = false; 
    67                         double qtold; 
    68                         double t = t_ini; 
    69                         int loopNumber = 0; 
    70          
    71                         double q0 = P.functionValue(); 
    72                         double qp0 = P.gradientNormValue(); 
    73          
    74                         qt_ = q0; 
    75                         qpt_ = (gradient_.Count == 0) ? qp0 : -Utils.DotProduct(gradient_,searchDirection_); 
    76          
    77                         // Initialize gradient 
    78                         gradient_ = new Vector(P.currentValue().Count); 
    79                         // Compute new point 
    80                         xtd_ = P.currentValue(); 
    81                         t = update(xtd_, searchDirection_, t, constraint); 
    82                         // Compute function value at the new point 
    83                         qt_ = P.value (xtd_); 
    84          
    85                         // Enter in the loop if the criterion is not satisfied 
    86                         if ((qt_-q0) > -alpha_ *t *qpt_) 
    87                         { 
    88                                 do 
    89                                 { 
    90                                         loopNumber++; 
    91                                         // Decrease step 
    92                                         t *= beta_; 
    93                                         // Store old value of the function 
    94                                         qtold = qt_; 
    95                                         // New point value 
    96                                         xtd_ = P.currentValue(); 
    97                                         t = update(xtd_, searchDirection_, t, constraint); 
    98          
    99                                         // Compute function value at the new point 
    100                                         qt_ = P.value (xtd_); 
    101                                         P.gradient (gradient_, xtd_); 
    102                                         // and it squared norm 
    103                                         maxIter = endCriteria.checkMaxIterations(loopNumber, ref ecType); 
    104                                 } while ((((qt_ - q0) > (-alpha_ * t * qpt_)) || ((qtold - q0) <= (-alpha_ * t * qpt_ / beta_))) && (!maxIter)); 
    105                         } 
    106          
    107                         if (maxIter) 
    108                                 succeed_ = false; 
    109          
    110                         // Compute new gradient 
    111                         P.gradient(gradient_, xtd_); 
    112                         // and it squared norm 
     53        //! Perform line search 
     54        public override double value(Problem P, ref EndCriteria.Type ecType, EndCriteria endCriteria, double t_ini) { 
     55            //OptimizationMethod& method = P.method(); 
     56            Constraint constraint = P.constraint(); 
     57            succeed_ = true; 
     58            bool maxIter = false; 
     59            double qtold; 
     60            double t = t_ini; 
     61            int loopNumber = 0; 
     62 
     63            double q0 = P.functionValue(); 
     64            double qp0 = P.gradientNormValue(); 
     65 
     66            qt_ = q0; 
     67            qpt_ = (gradient_.Count == 0) ? qp0 : -Utils.DotProduct(gradient_, searchDirection_); 
     68 
     69            // Initialize gradient 
     70            gradient_ = new Vector(P.currentValue().Count); 
     71            // Compute new point 
     72            xtd_ = (Vector)P.currentValue().Clone(); 
     73            t = update(ref xtd_, searchDirection_, t, constraint); 
     74            // Compute function value at the new point 
     75            qt_ = P.value(xtd_); 
     76 
     77            // Enter in the loop if the criterion is not satisfied 
     78            if ((qt_ - q0) > -alpha_ * t * qpt_) { 
     79                do { 
     80                    loopNumber++; 
     81                    // Decrease step 
     82                    t *= beta_; 
     83                    // Store old value of the function 
     84                    qtold = qt_; 
     85                    // New point value 
     86                    xtd_ = P.currentValue(); 
     87                    t = update(ref xtd_, searchDirection_, t, constraint); 
     88 
     89                    // Compute function value at the new point 
     90                    qt_ = P.value(xtd_); 
     91                    P.gradient(gradient_, xtd_); 
     92                    // and it squared norm 
     93                    maxIter = endCriteria.checkMaxIterations(loopNumber, ref ecType); 
     94                } while ((((qt_ - q0) > (-alpha_ * t * qpt_)) || ((qtold - q0) <= (-alpha_ * t * qpt_ / beta_))) && (!maxIter)); 
     95            } 
     96 
     97            if (maxIter) 
     98                succeed_ = false; 
     99 
     100            // Compute new gradient 
     101            P.gradient(gradient_, xtd_); 
     102            // and it squared norm 
    113103            qpt_ = Utils.DotProduct(gradient_, gradient_); 
    114          
    115                        // Return new step value 
    116                        return t; 
    117                } 
    118                private double alpha_; 
    119                private double beta_; 
    120        } 
     104 
     105            // Return new step value 
     106            return t; 
     107        } 
     108        private double alpha_; 
     109        private double beta_; 
     110    } 
    121111} 

  • trunk/QLNet/QLNet/Math/Optimization/LineSearch.cs

    r145 r204  
    11/* 
    22 Copyright (C) 2008 Toyin Akin (toyin_akin@hotmail.com) 
     3 Copyright (C) 2008 Siarhei Novik (snovik@gmail.com) 
    34 *  
    45 This file is part of QLNet Project http://www.qlnet.org 
    …… 
    2223using System.Text; 
    2324 
    24 namespace QLNet 
    25 { 
     25namespace QLNet { 
    2626 
    2727    //! Base class for line search 
    28     public abstract class LineSearch 
    29     { 
     28    public abstract class LineSearch { 
    3029        //! current values of the search direction 
    3130        protected Vector searchDirection_; 
    3231        //! new x and its gradient 
    3332        protected Vector xtd_; 
    34         protected Vector gradient_; 
     33        protected Vector gradient_ = new Vector(); 
    3534        //! cost function value and gradient norm corresponding to xtd_ 
    3635        protected double qt_; 
    …… 
    4039 
    4140        //! Default constructor 
    42         public LineSearch() 
    43             : this(0.0) 
    44         { 
    45         } 
    46         public LineSearch(double UnnamedParameter1) 
    47         { 
     41        public LineSearch() : this(0.0) { } 
     42        public LineSearch(double UnnamedParameter1) { 
    4843            qt_ = 0.0; 
    4944            qpt_ = 0.0; 
    …… 
    5247 
    5348        //! return last x value 
    54         public Vector lastX() 
    55         { 
    56             return xtd_; 
    57         } 
     49        public Vector lastX() { return xtd_; } 
    5850        //! return last cost function value 
    59         public double lastFunctionValue() 
    60         { 
    61             return qt_; 
    62         } 
     51        public double lastFunctionValue() { return qt_; } 
    6352        //! return last gradient 
    64         public Vector lastGradient() 
    65         { 
    66             return gradient_; 
    67         } 
     53        public Vector lastGradient() { return gradient_; } 
    6854        //! return square norm of last gradient 
    69         public double lastGradientNorm2() 
    70         { 
    71             return qpt_; 
    72         } 
     55        public double lastGradientNorm2() { return qpt_; } 
    7356 
    74         public bool succeed() 
    75         { 
    76             return succeed_; 
    77         } 
     57        public bool succeed() { return succeed_; } 
    7858 
    7959        //! Perform line search 
    8060        public abstract double value(Problem P, ref EndCriteria.Type ecType, EndCriteria NamelessParameter3, double t_ini); // initial value of line-search step 
    81         public double update(Vector data, Vector direction, double beta, Constraint constraint) 
    82         { 
     61        public double update(ref Vector data, Vector direction, double beta, Constraint constraint) { 
    8362 
    8463            double diff = beta; 
    …… 
    8665            bool valid = constraint.test(newParams); 
    8766            int icount = 0; 
    88             while (!valid) 
    89             { 
     67            while (!valid) { 
    9068                if (icount > 200) 
    9169                    throw new ApplicationException("can't update linesearch"); 
    …… 
    10078 
    10179        //! current value of the search direction 
    102         public Vector searchDirection 
    103         { 
     80        public Vector searchDirection { 
    10481            get { return searchDirection_; } 
    10582            set { searchDirection_ = value; } 

  • trunk/QLNet/QLNet/Math/Optimization/levenbergmarquardt.cs

    r202 r204  
    3838 
    3939        //public LevenbergMarquardt(double epsfcn = 1.0e-8, double xtol = 1.0e-8, double gtol = 1.0e-8) { 
     40        public LevenbergMarquardt() : this(1.0e-8, 1.0e-8, 1.0e-8) { } 
    4041        public LevenbergMarquardt(double epsfcn, double xtol, double gtol) { 
    4142            info_ = 0; 
    …… 
    7273            // call lmdif to minimize the sum of the squares of m functions 
    7374            // in n variables by the Levenberg-Marquardt algorithm. 
    74             //MINPACK::LmdifCostFunction lmdifCostFunction =   
    75             //    boost::bind(&LevenbergMarquardt::fcn, this, _1, _2, _3, _4, _5); 
    76             //MINPACK::lmdif(m, n, xx.get(), fvec.get(), 
    77             //                         static_cast<double>(endCriteria.functionEpsilon()), 
    78             //                         static_cast<double>(xtol_), 
    79             //                         static_cast<double>(gtol_), 
    80             //                         static_cast<int>(endCriteria.maxIterations()), 
    81             //                         static_cast<double>(epsfcn_), 
    82             //                         diag.get(), mode, factor, 
    83             //                         nprint, &info, &nfev, fjac.get(), 
    84             //                         ldfjac, ipvt.get(), qtf.get(), 
    85             //                         wa1.get(), wa2.get(), wa3.get(), wa4.get(), 
    86             //                         lmdifCostFunction); 
    87             //info_ = info; 
    88             //// check requirements & endCriteria evaluation 
    89             //QL_REQUIRE(info != 0, "MINPACK: improper input parameters"); 
    90             ////QL_REQUIRE(info != 6, "MINPACK: ftol is too small. no further " 
    91             ////                               "reduction in the sum of squares " 
    92             ////                               "is possible."); 
    93             //if (info != 6) ecType = EndCriteria.Type.StationaryFunctionValue; 
    94             ////QL_REQUIRE(info != 5, "MINPACK: number of calls to fcn has " 
    95             ////                               "reached or exceeded maxfev."); 
    96             //endCriteria.checkMaxIterations(nfev, ecType); 
    97             //QL_REQUIRE(info != 7, "MINPACK: xtol is too small. no further " 
    98             //                               "improvement in the approximate " 
    99             //                               "solution x is possible."); 
    100             //QL_REQUIRE(info != 8, "MINPACK: gtol is too small. fvec is " 
    101             //                               "orthogonal to the columns of the " 
    102             //                               "jacobian to machine precision."); 
    103             //// set problem 
    104             //std::copy(xx.get(), xx.get()+n, x_.begin()); 
    105             //P.setCurrentValue(x_); 
     75            MINPACK.lmdif(m, n, xx, ref fvec, 
     76                                     endCriteria.functionEpsilon(), 
     77                                     xtol_, 
     78                                     gtol_, 
     79                                     endCriteria.maxIterations(), 
     80                                     epsfcn_, 
     81                                     diag, mode, factor, 
     82                                     nprint, ref info, ref nfev, ref fjac, 
     83                                     ldfjac, ref ipvt, ref qtf, 
     84                                     wa1, wa2, wa3, wa4, 
     85                                     fcn); 
     86            info_ = info; 
     87            // check requirements & endCriteria evaluation 
     88            if(info == 0) throw new ApplicationException("MINPACK: improper input parameters"); 
     89            //if(info == 6) throw new ApplicationException("MINPACK: ftol is too small. no further " + 
     90            //                                             "reduction in the sum of squares is possible."); 
     91 
     92            if (info != 6) ecType = EndCriteria.Type.StationaryFunctionValue; 
     93            //QL_REQUIRE(info != 5, "MINPACK: number of calls to fcn has reached or exceeded maxfev."); 
     94            endCriteria.checkMaxIterations(nfev, ref ecType); 
     95            if(info == 7) throw new ApplicationException("MINPACK: xtol is too small. no further " + 
     96                                           "improvement in the approximate " + 
     97                                           "solution x is possible."); 
     98            if(info == 8) throw new ApplicationException("MINPACK: gtol is too small. fvec is " + 
     99                                           "orthogonal to the columns of the " + 
     100                                           "jacobian to machine precision."); 
     101            // set problem 
     102            x_ = new Vector(xx.GetRange(0, n)); 
     103            P.setCurrentValue(x_); 
    106104 
    107105            return ecType; 
    108106        } 
    109107 
    110         public void fcn(int m, int n, Vector x, out Vector fvec, int iflag) { 
     108        public Vector fcn(int m, int n, Vector x, int iflag) { 
    111109            Vector xt = new Vector(x); 
    112  
     110            Vector fvec; 
    113111            // constraint handling needs some improvement in the future: 
    114112            // starting point should not be close to a constraint violation 
    …… 
    118116                fvec = new Vector(initCostValues_); 
    119117            } 
     118            return fvec; 
    120119        } 
    121120    } 

  • trunk/QLNet/QLNet/Math/Optimization/lmdif.cs

    r202 r204  
    8383 
    8484namespace QLNet { 
    85     public class MINPACK { 
     85    public static class MINPACK { 
     86        /* resolution of arithmetic */ 
     87        const double MACHEP = 1.2e-16; 
     88        /* smallest nonzero number */ 
     89        const double DWARF = 1.0e-38; 
     90 
    8691        /* 
    8792        *     ********** 
    …… 
    259264        *   minpack-supplied ... dpmpar,enorm,fdjac2,lmpar,qrfac 
    260265        * 
    261         *   fortran-supplied ... dabs,dmax1,dmin1,dsqrt,mod 
    262         * 
    263266        *     argonne national laboratory. minpack project. march 1980. 
    264267        *     burton s. garbow, kenneth e. hillstrom, jorge j. more 
    …… 
    266269        *     ********** 
    267270        */ 
    268         void lmdif(int m, int n, Vector x, Vector fvec, double ftol, 
    269                    double xtol,double gtol,int maxfev,double epsfcn, 
    270                    Vector diag, int mode, double factor, 
    271                    int nprint, out int info, out int nfev, Vector fjac, 
    272                    int ldfjac, List<int> ipvt, Vector qtf, 
    273                    Vector wa1, Vector wa2, Vector wa3, Vector wa4, 
    274                    Action<int, int, Vector, Vector, int> fcn) { 
    275  
    276             int i,iflag,ij,jj,iter,j,l; 
    277             double actred,delta=0,dirder,fnorm,fnorm1,gnorm; 
    278             double par,pnorm,prered,ratio; 
    279             double sum,temp,temp1,temp2,temp3,xnorm=0; 
    280  
    281             double one = 1.0; 
    282             double p1 = 0.1; 
    283             double p5 = 0.5; 
    284             double p25 = 0.25; 
    285             double p75 = 0.75; 
    286             double p0001 = 1.0e-4; 
    287             double zero = 0.0; 
     271        public static void lmdif(int m, int n, Vector x, ref Vector fvec, double ftol, 
     272                                 double xtol, double gtol, int maxfev, double epsfcn, 
     273                                 Vector diag, int mode, double factor, 
     274                                 int nprint, ref int info, ref int nfev, ref Vector fjac, 
     275                                 int ldfjac, ref List<int> ipvt, ref Vector qtf, 
     276                                 Vector wa1, Vector wa2, Vector wa3, Vector wa4, 
     277                                 Func<int, int, Vector, int, Vector> fcn) { 
     278 
     279            int i, iflag, ij, jj, iter, j, l; 
     280            double actred, delta = 0, dirder, fnorm, fnorm1, gnorm; 
     281            double par, pnorm, prered, ratio; 
     282            double sum, temp, temp1, temp2, temp3, xnorm = 0; 
     283 
     284            const double one = 1.0; 
     285            const double p1 = 0.1; 
     286            const double p5 = 0.5; 
     287            const double p25 = 0.25; 
     288            const double p75 = 0.75; 
     289            const double p0001 = 1.0e-4; 
     290            const double zero = 0.0; 
    288291 
    289292            info = 0; 
    …… 
    294297            *     check the input parameters for errors. 
    295298            */ 
    296             if( (n <= 0) || (m < n) || (ldfjac < m) || (ftol < zero) 
     299            if ((n <= 0) || (m < n) || (ldfjac < m) || (ftol < zero) 
    297300                || (xtol < zero) || (gtol < zero) || (maxfev <= 0) 
    298                 || (factor <= zero) ) 
     301                || (factor <= zero)) 
    299302                goto L300; 
    300303 
    301             if( mode == 2 ) 
    302                 { /* scaling by diag[] */ 
    303                 for( j=0; j<n; j++ ) 
    304                     { 
    305                     if( diag[j] <= 0.0 ) 
     304            if (mode == 2) { /* scaling by diag[] */ 
     305                for (j = 0; j < n; j++) { 
     306                    if (diag[j] <= 0.0) 
    306307                        goto L300; 
    307                     } 
    308                 } 
     308                } 
     309            } 
    309310            /* 
    310311            *     evaluate the function at the starting point 
    …… 
    312313            */ 
    313314            iflag = 1; 
    314             fcn(m,n,x, ref fvec,iflag); 
     315            fvec = fcn(m, n, x, iflag); 
    315316            nfev = 1; 
    316             if(iflag < 0) 
     317            if (iflag < 0) 
    317318                goto L300; 
    318             fnorm = enorm(m,fvec); 
     319            fnorm = enorm(m, fvec); 
    319320            /* 
    320321            *     initialize levenberg-marquardt parameter and iteration counter. 
    …… 
    322323            par = zero; 
    323324            iter = 1; 
    324             /* 
    325             *     beginning of the outer loop. 
    326             */ 
     325        /* 
     326        *     beginning of the outer loop. 
     327        */ 
    327328 
    328329            L30: 
    …… 
    332333            */ 
    333334            iflag = 2; 
    334             fdjac2(m,n,x,fvec,fjac,ldfjac,iflag,epsfcn,wa4, fcn); 
     335            fdjac2(m, n, x, fvec, fjac, ldfjac, iflag, epsfcn, ref wa4, fcn); 
    335336            nfev += n; 
    336             if(iflag < 0) 
     337            if (iflag < 0) 
    337338                goto L300; 
    338339            /* 
    339340            *    if requested, call fcn to enable printing of iterates. 
    340341            */ 
    341             if( nprint > 0 ) 
    342                 { 
     342            if (nprint > 0) { 
    343343                iflag = 0; 
    344                 if(mod(iter-1,nprint) == 0) 
    345                     { 
    346                     fcn(m,n,x,fvec,iflag); 
    347                     if(iflag < 0) 
     344                if (mod(iter - 1, nprint) == 0) { 
     345                    fvec = fcn(m, n, x, iflag); 
     346                    if (iflag < 0) 
    348347                        goto L300; 
    349                     } 
    350                 } 
     348                } 
     349            } 
    351350            /* 
    352351            *    compute the qr factorization of the jacobian. 
    353352            */ 
    354             qrfac(m,n,fjac,ldfjac,1,ipvt,n,wa1,wa2,wa3); 
     353            qrfac(m, n, fjac, ldfjac, 1, ref ipvt, n, ref wa1, ref wa2, wa3); 
    355354            /* 
    356355            *    on the first iteration and if mode is 1, scale according 
    357356            *    to the norms of the columns of the initial jacobian. 
    358357            */ 
    359             if(iter == 1) 
    360                 { 
    361                 if(mode != 2) 
    362                     { 
    363                     for( j=0; j<n; j++ ) 
    364                         { 
     358            if (iter == 1) { 
     359                if (mode != 2) { 
     360                    for (j = 0; j < n; j++) { 
    365361                        diag[j] = wa2[j]; 
    366                         if( wa2[j] == zero ) 
     362                        if (wa2[j] == zero) 
    367363                            diag[j] = one; 
    368                         } 
    369364                    } 
    370  
    371             /* 
    372             *    on the first iteration, calculate the norm of the scaled x 
    373             *    and initialize the step bound delta. 
    374             */ 
    375                 for( j=0; j<n; j++ ) 
     365                } 
     366 
     367                /* 
     368                *    on the first iteration, calculate the norm of the scaled x 
     369                *    and initialize the step bound delta. 
     370                */ 
     371                for (j = 0; j < n; j++) 
    376372                    wa3[j] = diag[j] * x[j]; 
    377373 
    378                 xnorm = enorm(n,wa3); 
    379                 delta = factor*xnorm; 
    380                 if(delta == zero) 
     374                xnorm = enorm(n, wa3); 
     375                delta = factor * xnorm; 
     376                if (delta == zero) 
    381377                    delta = factor; 
    382                 } 
     378            } 
    383379 
    384380            /* 
    …… 
    386382            *    qtf. 
    387383            */ 
    388             for( i=0; i<m; i++ ) 
     384            for (i = 0; i < m; i++) 
    389385                wa4[i] = fvec[i]; 
    390386            jj = 0; 
    391             for( j=0; j<n; j++ ) 
    392                 { 
     387            for (j = 0; j < n; j++) { 
    393388                temp3 = fjac[jj]; 
    394                 if(temp3 != zero) 
    395                     { 
     389                if (temp3 != zero) { 
    396390                    sum = zero; 
    397391                    ij = jj; 
    398                     for( i=j; i<m; i++ ) 
    399                         { 
     392                    for (i = j; i < m; i++) { 
    400393                        sum += fjac[ij] * wa4[i]; 
    401394                        ij += 1;    /* fjac[i+m*j] */ 
    402                         } 
     395                    } 
    403396                    temp = -sum / temp3; 
    404397                    ij = jj; 
    405                     for( i=j; i<m; i++ ) 
    406                         { 
     398                    for (i = j; i < m; i++) { 
    407399                        wa4[i] += fjac[ij] * temp; 
    408400                        ij += 1;    /* fjac[i+m*j] */ 
    409                         } 
    410401                    } 
     402                } 
    411403                fjac[jj] = wa1[j]; 
    412                 jj += m+1;  /* fjac[j+m*j] */ 
     404                jj += m + 1;  /* fjac[j+m*j] */ 
    413405                qtf[j] = wa4[j]; 
    414                 } 
     406            } 
    415407 
    416408            /* 
    417409            *    compute the norm of the scaled gradient. 
    418410            */ 
    419              gnorm = zero; 
    420              if(fnorm != zero) 
    421                 { 
     411            gnorm = zero; 
     412            if (fnorm != zero) { 
    422413                jj = 0; 
    423                 for( j=0; j<n; j++ ) 
    424                     { 
     414                for (j = 0; j < n; j++) { 
    425415                    l = ipvt[j]; 
    426                     if(wa2[l] != zero) 
    427                         { 
     416                    if (wa2[l] != zero) { 
    428417                        sum = zero; 
    429418                        ij = jj; 
    430                         for( i=0; i<=j; i++ ) 
    431                             { 
    432                             sum += fjac[ij]*(qtf[i]/fnorm); 
     419                        for (i = 0; i <= j; i++) { 
     420                            sum += fjac[ij] * (qtf[i] / fnorm); 
    433421                            ij += 1; /* fjac[i+m*j] */ 
    434                             } 
    435                         gnorm = dmax1(gnorm,Math.Abs(sum/wa2[l])); 
    436422                        } 
     423                        gnorm = dmax1(gnorm, Math.Abs(sum / wa2[l])); 
     424                    } 
    437425                    jj += m; 
    438                     } 
    439                 } 
     426                } 
     427            } 
    440428 
    441429            /* 
    442430            *    test for convergence of the gradient norm. 
    443431            */ 
    444             if(gnorm <= gtol) 
     432            if (gnorm <= gtol) 
    445433                info = 4; 
    446             if( info != 0) 
     434            if (info != 0) 
    447435                goto L300; 
    448436            /* 
    449437            *    rescale if necessary. 
    450438            */ 
    451              if(mode != 2) 
    452                 { 
    453                 for( j=0; j<n; j++ ) 
    454                     diag[j] = dmax1(diag[j],wa2[j]); 
    455                 } 
     439            if (mode != 2) { 
     440                for (j = 0; j < n; j++) 
     441                    diag[j] = dmax1(diag[j], wa2[j]); 
     442            } 
    456443 
    457444            /* 
    458445            *    beginning of the inner loop. 
    459446            */ 
    460             L200: 
     447        L200: 
    461448            /* 
    462449            *       determine the levenberg-marquardt parameter. 
    463450            */ 
    464             lmpar(n,fjac,ldfjac,ipvt,diag,qtf,delta,&par,wa1,wa2,wa3,wa4); 
     451            lmpar(n, fjac, ldfjac, ipvt, diag, qtf, delta, par, wa1, wa2, wa3, wa4); 
    465452            /* 
    466453            *       store the direction p and x + p. calculate the norm of p. 
    467454            */ 
    468             for( j=0; j<n; j++ ) 
    469                 { 
    470                    wa1[j] = -wa1[j]; 
    471                    wa2[j] = x[j] + wa1[j]; 
    472                    wa3[j] = diag[j]*wa1[j]; 
    473                 } 
    474             pnorm = enorm(n,wa3); 
     455            for (j = 0; j < n; j++) { 
     456                wa1[j] = -wa1[j]; 
     457                wa2[j] = x[j] + wa1[j]; 
     458                wa3[j] = diag[j] * wa1[j]; 
     459            } 
     460            pnorm = enorm(n, wa3); 
    475461            /* 
    476462            *       on the first iteration, adjust the initial step bound. 
    477463            */ 
    478             if(iter == 1) 
    479                 delta = dmin1(delta,pnorm); 
     464            if (iter == 1) 
     465                delta = dmin1(delta, pnorm); 
    480466            /* 
    481467            *       evaluate the function at x + p and calculate its norm. 
    482468            */ 
    483469            iflag = 1; 
    484             fcn(m,n,wa2,wa4,&iflag); 
     470            wa4 = fcn(m, n, wa2, iflag); 
    485471            nfev += 1; 
    486             if(iflag < 0) 
     472            if (iflag < 0) 
    487473                goto L300; 
    488             fnorm1 = enorm(m,wa4); 
    489             #if BUG 
    490             printf( "pnorm %.10e  fnorm1 %.10e\n", pnorm, fnorm1 ); 
    491             #endif 
     474          &nb