Assembla home | Assembla project page
 

Changeset 146

Show
Ignore:
Timestamp:
05/05/08 10:31:44 (2 months ago)
Author:
ToyinA
Message:

Bug fix for Simplex

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

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

    r145 r146  
    4141        } 
    4242 
    43         //! Multi-dimensional simplex class 
    44         public class Simplex : OptimizationMethod 
    45         { 
    46             //! Constructor taking as input the characteristic length  
    47             public Simplex(double lambda) 
    48             { 
    49                 lambda_ = lambda; 
     43    } 
     44 
     45    //! Multi-dimensional simplex class 
     46    public class Simplex : OptimizationMethod 
     47    { 
     48        //! Constructor taking as input the characteristic length  
     49        public Simplex(double lambda) 
     50        { 
     51            lambda_ = lambda; 
     52        } 
     53        public override EndCriteria.Type minimize(Problem P, EndCriteria endCriteria) 
     54        { 
     55            // set up of the problem 
     56            //double ftol = endCriteria.functionEpsilon();    // end criteria on f(x) (see Numerical Recipes in C++, p.410) 
     57            double xtol = endCriteria.rootEpsilon(); // end criteria on x (see GSL v. 1.9, http://www.gnu.org/software/gsl/) 
     58            int maxStationaryStateIterations_ = endCriteria.maxStationaryStateIterations(); 
     59            EndCriteria.Type ecType = EndCriteria.Type.None; 
     60            P.reset(); 
     61            Vector x_ = P.currentValue(); 
     62            int iterationNumber_ = 0; 
     63 
     64            // Initialize vertices of the simplex 
     65            bool end = false; 
     66            int n = x_.Count; 
     67            vertices_ = new InitializedList<Vector>(n + 1, x_); 
     68            for (int i = 0; i < n; i++) 
     69            { 
     70                Vector direction = new Vector(n, 0.0); 
     71                direction[i] = 1.0; 
     72                P.constraint().update(vertices_[i + 1], direction, lambda_); 
    5073            } 
    51             public override EndCriteria.Type minimize(Problem P, EndCriteria endCriteria) 
    52             { 
    53                 // set up of the problem 
    54                 //double ftol = endCriteria.functionEpsilon();    // end criteria on f(x) (see Numerical Recipes in C++, p.410) 
    55                 double xtol = endCriteria.rootEpsilon(); // end criteria on x (see GSL v. 1.9, http://www.gnu.org/software/gsl/) 
    56                 int maxStationaryStateIterations_ = endCriteria.maxStationaryStateIterations(); 
    57                 EndCriteria.Type ecType = EndCriteria.Type.None; 
    58                 P.reset(); 
    59                 Vector x_ = P.currentValue(); 
    60                 int iterationNumber_ = 0; 
    61  
    62                 // Initialize vertices of the simplex 
    63                 bool end = false; 
    64                 int n = x_.Count; 
    65                 vertices_ = new InitializedList<Vector>(n + 1, x_); 
    66                 for (int i = 0; i < n; i++) 
    67                 { 
    68                     Vector direction = new Vector(n, 0.0); 
    69                     direction[i] = 1.0; 
    70                     P.constraint().update(vertices_[i + 1], direction, lambda_); 
    71                 } 
    72                 // Initialize function values at the vertices of the simplex 
    73                 values_ = new Vector(n + 1, 0.0); 
     74            // Initialize function values at the vertices of the simplex 
     75            values_ = new Vector(n + 1, 0.0); 
     76            for (int i = 0; i <= n; i++) 
     77                values_[i] = P.value(vertices_[i]); 
     78            // Loop looking for minimum 
     79            do 
     80            { 
     81                sum_ = new Vector(n, 0.0); 
    7482                for (int i = 0; i <= n; i++) 
    75                     values_[i] = P.value(vertices_[i]); 
    76                 // Loop looking for minimum 
    77                 do 
    78                 { 
    79                     sum_ = new Vector(n, 0.0); 
    80                     for (int i = 0; i <= n; i++) 
    81                         sum_ += vertices_[i]; 
    82                     // Determine the best (iLowest), worst (iHighest) 
    83                     // and 2nd worst (iNextHighest) vertices 
    84                     int iLowest = 0; 
    85                     int iHighest; 
    86                     int iNextHighest; 
    87                     if (values_[0] < values_[1]) 
     83                    sum_ += vertices_[i]; 
     84                // Determine the best (iLowest), worst (iHighest) 
     85                // and 2nd worst (iNextHighest) vertices 
     86                int iLowest = 0; 
     87                int iHighest; 
     88                int iNextHighest; 
     89                if (values_[0] < values_[1]) 
     90                { 
     91                    iHighest = 1; 
     92                    iNextHighest = 0; 
     93                } 
     94                else 
     95                { 
     96                    iHighest = 0; 
     97                    iNextHighest = 1; 
     98                } 
     99                for (int i = 1; i <= n; i++) 
     100                { 
     101                    if (values_[i] > values_[iHighest]) 
    88102                    { 
    89                         iHighest = 1; 
    90                         iNextHighest = 0; 
     103                        iNextHighest = iHighest; 
     104                        iHighest = i; 
    91105                    } 
    92106                    else 
    93107                    { 
    94                         iHighest = 0; 
    95                         iNextHighest = 1; 
     108                        if ((values_[i] > values_[iNextHighest]) && i != iHighest) 
     109                            iNextHighest = i; 
    96110                    } 
    97                     for (int i = 1; i <= n; i++) 
     111                    if (values_[i] < values_[iLowest]) 
     112                        iLowest = i; 
     113                } 
     114                // Now compute accuracy, update iteration number and check end criteria 
     115                //// Numerical Recipes exit strategy on fx (see NR in C++, p.410) 
     116                //double low = values_[iLowest]; 
     117                //double high = values_[iHighest]; 
     118                //double rtol = 2.0*std::fabs(high - low)/ 
     119                //    (std::fabs(high) + std::fabs(low) + QL_EPSILON); 
     120                //++iterationNumber_; 
     121                //if (rtol < ftol || 
     122                //    endCriteria.checkMaxIterations(iterationNumber_, ecType)) { 
     123                // GSL exit strategy on x (see GSL v. 1.9, http://www.gnu.org/software/gsl 
     124                double simplexSize = Utils.computeSimplexSize(vertices_); 
     125                ++iterationNumber_; 
     126                if (simplexSize < xtol || endCriteria.checkMaxIterations(iterationNumber_, ref ecType)) 
     127                { 
     128                    endCriteria.checkStationaryPoint(0.0, 0.0, ref maxStationaryStateIterations_, ref ecType); 
     129                    endCriteria.checkMaxIterations(iterationNumber_, ref ecType); 
     130                    x_ = vertices_[iLowest]; 
     131                    double low = values_[iLowest]; 
     132                    P.setFunctionValue(low); 
     133                    P.setCurrentValue(x_); 
     134                    return ecType; 
     135                } 
     136                // If end criteria is not met, continue 
     137                double factor = -1.0; 
     138                double vTry = extrapolate(ref P, iHighest, ref factor); 
     139                if ((vTry <= values_[iLowest]) && (factor == -1.0)) 
     140                { 
     141                    factor = 2.0; 
     142                    extrapolate(ref P, iHighest, ref factor); 
     143                } 
     144                else 
     145                { 
     146                    if (vTry >= values_[iNextHighest]) 
    98147                    { 
    99                         if (values_[i] > values_[iHighest]) 
     148                        double vSave = values_[iHighest]; 
     149                        factor = 0.5; 
     150                        vTry = extrapolate(ref P, iHighest, ref factor); 
     151                        if (vTry >= vSave) 
    100152                        { 
    101                             iNextHighest = iHighest; 
    102                             iHighest = i; 
    103                         } 
    104                         else 
    105                         { 
    106                             if ((values_[i] > values_[iNextHighest]) && i != iHighest) 
    107                                 iNextHighest = i; 
    108                         } 
    109                         if (values_[i] < values_[iLowest]) 
    110                             iLowest = i; 
    111                     } 
    112                     // Now compute accuracy, update iteration number and check end criteria 
    113                     //// Numerical Recipes exit strategy on fx (see NR in C++, p.410) 
    114                     //double low = values_[iLowest]; 
    115                     //double high = values_[iHighest]; 
    116                     //double rtol = 2.0*std::fabs(high - low)/ 
    117                     //    (std::fabs(high) + std::fabs(low) + QL_EPSILON); 
    118                     //++iterationNumber_; 
    119                     //if (rtol < ftol || 
    120                     //    endCriteria.checkMaxIterations(iterationNumber_, ecType)) { 
    121                     // GSL exit strategy on x (see GSL v. 1.9, http://www.gnu.org/software/gsl 
    122                     double simplexSize = Utils.computeSimplexSize(vertices_); 
    123                     ++iterationNumber_; 
    124                     if (simplexSize < xtol || endCriteria.checkMaxIterations(iterationNumber_, ref ecType)) 
    125                     { 
    126                         endCriteria.checkStationaryPoint(0.0, 0.0, ref maxStationaryStateIterations_, ref ecType); 
    127                         endCriteria.checkMaxIterations(iterationNumber_, ref ecType); 
    128                         x_ = vertices_[iLowest]; 
    129                         double low = values_[iLowest]; 
    130                         P.setFunctionValue(low); 
    131                         P.setCurrentValue(x_); 
    132                         return ecType; 
    133                     } 
    134                     // If end criteria is not met, continue 
    135                     double factor = -1.0; 
    136                     double vTry = extrapolate(ref P, iHighest, ref factor); 
    137                     if ((vTry <= values_[iLowest]) && (factor == -1.0)) 
    138                     { 
    139                         factor = 2.0; 
    140                         extrapolate(ref P, iHighest, ref factor); 
    141                     } 
    142                     else 
    143                     { 
    144                         if (vTry >= values_[iNextHighest]) 
    145                         { 
    146                             double vSave = values_[iHighest]; 
    147                             factor = 0.5; 
    148                             vTry = extrapolate(ref P, iHighest, ref factor); 
    149                             if (vTry >= vSave) 
     153                            for (int i = 0; i <= n; i++) 
    150154                            { 
    151                                 for (int i = 0; i <= n; i++) 
     155                                if (i != iLowest) 
    152156                                { 
    153                                     if (i != iLowest) 
    154                                     { 
    155 //#if QL_ARRAY_EXPRESSIONS 
    156                                                                                 vertices_[i] = 0.5*(vertices_[i] + vertices_[iLowest]); 
    157 //#else 
    158 //                                        vertices_[i] += vertices_[iLowest]; 
    159 //                                        vertices_[i] *= 0.5; 
    160 //#endif 
    161 //                                        values_[i] = P.value(vertices_[i]); 
    162                                     } 
     157                                    //#if QL_ARRAY_EXPRESSIONS 
     158                                    vertices_[i] = 0.5 * (vertices_[i] + vertices_[iLowest]); 
     159                                    //#else 
     160                                    //                                        vertices_[i] += vertices_[iLowest]; 
     161                                    //                                        vertices_[i] *= 0.5; 
     162                                    //#endif 
     163                                    //                                        values_[i] = P.value(vertices_[i]); 
    163164                                } 
    164165                            } 
    165166                        } 
    166167                    } 
    167                 } while (end == false); 
    168                 throw new ApplicationException("optimization failed: unexpected behaviour"); 
     168                } 
     169            } while (end == false); 
     170            throw new ApplicationException("optimization failed: unexpected behaviour"); 
     171        } 
     172 
     173        private double extrapolate(ref Problem P, int iHighest, ref double factor) 
     174        { 
     175 
     176            Vector pTry; 
     177            do 
     178            { 
     179                int dimensions = values_.Count - 1; 
     180                double factor1 = (1.0 - factor) / dimensions; 
     181                double factor2 = factor1 - factor; 
     182                // #if QL_ARRAY_EXPRESSIONS 
     183                pTry = sum_ * factor1 - vertices_[iHighest] * factor2; 
     184                //#else 
     185                //                    // composite expressions fail to compile with gcc 3.4 on windows 
     186                //                    pTry = sum_ * factor1; 
     187                //                    pTry -= vertices_[iHighest] * factor2; 
     188                //#endif 
     189                factor *= 0.5; 
     190            } while (!P.constraint().test(pTry)); 
     191            factor *= 2.0; 
     192            double vTry = P.value(pTry); 
     193            if (vTry < values_[iHighest]) 
     194            { 
     195                values_[iHighest] = vTry; 
     196                //#if QL_ARRAY_EXPRESSIONS 
     197                sum_ += pTry - vertices_[iHighest]; 
     198                //#else 
     199                //                    sum_ += pTry; 
     200                //                    sum_ -= vertices_[iHighest]; 
     201                //#endif 
     202                vertices_[iHighest] = pTry; 
    169203            } 
    170  
    171             private double extrapolate(ref Problem P, int iHighest, ref double factor) 
    172             { 
    173  
    174                 Vector pTry; 
    175                 do 
    176                 { 
    177                     int dimensions = values_.Count - 1; 
    178                     double factor1 = (1.0 - factor) / dimensions; 
    179                     double factor2 = factor1 - factor; 
    180 // #if QL_ARRAY_EXPRESSIONS 
    181                                     pTry = sum_ *factor1 - vertices_[iHighest]*factor2; 
    182 //#else 
    183 //                    // composite expressions fail to compile with gcc 3.4 on windows 
    184 //                    pTry = sum_ * factor1; 
    185 //                    pTry -= vertices_[iHighest] * factor2; 
    186 //#endif 
    187                     factor *= 0.5; 
    188                 } while (!P.constraint().test(pTry)); 
    189                 factor *= 2.0; 
    190                 double vTry = P.value(pTry); 
    191                 if (vTry < values_[iHighest]) 
    192                 { 
    193                     values_[iHighest] = vTry; 
    194 //#if QL_ARRAY_EXPRESSIONS 
    195                                     sum_ += pTry - vertices_[iHighest]; 
    196 //#else 
    197 //                    sum_ += pTry; 
    198 //                    sum_ -= vertices_[iHighest]; 
    199 //#endif 
    200                     vertices_[iHighest] = pTry; 
    201                 } 
    202                 return vTry; 
    203  
    204             } 
    205             private double lambda_; 
    206             private InitializedList<Vector> vertices_; 
    207             private Vector values_; 
    208             private Vector sum_; 
    209         } 
     204            return vTry; 
     205 
     206        } 
     207        private double lambda_; 
     208        private InitializedList<Vector> vertices_; 
     209        private Vector values_; 
     210        private Vector sum_; 
    210211    } 
    211212}