Changeset 9291

Timestamp:

03/06/13 17:33:38 (12 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:

• CMAOperators/CMAMutator.cs (modified) (5 diffs)
• CMAOperators/CMAUpdater.cs (modified) (3 diffs)
• HeuristicLab.Algorithms.CMAEvolutionStrategy-3.3.csproj (modified) (1 diff)
• Plugin.cs.frame (modified) (1 diff)
• Properties/AssemblyInfo.cs.frame (modified) (2 diffs)

branches/CMAES/HeuristicLab.Algorithms.CMAEvolutionStrategy/3.3/CMAOperators/CMAMutator.cs

@@ -20 +20 @@
 #endregion
 
-using System;
 using HeuristicLab.Common;
 using HeuristicLab.Core;
@@ -29 +28 @@
 using HeuristicLab.Parameters;
 using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
+using HeuristicLab.Random;
+using System;
 
 namespace HeuristicLab.Algorithms.CMAEvolutionStrategy {
@@ -113 +114 @@
         RealVectorParameter.ActualValue = arx;
       }
-      alglib.hqrndstate state;
-      alglib.hqrndseed(random.Next(), random.Next(), out state);
+      var nd = new NormalDistributedRandom(random, 0, 1);
+      //alglib.hqrndstate state;
+      //alglib.hqrndseed(random.Next(), random.Next(), out state);
 
       for (int i = 0; i < lambda; i++) {
@@ -122 +124 @@
           for (int k = 0; k < arx[i].Length; k++) {
             do {
-              arx[i][k] = xmean[k] + sp.Sigma.Value * sp.D[k] * alglib.hqrndnormal(state);
+              arx[i][k] = xmean[k] + sp.Sigma.Value * sp.D[k] * nd.NextDouble(); //alglib.hqrndnormal(state);
               inRange = bounds[k % bounds.Rows, 0] <= arx[i][k] && arx[i][k] <= bounds[k % bounds.Rows, 1];
               if (!inRange) tries++;
@@ -138 +140 @@
             var artmp = new double[arx[0].Length];
             for (int k = 0; k < arx[0].Length; ++k) {
-              artmp[k] = sp.D[k] * alglib.hqrndnormal(state);
+              artmp[k] = sp.D[k] * nd.NextDouble(); //alglib.hqrndnormal(state);
             }
 

branches/CMAES/HeuristicLab.Algorithms.CMAEvolutionStrategy/3.3/CMAOperators/CMAUpdater.cs

@@ -237 +237 @@
         DegenerateStateParameter.ActualValue = new BoolValue(false);
       } else {
-        var c = new double[sp.C.Rows, sp.C.Columns];
+        /*var c = new double[sp.C.Rows, sp.C.Columns];
         for (int i = 0; i < sp.C.Rows; i++)
           for (int j = 0; j < sp.C.Columns; j++)
@@ -243 +243 @@
 
         double[] eVal;
-        double[,] eVec;
-        var success = alglib.smatrixevd(c, N, 1, false, out eVal, out eVec);
+        double[,] eVec;*/
+
+        // set B <- C
+        for (int i = 0; i < N; ++i)
+          for (int j = 0; j <= i; ++j)
+            sp.B[i, j] = sp.B[j, i] = sp.C[i, j];
+
+        var success = Eigendecomposition(N, sp.B, sp.D); // alglib.smatrixevd(c, N, 1, false, out eVal, out eVec);
         DegenerateStateParameter.ActualValue = new BoolValue(!success);
-        sp.B = new DoubleMatrix(eVec);
-
-        for (int i = 0; i < sp.D.Length; i++) {
+        //sp.B = new DoubleMatrix(eVec);
+
+        /*for (int i = 0; i < sp.D.Length; i++) {
           if (eVal[i] < 0) {
             DegenerateStateParameter.ActualValue.Value = true;
             sp.D[i] = 0;
           } else sp.D[i] = Math.Sqrt(eVal[i]);
-        }
+        }*/
 
         if (sp.D.Min() == 0.0) sp.AxisRatio.Value = double.PositiveInfinity;
@@ -260 +266 @@
       return base.Apply();
     }
+
+    private bool Eigendecomposition(int N, DoubleMatrix B, DoubleArray diagD) {
+      bool result = true;
+      // eigendecomposition
+      var offdiag = new double[N];
+      tred2(N, B, diagD, offdiag);
+      tql2(N, diagD, offdiag, B);
+
+      // assign diagD to eigenvalue square roots
+      for (int i = 0; i < N; ++i) {
+        if (diagD[i] < 0) // numerical problem?
+          result = false;
+        else diagD[i] = Math.Sqrt(diagD[i]);
+      }
+
+      return result;
+    } // eigendecomposition
+
+
+    // Symmetric Householder reduction to tridiagonal form, taken from JAMA package.
+    private void tred2(int n, DoubleMatrix V, DoubleArray d, double[] e) {
+
+      //  This is derived from the Algol procedures tred2 by
+      //  Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
+      //  Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
+      //  Fortran subroutine in EISPACK.
+
+      for (int j = 0; j < n; j++) {
+        d[j] = V[n - 1, j];
+      }
+
+      // Householder reduction to tridiagonal form.
+
+      for (int i = n - 1; i > 0; i--) {
+
+        // Scale to avoid under/overflow.
+
+        double scale = 0.0;
+        double h = 0.0;
+        for (int k = 0; k < i; k++) {
+          scale = scale + Math.Abs(d[k]);
+        }
+        if (scale == 0.0) {
+          e[i] = d[i - 1];
+          for (int j = 0; j < i; j++) {
+            d[j] = V[i - 1, j];
+            V[i, j] = 0.0;
+            V[j, i] = 0.0;
+          }
+        } else {
+
+          // Generate Householder vector.
+
+          for (int k = 0; k < i; k++) {
+            d[k] /= scale;
+            h += d[k] * d[k];
+          }
+          double f = d[i - 1];
+          double g = Math.Sqrt(h);
+          if (f > 0) {
+            g = -g;
+          }
+          e[i] = scale * g;
+          h = h - f * g;
+          d[i - 1] = f - g;
+          for (int j = 0; j < i; j++) {
+            e[j] = 0.0;
+          }
+
+          // Apply similarity transformation to remaining columns.
+
+          for (int j = 0; j < i; j++) {
+            f = d[j];
+            V[j, i] = f;
+            g = e[j] + V[j, j] * f;
+            for (int k = j + 1; k <= i - 1; k++) {
+              g += V[k, j] * d[k];
+              e[k] += V[k, j] * f;
+            }
+            e[j] = g;
+          }
+          f = 0.0;
+          for (int j = 0; j < i; j++) {
+            e[j] /= h;
+            f += e[j] * d[j];
+          }
+          double hh = f / (h + h);
+          for (int j = 0; j < i; j++) {
+            e[j] -= hh * d[j];
+          }
+          for (int j = 0; j < i; j++) {
+            f = d[j];
+            g = e[j];
+            for (int k = j; k <= i - 1; k++) {
+              V[k, j] -= (f * e[k] + g * d[k]);
+            }
+            d[j] = V[i - 1, j];
+            V[i, j] = 0.0;
+          }
+        }
+        d[i] = h;
+      }
+
+      // Accumulate transformations.
+
+      for (int i = 0; i < n - 1; i++) {
+        V[n - 1, i] = V[i, i];
+        V[i, i] = 1.0;
+        double h = d[i + 1];
+        if (h != 0.0) {
+          for (int k = 0; k <= i; k++) {
+            d[k] = V[k, i + 1] / h;
+          }
+          for (int j = 0; j <= i; j++) {
+            double g = 0.0;
+            for (int k = 0; k <= i; k++) {
+              g += V[k, i + 1] * V[k, j];
+            }
+            for (int k = 0; k <= i; k++) {
+              V[k, j] -= g * d[k];
+            }
+          }
+        }
+        for (int k = 0; k <= i; k++) {
+          V[k, i + 1] = 0.0;
+        }
+      }
+      for (int j = 0; j < n; j++) {
+        d[j] = V[n - 1, j];
+        V[n - 1, j] = 0.0;
+      }
+      V[n - 1, n - 1] = 1.0;
+      e[0] = 0.0;
+    }
+
+    // Symmetric tridiagonal QL algorithm, taken from JAMA package.
+    private void tql2(int n, DoubleArray d, double[] e, DoubleMatrix V) {
+
+      //  This is derived from the Algol procedures tql2, by
+      //  Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
+      //  Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
+      //  Fortran subroutine in EISPACK.
+
+      for (int i = 1; i < n; i++) {
+        e[i - 1] = e[i];
+      }
+      e[n - 1] = 0.0;
+
+      double f = 0.0;
+      double tst1 = 0.0;
+      double eps = Math.Pow(2.0, -52.0);
+      for (int l = 0; l < n; l++) {
+
+        // Find small subdiagonal element
+
+        tst1 = Math.Max(tst1, Math.Abs(d[l]) + Math.Abs(e[l]));
+        int m = l;
+        while (m < n) {
+          if (Math.Abs(e[m]) <= eps * tst1) {
+            break;
+          }
+          m++;
+        }
+
+        // If m == l, d[l] is an eigenvalue,
+        // otherwise, iterate.
+
+        if (m > l) {
+          int iter = 0;
+          do {
+            iter = iter + 1;  // (Could check iteration count here.)
+
+            // Compute implicit shift
+
+            double g = d[l];
+            double p = (d[l + 1] - g) / (2.0 * e[l]);
+            double r = hypot(p, 1.0);
+            if (p < 0) {
+              r = -r;
+            }
+            d[l] = e[l] / (p + r);
+            d[l + 1] = e[l] * (p + r);
+            double dl1 = d[l + 1];
+            double h = g - d[l];
+            for (int i = l + 2; i < n; i++) {
+              d[i] -= h;
+            }
+            f = f + h;
+
+            // Implicit QL transformation.
+
+            p = d[m];
+            double c = 1.0;
+            double c2 = c;
+            double c3 = c;
+            double el1 = e[l + 1];
+            double s = 0.0;
+            double s2 = 0.0;
+            for (int i = m - 1; i >= l; i--) {
+              c3 = c2;
+              c2 = c;
+              s2 = s;
+              g = c * e[i];
+              h = c * p;
+              r = hypot(p, e[i]);
+              e[i + 1] = s * r;
+              s = e[i] / r;
+              c = p / r;
+              p = c * d[i] - s * g;
+              d[i + 1] = h + s * (c * g + s * d[i]);
+
+              // Accumulate transformation.
+
+              for (int k = 0; k < n; k++) {
+                h = V[k, i + 1];
+                V[k, i + 1] = s * V[k, i] + c * h;
+                V[k, i] = c * V[k, i] - s * h;
+              }
+            }
+            p = -s * s2 * c3 * el1 * e[l] / dl1;
+            e[l] = s * p;
+            d[l] = c * p;
+
+            // Check for convergence.
+
+          } while (Math.Abs(e[l]) > eps * tst1);
+        }
+        d[l] = d[l] + f;
+        e[l] = 0.0;
+      }
+
+      // Sort eigenvalues and corresponding vectors.
+
+      for (int i = 0; i < n - 1; i++) {
+        int k = i;
+        double p = d[i];
+        for (int j = i + 1; j < n; j++) {
+          if (d[j] < p) { // NH find smallest k>i
+            k = j;
+            p = d[j];
+          }
+        }
+        if (k != i) {
+          d[k] = d[i]; // swap k and i 
+          d[i] = p;
+          for (int j = 0; j < n; j++) {
+            p = V[j, i];
+            V[j, i] = V[j, k];
+            V[j, k] = p;
+          }
+        }
+      }
+    }
+
+    /** sqrt(a^2 + b^2) without under/overflow. **/
+    private double hypot(double a, double b) {
+      double r = 0;
+      if (Math.Abs(a) > Math.Abs(b)) {
+        r = b / a;
+        r = Math.Abs(a) * Math.Sqrt(1 + r * r);
+      } else if (b != 0) {
+        r = a / b;
+        r = Math.Abs(b) * Math.Sqrt(1 + r * r);
+      }
+      return r;
+    }
   }
 }

branches/CMAES/HeuristicLab.Algorithms.CMAEvolutionStrategy/3.3/HeuristicLab.Algorithms.CMAEvolutionStrategy-3.3.csproj

@@ -93 +93 @@
   </PropertyGroup>
   <ItemGroup>
-    <Reference Include="ALGLIB-3.6.0, Version=3.6.0.0, Culture=neutral, PublicKeyToken=ba48961d6f65dcec, processorArchitecture=MSIL">
-      <SpecificVersion>False</SpecificVersion>
-      <HintPath>..\..\..\..\trunk\sources\bin\ALGLIB-3.6.0.dll</HintPath>
-      <Private>False</Private>
-    </Reference>
     <Reference Include="HeuristicLab.Analysis-3.3, Version=3.3.0.0, Culture=neutral, PublicKeyToken=ba48961d6f65dcec, processorArchitecture=MSIL">
       <SpecificVersion>False</SpecificVersion>

branches/CMAES/HeuristicLab.Algorithms.CMAEvolutionStrategy/3.3/Plugin.cs.frame

@@ -33 +33 @@
   [PluginDependency("HeuristicLab.Core", "3.3")]
   [PluginDependency("HeuristicLab.Data", "3.3")]
+  [PluginDependency("HeuristicLab.Encodings.RealVectorEncoding", "3.3")]
   [PluginDependency("HeuristicLab.Operators", "3.3")]
   [PluginDependency("HeuristicLab.Optimization", "3.3")]

branches/CMAES/HeuristicLab.Algorithms.CMAEvolutionStrategy/3.3/Properties/AssemblyInfo.cs.frame

@@ -1 +1 @@
 #region License Information
 /* HeuristicLab
- * Copyright (C) 2002-2012 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
+ * Copyright (C) 2002-2013 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
  *
  * This file is part of HeuristicLab.
@@ -31 +31 @@
 [assembly: AssemblyCompany("")]
 [assembly: AssemblyProduct("HeuristicLab")]
-[assembly: AssemblyCopyright("(c) 2002-2012 HEAL")]
+[assembly: AssemblyCopyright("(c) 2002-2013 HEAL")]
 [assembly: AssemblyTrademark("")]
 [assembly: AssemblyCulture("")]