Changeset 17213 for branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Regression.Symbolic.Extensions/ConstrainedNLSInternal.cs

Timestamp:

08/13/19 19:01:54 (5 years ago)

Author:

gkronber

Message:

#2994: fixed a bug caused by cloning of trees, support other NLOpt solvers, implement idisposable, experiment with preconditioning (still not working)

File:

• branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Regression.Symbolic.Extensions/ConstrainedNLSInternal.cs (modified) (17 diffs)

branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Regression.Symbolic.Extensions/ConstrainedNLSInternal.cs

@@ -7 +7 @@
 using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
 
-namespace HeuristicLab.Problems.DataAnalysis.Symbolic.Regression.Extensions {
-  internal class ConstrainedNLSInternal {
+namespace HeuristicLab.Problems.DataAnalysis.Symbolic.Regression {
+  internal class ConstrainedNLSInternal : IDisposable {
     private readonly int maxIterations;
     public int MaxIterations => maxIterations;
@@ -58 +58 @@
     private readonly List<ConstantTreeNode>[] allThetaNodes;
     public List<ISymbolicExpressionTree> constraintTrees;    // TODO make local in ctor (public for debugging)
-    
+
     private readonly double[] fi_eval;
     private readonly double[,] jac_eval;
@@ -64 +64 @@
     private readonly VectorAutoDiffEvaluator autoDiffEval;
     private readonly VectorEvaluator eval;
+    private readonly bool invalidProblem = false;
 
     // end internal state
@@ -97 +98 @@
       var pred = interpreter.GetSymbolicExpressionTreeValues(expr, problemData.Dataset, trainingRows).ToArray();
 
-      if (pred.Any(pi => double.IsInfinity(pi) || double.IsNaN(pi))) throw new ArgumentException("The expression produces NaN or infinite values.");
+      if (pred.Any(pi => double.IsInfinity(pi) || double.IsNaN(pi))) {
+        bestError = targetVariance;
+        invalidProblem = true;
+      }
 
       #region linear scaling
       var predStDev = pred.StandardDeviationPop();
-      if (predStDev == 0) throw new ArgumentException("The expression is constant.");
+      if (predStDev == 0) {
+        bestError = targetVariance;
+        invalidProblem = true;
+      }
       var predMean = pred.Average();
 
@@ -118 +125 @@
       constraintTrees = new List<ISymbolicExpressionTree>();
       foreach (var constraint in intervalConstraints.Constraints) {
+        if (!constraint.Enabled) continue;
         if (constraint.IsDerivation) {
           if (!problemData.AllowedInputVariables.Contains(constraint.Variable))
@@ -130 +138 @@
             // convert variables named theta back to constants
             var df_prepared = ReplaceVarWithConst(df_smaller_upper, thetaNames, thetaValues, allThetaNodes);
-            constraintTrees.Add((ISymbolicExpressionTree)df_prepared.Clone());
+            constraintTrees.Add(df_prepared);
           }
           if (constraint.Interval.LowerBound > double.NegativeInfinity) {
             var df_larger_lower = Subtract(CreateConstant(constraint.Interval.LowerBound), (ISymbolicExpressionTree)df.Clone());
             // convert variables named theta back to constants
-            var df_prepared = ReplaceVarWithConst(df_larger_lower, thetaNames, thetaValues, allThetaNodes);            
-            constraintTrees.Add((ISymbolicExpressionTree)df_prepared.Clone());
+            var df_prepared = ReplaceVarWithConst(df_larger_lower, thetaNames, thetaValues, allThetaNodes);
+            constraintTrees.Add(df_prepared);
           }
         } else {
@@ -143 +151 @@
             // convert variables named theta back to constants
             var df_prepared = ReplaceVarWithConst(f_smaller_upper, thetaNames, thetaValues, allThetaNodes);
-            constraintTrees.Add((ISymbolicExpressionTree)df_prepared.Clone());
+            constraintTrees.Add(df_prepared);
           }
           if (constraint.Interval.LowerBound > double.NegativeInfinity) {
@@ -149 +157 @@
             // convert variables named theta back to constants
             var df_prepared = ReplaceVarWithConst(f_larger_lower, thetaNames, thetaValues, allThetaNodes);
-            constraintTrees.Add((ISymbolicExpressionTree)df_prepared.Clone());
+            constraintTrees.Add(df_prepared);
           }
         }
@@ -173 +181 @@
       NLOpt.nlopt_set_min_objective(nlopt, calculateObjectiveDelegate, IntPtr.Zero); // --> without preconditioning
 
-      // preconditionDelegate = new NLOpt.nlopt_precond(PreconditionObjective);
-      // NLOpt.nlopt_set_precond_min_objective(nlopt, calculateObjectiveDelegate, preconditionDelegate, IntPtr.Zero);
+      //preconditionDelegate = new NLOpt.nlopt_precond(PreconditionObjective);
+      //NLOpt.nlopt_set_precond_min_objective(nlopt, calculateObjectiveDelegate, preconditionDelegate, IntPtr.Zero);
 
 
@@ -185 +193 @@
         calculateConstraintDelegates[i] = new NLOpt.nlopt_func(CalculateConstraint);
         NLOpt.nlopt_add_inequality_constraint(nlopt, calculateConstraintDelegates[i], constraintDataPtr[i], 1e-8);
-        //NLOpt.nlopt_add_precond_inequality_constraint(nlopt, calculateConstraintDelegates[i], preconditionDelegate, constraintDataPtr[i], 1e-8);
+        // NLOpt.nlopt_add_precond_inequality_constraint(nlopt, calculateConstraintDelegates[i], preconditionDelegate, constraintDataPtr[i], 1e-8);
       }
 
@@ -196 +204 @@
 
     ~ConstrainedNLSInternal() {
-      if (nlopt != IntPtr.Zero)
-        NLOpt.nlopt_destroy(nlopt);
-      if (constraintDataPtr != null) {
-        for (int i = 0; i < constraintDataPtr.Length; i++)
-          if (constraintDataPtr[i] != IntPtr.Zero)
-            Marshal.FreeHGlobal(constraintDataPtr[i]);
-      }
+      Dispose();
     }
 
 
     internal void Optimize() {
+      if (invalidProblem) return;
       var x = thetaValues.ToArray();  /* initial guess */
       double minf = double.MaxValue; /* minimum objective value upon return */
@@ -219 +222 @@
         double[] _ = new double[x.Length];
         bestConstraintValues = new double[calculateConstraintDelegates.Length];
-        for(int i=0;i<calculateConstraintDelegates.Length;i++) {
+        for (int i = 0; i < calculateConstraintDelegates.Length; i++) {
           bestConstraintValues[i] = calculateConstraintDelegates[i].Invoke((uint)x.Length, x, _, constraintDataPtr[i]);
         }
@@ -300 +303 @@
     }
 
-    // // NOT WORKING YET? WHY is det(H) always zero?!
-    // // TODO
+    // TODO
     // private void PreconditionObjective(uint n, double[] x, double[] v, double[] vpre, IntPtr data) {
     //   UpdateThetaValues(x); // calc H(x)
@@ -311 +313 @@
     //   
     //   // calc residuals and scale jac_eval
-    //   var f = -2.0 / k;
-    //   for (int i = 0;i<k;i++) {
-    //     var r = target[i] - fi_eval[i];
-    //     for(int j = 0;j<n;j++) {
-    //       jac_eval[i, j] *= f * r;
-    //     }
-    //   }
-    //   
+    //   var f = 2.0 / (k*k);
+    // 
     //   // approximate hessian H(x) = J(x)^T * J(x)
     //   alglib.rmatrixgemm((int)n, (int)n, k,
-    //     1.0, jac_eval, 0, 0, 1,  // transposed
+    //     f, jac_eval, 0, 0, 1,  // transposed
     //     jac_eval, 0, 0, 0, 
     //     0.0, ref h, 0, 0,
@@ -327 +323 @@
     //     );
     //   
+    // 
     //   // scale v
     //   alglib.rmatrixmv((int)n, (int)n, h, 0, 0, 0, v, 0, ref vpre, 0, alglib.serial);
-    //   var det = alglib.matdet.rmatrixdet(h, (int)n, alglib.serial);
+    // 
+    // 
+    //   alglib.spdmatrixcholesky(ref h, (int)n, true);
+    // 
+    //   var det = alglib.matdet.spdmatrixcholeskydet(h, (int)n, alglib.serial);
     // }
 
@@ -446 +447 @@
       if (solver.Contains("CCSAQ")) return NLOpt.nlopt_algorithm.NLOPT_LD_CCSAQ;
       if (solver.Contains("ISRES")) return NLOpt.nlopt_algorithm.NLOPT_GN_ISRES;
+
+      if (solver.Contains("DIRECT_G")) return NLOpt.nlopt_algorithm.NLOPT_GN_DIRECT;
+      if (solver.Contains("NLOPT_GN_DIRECT_L")) return NLOpt.nlopt_algorithm.NLOPT_GN_DIRECT_L;
+      if (solver.Contains("NLOPT_GN_DIRECT_L_RAND")) return NLOpt.nlopt_algorithm.NLOPT_GN_DIRECT_L_RAND;
+      if (solver.Contains("NLOPT_GN_ORIG_DIRECT")) return NLOpt.nlopt_algorithm.NLOPT_GN_DIRECT;
+      if (solver.Contains("NLOPT_GN_ORIG_DIRECT_L")) return NLOpt.nlopt_algorithm.NLOPT_GN_ORIG_DIRECT_L;
+      if (solver.Contains("NLOPT_GD_STOGO")) return NLOpt.nlopt_algorithm.NLOPT_GD_STOGO;
+      if (solver.Contains("NLOPT_GD_STOGO_RAND")) return NLOpt.nlopt_algorithm.NLOPT_GD_STOGO_RAND;
+      if (solver.Contains("NLOPT_LD_LBFGS_NOCEDAL")) return NLOpt.nlopt_algorithm.NLOPT_LD_LBFGS_NOCEDAL;
+      if (solver.Contains("NLOPT_LD_LBFGS")) return NLOpt.nlopt_algorithm.NLOPT_LD_LBFGS;
+      if (solver.Contains("NLOPT_LN_PRAXIS")) return NLOpt.nlopt_algorithm.NLOPT_LN_PRAXIS;
+      if (solver.Contains("NLOPT_LD_VAR1")) return NLOpt.nlopt_algorithm.NLOPT_LD_VAR1;
+      if (solver.Contains("NLOPT_LD_VAR2")) return NLOpt.nlopt_algorithm.NLOPT_LD_VAR2;
+      if (solver.Contains("NLOPT_LD_TNEWTON")) return NLOpt.nlopt_algorithm.NLOPT_LD_TNEWTON;
+      if (solver.Contains("NLOPT_LD_TNEWTON_RESTART")) return NLOpt.nlopt_algorithm.NLOPT_LD_TNEWTON_RESTART;
+      if (solver.Contains("NLOPT_LD_TNEWTON_PRECOND")) return NLOpt.nlopt_algorithm.NLOPT_LD_TNEWTON_PRECOND;
+      if (solver.Contains("NLOPT_LD_TNEWTON_PRECOND_RESTART")) return NLOpt.nlopt_algorithm.NLOPT_LD_TNEWTON_PRECOND_RESTART;
+      if (solver.Contains("NLOPT_GN_CRS2_LM")) return NLOpt.nlopt_algorithm.NLOPT_GN_CRS2_LM;
+      if (solver.Contains("NLOPT_GN_MLSL")) return NLOpt.nlopt_algorithm.NLOPT_GN_MLSL;
+      if (solver.Contains("NLOPT_GD_MLSL")) return NLOpt.nlopt_algorithm.NLOPT_GD_MLSL;
+      if (solver.Contains("NLOPT_GN_MLSL_LDS")) return NLOpt.nlopt_algorithm.NLOPT_GN_MLSL_LDS;
+      if (solver.Contains("NLOPT_GD_MLSL_LDS")) return NLOpt.nlopt_algorithm.NLOPT_GD_MLSL_LDS;
+      if (solver.Contains("NLOPT_LN_NEWUOA")) return NLOpt.nlopt_algorithm.NLOPT_LN_NEWUOA;
+      if (solver.Contains("NLOPT_LN_NEWUOA_BOUND")) return NLOpt.nlopt_algorithm.NLOPT_LN_NEWUOA_BOUND;
+      if (solver.Contains("NLOPT_LN_NELDERMEAD")) return NLOpt.nlopt_algorithm.NLOPT_LN_NELDERMEAD;
+      if (solver.Contains("NLOPT_LN_SBPLX")) return NLOpt.nlopt_algorithm.NLOPT_LN_SBPLX;
+      if (solver.Contains("NLOPT_LN_AUGLAG")) return NLOpt.nlopt_algorithm.NLOPT_LN_AUGLAG;
+      if (solver.Contains("NLOPT_LD_AUGLAG")) return NLOpt.nlopt_algorithm.NLOPT_LD_AUGLAG;
+      if (solver.Contains("NLOPT_LN_BOBYQA")) return NLOpt.nlopt_algorithm.NLOPT_LN_BOBYQA;
+      if (solver.Contains("NLOPT_AUGLAG")) return NLOpt.nlopt_algorithm.NLOPT_AUGLAG;
+      if (solver.Contains("NLOPT_LD_SLSQP")) return NLOpt.nlopt_algorithm.NLOPT_LD_SLSQP;
+      if (solver.Contains("NLOPT_LD_CCSAQ))")) return NLOpt.nlopt_algorithm.NLOPT_LD_CCSAQ;
+      if (solver.Contains("NLOPT_GN_ESCH")) return NLOpt.nlopt_algorithm.NLOPT_GN_ESCH;
+      if (solver.Contains("NLOPT_GN_AGS")) return NLOpt.nlopt_algorithm.NLOPT_GN_AGS;
+
       throw new ArgumentException($"Unknown solver {solver}");
     }
@@ -559 +595 @@
       return node;
     }
+
+    public void Dispose() {
+      if (nlopt != IntPtr.Zero) {
+        NLOpt.nlopt_destroy(nlopt);
+        nlopt = IntPtr.Zero;
+      }
+      if (constraintDataPtr != null) {
+        for (int i = 0; i < constraintDataPtr.Length; i++)
+          if (constraintDataPtr[i] != IntPtr.Zero) {
+            Marshal.FreeHGlobal(constraintDataPtr[i]);
+            constraintDataPtr[i] = IntPtr.Zero;
+          }
+      }
+    }
     #endregion
   }