Changeset 9291


Ignore:
Timestamp:
03/06/13 17:33:38 (7 years ago)
Author:
abeham
Message:

#1961: Removed dependency to ALGLIB in CMA-ES implementation (code is commented)

Location:
branches/CMAES/HeuristicLab.Algorithms.CMAEvolutionStrategy/3.3
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • branches/CMAES/HeuristicLab.Algorithms.CMAEvolutionStrategy/3.3/CMAOperators/CMAMutator.cs

    r9245 r9291  
    2020#endregion
    2121
    22 using System;
    2322using HeuristicLab.Common;
    2423using HeuristicLab.Core;
     
    2928using HeuristicLab.Parameters;
    3029using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
     30using HeuristicLab.Random;
     31using System;
    3132
    3233namespace HeuristicLab.Algorithms.CMAEvolutionStrategy {
     
    113114        RealVectorParameter.ActualValue = arx;
    114115      }
    115       alglib.hqrndstate state;
    116       alglib.hqrndseed(random.Next(), random.Next(), out state);
     116      var nd = new NormalDistributedRandom(random, 0, 1);
     117      //alglib.hqrndstate state;
     118      //alglib.hqrndseed(random.Next(), random.Next(), out state);
    117119
    118120      for (int i = 0; i < lambda; i++) {
     
    122124          for (int k = 0; k < arx[i].Length; k++) {
    123125            do {
    124               arx[i][k] = xmean[k] + sp.Sigma.Value * sp.D[k] * alglib.hqrndnormal(state);
     126              arx[i][k] = xmean[k] + sp.Sigma.Value * sp.D[k] * nd.NextDouble(); //alglib.hqrndnormal(state);
    125127              inRange = bounds[k % bounds.Rows, 0] <= arx[i][k] && arx[i][k] <= bounds[k % bounds.Rows, 1];
    126128              if (!inRange) tries++;
     
    138140            var artmp = new double[arx[0].Length];
    139141            for (int k = 0; k < arx[0].Length; ++k) {
    140               artmp[k] = sp.D[k] * alglib.hqrndnormal(state);
     142              artmp[k] = sp.D[k] * nd.NextDouble(); //alglib.hqrndnormal(state);
    141143            }
    142144
  • branches/CMAES/HeuristicLab.Algorithms.CMAEvolutionStrategy/3.3/CMAOperators/CMAUpdater.cs

    r9148 r9291  
    237237        DegenerateStateParameter.ActualValue = new BoolValue(false);
    238238      } else {
    239         var c = new double[sp.C.Rows, sp.C.Columns];
     239        /*var c = new double[sp.C.Rows, sp.C.Columns];
    240240        for (int i = 0; i < sp.C.Rows; i++)
    241241          for (int j = 0; j < sp.C.Columns; j++)
     
    243243
    244244        double[] eVal;
    245         double[,] eVec;
    246         var success = alglib.smatrixevd(c, N, 1, false, out eVal, out eVec);
     245        double[,] eVec;*/
     246
     247        // set B <- C
     248        for (int i = 0; i < N; ++i)
     249          for (int j = 0; j <= i; ++j)
     250            sp.B[i, j] = sp.B[j, i] = sp.C[i, j];
     251
     252        var success = Eigendecomposition(N, sp.B, sp.D); // alglib.smatrixevd(c, N, 1, false, out eVal, out eVec);
    247253        DegenerateStateParameter.ActualValue = new BoolValue(!success);
    248         sp.B = new DoubleMatrix(eVec);
    249 
    250         for (int i = 0; i < sp.D.Length; i++) {
     254        //sp.B = new DoubleMatrix(eVec);
     255
     256        /*for (int i = 0; i < sp.D.Length; i++) {
    251257          if (eVal[i] < 0) {
    252258            DegenerateStateParameter.ActualValue.Value = true;
    253259            sp.D[i] = 0;
    254260          } else sp.D[i] = Math.Sqrt(eVal[i]);
    255         }
     261        }*/
    256262
    257263        if (sp.D.Min() == 0.0) sp.AxisRatio.Value = double.PositiveInfinity;
     
    260266      return base.Apply();
    261267    }
     268
     269    private bool Eigendecomposition(int N, DoubleMatrix B, DoubleArray diagD) {
     270      bool result = true;
     271      // eigendecomposition
     272      var offdiag = new double[N];
     273      tred2(N, B, diagD, offdiag);
     274      tql2(N, diagD, offdiag, B);
     275
     276      // assign diagD to eigenvalue square roots
     277      for (int i = 0; i < N; ++i) {
     278        if (diagD[i] < 0) // numerical problem?
     279          result = false;
     280        else diagD[i] = Math.Sqrt(diagD[i]);
     281      }
     282
     283      return result;
     284    } // eigendecomposition
     285
     286
     287    // Symmetric Householder reduction to tridiagonal form, taken from JAMA package.
     288    private void tred2(int n, DoubleMatrix V, DoubleArray d, double[] e) {
     289
     290      //  This is derived from the Algol procedures tred2 by
     291      //  Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
     292      //  Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
     293      //  Fortran subroutine in EISPACK.
     294
     295      for (int j = 0; j < n; j++) {
     296        d[j] = V[n - 1, j];
     297      }
     298
     299      // Householder reduction to tridiagonal form.
     300
     301      for (int i = n - 1; i > 0; i--) {
     302
     303        // Scale to avoid under/overflow.
     304
     305        double scale = 0.0;
     306        double h = 0.0;
     307        for (int k = 0; k < i; k++) {
     308          scale = scale + Math.Abs(d[k]);
     309        }
     310        if (scale == 0.0) {
     311          e[i] = d[i - 1];
     312          for (int j = 0; j < i; j++) {
     313            d[j] = V[i - 1, j];
     314            V[i, j] = 0.0;
     315            V[j, i] = 0.0;
     316          }
     317        } else {
     318
     319          // Generate Householder vector.
     320
     321          for (int k = 0; k < i; k++) {
     322            d[k] /= scale;
     323            h += d[k] * d[k];
     324          }
     325          double f = d[i - 1];
     326          double g = Math.Sqrt(h);
     327          if (f > 0) {
     328            g = -g;
     329          }
     330          e[i] = scale * g;
     331          h = h - f * g;
     332          d[i - 1] = f - g;
     333          for (int j = 0; j < i; j++) {
     334            e[j] = 0.0;
     335          }
     336
     337          // Apply similarity transformation to remaining columns.
     338
     339          for (int j = 0; j < i; j++) {
     340            f = d[j];
     341            V[j, i] = f;
     342            g = e[j] + V[j, j] * f;
     343            for (int k = j + 1; k <= i - 1; k++) {
     344              g += V[k, j] * d[k];
     345              e[k] += V[k, j] * f;
     346            }
     347            e[j] = g;
     348          }
     349          f = 0.0;
     350          for (int j = 0; j < i; j++) {
     351            e[j] /= h;
     352            f += e[j] * d[j];
     353          }
     354          double hh = f / (h + h);
     355          for (int j = 0; j < i; j++) {
     356            e[j] -= hh * d[j];
     357          }
     358          for (int j = 0; j < i; j++) {
     359            f = d[j];
     360            g = e[j];
     361            for (int k = j; k <= i - 1; k++) {
     362              V[k, j] -= (f * e[k] + g * d[k]);
     363            }
     364            d[j] = V[i - 1, j];
     365            V[i, j] = 0.0;
     366          }
     367        }
     368        d[i] = h;
     369      }
     370
     371      // Accumulate transformations.
     372
     373      for (int i = 0; i < n - 1; i++) {
     374        V[n - 1, i] = V[i, i];
     375        V[i, i] = 1.0;
     376        double h = d[i + 1];
     377        if (h != 0.0) {
     378          for (int k = 0; k <= i; k++) {
     379            d[k] = V[k, i + 1] / h;
     380          }
     381          for (int j = 0; j <= i; j++) {
     382            double g = 0.0;
     383            for (int k = 0; k <= i; k++) {
     384              g += V[k, i + 1] * V[k, j];
     385            }
     386            for (int k = 0; k <= i; k++) {
     387              V[k, j] -= g * d[k];
     388            }
     389          }
     390        }
     391        for (int k = 0; k <= i; k++) {
     392          V[k, i + 1] = 0.0;
     393        }
     394      }
     395      for (int j = 0; j < n; j++) {
     396        d[j] = V[n - 1, j];
     397        V[n - 1, j] = 0.0;
     398      }
     399      V[n - 1, n - 1] = 1.0;
     400      e[0] = 0.0;
     401    }
     402
     403    // Symmetric tridiagonal QL algorithm, taken from JAMA package.
     404    private void tql2(int n, DoubleArray d, double[] e, DoubleMatrix V) {
     405
     406      //  This is derived from the Algol procedures tql2, by
     407      //  Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
     408      //  Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
     409      //  Fortran subroutine in EISPACK.
     410
     411      for (int i = 1; i < n; i++) {
     412        e[i - 1] = e[i];
     413      }
     414      e[n - 1] = 0.0;
     415
     416      double f = 0.0;
     417      double tst1 = 0.0;
     418      double eps = Math.Pow(2.0, -52.0);
     419      for (int l = 0; l < n; l++) {
     420
     421        // Find small subdiagonal element
     422
     423        tst1 = Math.Max(tst1, Math.Abs(d[l]) + Math.Abs(e[l]));
     424        int m = l;
     425        while (m < n) {
     426          if (Math.Abs(e[m]) <= eps * tst1) {
     427            break;
     428          }
     429          m++;
     430        }
     431
     432        // If m == l, d[l] is an eigenvalue,
     433        // otherwise, iterate.
     434
     435        if (m > l) {
     436          int iter = 0;
     437          do {
     438            iter = iter + 1;  // (Could check iteration count here.)
     439
     440            // Compute implicit shift
     441
     442            double g = d[l];
     443            double p = (d[l + 1] - g) / (2.0 * e[l]);
     444            double r = hypot(p, 1.0);
     445            if (p < 0) {
     446              r = -r;
     447            }
     448            d[l] = e[l] / (p + r);
     449            d[l + 1] = e[l] * (p + r);
     450            double dl1 = d[l + 1];
     451            double h = g - d[l];
     452            for (int i = l + 2; i < n; i++) {
     453              d[i] -= h;
     454            }
     455            f = f + h;
     456
     457            // Implicit QL transformation.
     458
     459            p = d[m];
     460            double c = 1.0;
     461            double c2 = c;
     462            double c3 = c;
     463            double el1 = e[l + 1];
     464            double s = 0.0;
     465            double s2 = 0.0;
     466            for (int i = m - 1; i >= l; i--) {
     467              c3 = c2;
     468              c2 = c;
     469              s2 = s;
     470              g = c * e[i];
     471              h = c * p;
     472              r = hypot(p, e[i]);
     473              e[i + 1] = s * r;
     474              s = e[i] / r;
     475              c = p / r;
     476              p = c * d[i] - s * g;
     477              d[i + 1] = h + s * (c * g + s * d[i]);
     478
     479              // Accumulate transformation.
     480
     481              for (int k = 0; k < n; k++) {
     482                h = V[k, i + 1];
     483                V[k, i + 1] = s * V[k, i] + c * h;
     484                V[k, i] = c * V[k, i] - s * h;
     485              }
     486            }
     487            p = -s * s2 * c3 * el1 * e[l] / dl1;
     488            e[l] = s * p;
     489            d[l] = c * p;
     490
     491            // Check for convergence.
     492
     493          } while (Math.Abs(e[l]) > eps * tst1);
     494        }
     495        d[l] = d[l] + f;
     496        e[l] = 0.0;
     497      }
     498
     499      // Sort eigenvalues and corresponding vectors.
     500
     501      for (int i = 0; i < n - 1; i++) {
     502        int k = i;
     503        double p = d[i];
     504        for (int j = i + 1; j < n; j++) {
     505          if (d[j] < p) { // NH find smallest k>i
     506            k = j;
     507            p = d[j];
     508          }
     509        }
     510        if (k != i) {
     511          d[k] = d[i]; // swap k and i
     512          d[i] = p;
     513          for (int j = 0; j < n; j++) {
     514            p = V[j, i];
     515            V[j, i] = V[j, k];
     516            V[j, k] = p;
     517          }
     518        }
     519      }
     520    }
     521
     522    /** sqrt(a^2 + b^2) without under/overflow. **/
     523    private double hypot(double a, double b) {
     524      double r = 0;
     525      if (Math.Abs(a) > Math.Abs(b)) {
     526        r = b / a;
     527        r = Math.Abs(a) * Math.Sqrt(1 + r * r);
     528      } else if (b != 0) {
     529        r = a / b;
     530        r = Math.Abs(b) * Math.Sqrt(1 + r * r);
     531      }
     532      return r;
     533    }
    262534  }
    263535}
  • branches/CMAES/HeuristicLab.Algorithms.CMAEvolutionStrategy/3.3/HeuristicLab.Algorithms.CMAEvolutionStrategy-3.3.csproj

    r9260 r9291  
    9393  </PropertyGroup>
    9494  <ItemGroup>
    95     <Reference Include="ALGLIB-3.6.0, Version=3.6.0.0, Culture=neutral, PublicKeyToken=ba48961d6f65dcec, processorArchitecture=MSIL">
    96       <SpecificVersion>False</SpecificVersion>
    97       <HintPath>..\..\..\..\trunk\sources\bin\ALGLIB-3.6.0.dll</HintPath>
    98       <Private>False</Private>
    99     </Reference>
    10095    <Reference Include="HeuristicLab.Analysis-3.3, Version=3.3.0.0, Culture=neutral, PublicKeyToken=ba48961d6f65dcec, processorArchitecture=MSIL">
    10196      <SpecificVersion>False</SpecificVersion>
  • branches/CMAES/HeuristicLab.Algorithms.CMAEvolutionStrategy/3.3/Plugin.cs.frame

    r9115 r9291  
    3333  [PluginDependency("HeuristicLab.Core", "3.3")]
    3434  [PluginDependency("HeuristicLab.Data", "3.3")]
     35  [PluginDependency("HeuristicLab.Encodings.RealVectorEncoding", "3.3")]
    3536  [PluginDependency("HeuristicLab.Operators", "3.3")]
    3637  [PluginDependency("HeuristicLab.Optimization", "3.3")]
  • branches/CMAES/HeuristicLab.Algorithms.CMAEvolutionStrategy/3.3/Properties/AssemblyInfo.cs.frame

    r9129 r9291  
    11#region License Information
    22/* HeuristicLab
    3  * Copyright (C) 2002-2012 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
     3 * Copyright (C) 2002-2013 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
    44 *
    55 * This file is part of HeuristicLab.
     
    3131[assembly: AssemblyCompany("")]
    3232[assembly: AssemblyProduct("HeuristicLab")]
    33 [assembly: AssemblyCopyright("(c) 2002-2012 HEAL")]
     33[assembly: AssemblyCopyright("(c) 2002-2013 HEAL")]
    3434[assembly: AssemblyTrademark("")]
    3535[assembly: AssemblyCulture("")]
Note: See TracChangeset for help on using the changeset viewer.