From a09e8d5fc79fa73883a3eaddf867ae6ae02d7a10 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Fran=C3=A7ois=20Henrotte?= <francois.henrotte@uclouvain.be>
Date: Wed, 17 Oct 2018 09:19:23 +0200
Subject: [PATCH] mma and scipy as options on the same footing

---
 Ccore/ccore.pro        |  10 +--
 Ccore/ccore_common.pro |   4 ++
 Ccore/shp.py           | 151 +++++++++++++++++++----------------------
 3 files changed, 78 insertions(+), 87 deletions(-)

diff --git a/Ccore/ccore.pro b/Ccore/ccore.pro
index 620a355..020d7af 100644
--- a/Ccore/ccore.pro
+++ b/Ccore/ccore.pro
@@ -387,7 +387,7 @@ PostOperation {
 
       CreateDir["res"];
       //Print[js, OnElementsOf Vol_S_Mag, File "res/js.pos"];
-      Print[b, OnElementsOf Vol_Mag, File "res/b.pos"];
+      //Print[b, OnElementsOf Vol_Mag, File "res/b.pos"];
 
       Print[az, OnElementsOf Vol_Mag, File "res/az.pos"];
       Echo[ Str["k= PostProcessing.NbViews-1;",
@@ -501,10 +501,10 @@ PostOperation{
 	    File "tmp.geo", LastTimeStepOnly] ;
 
       //Print[az, OnElementsOf Vol_Mag, File "res/azstar.pos"];
-      Echo[ Str["k= PostProcessing.NbViews-1;",
-		"View[k].IntervalsType = 3;",
-		"View[0].NbIso = 30;"], 
-	    File "tmp.geo", LastTimeStepOnly] ;
+      // Echo[ Str["k= PostProcessing.NbViews-1;",
+      // 		"View[k].IntervalsType = 3;",
+      // 		"View[0].NbIso = 30;"], 
+      // 	    File "tmp.geo", LastTimeStepOnly] ;
 
       Print[ Lie_az, OnElementsOf Vol_Mag, File "res/Lie_az.pos" ];
       Echo[ Str["k= PostProcessing.NbViews-1;",
diff --git a/Ccore/ccore_common.pro b/Ccore/ccore_common.pro
index 6ef1cab..c780494 100644
--- a/Ccore/ccore_common.pro
+++ b/Ccore/ccore_common.pro
@@ -6,6 +6,10 @@ ParamIndex =
 	       Choices {0="Optimize on D and E", 
 			1="D: air gap length", 
 	       		2="E: C-core thickness"} ];
+Optimizer = 
+  DefineNumber[0, Name "Optimization/Optimizer",
+	       Choices {0="scipy", 
+			1="MMA"}];
 StepsInRange = 
   DefineNumber[20, Name "Optimization/Steps in range", Visible ParamIndex];
 
diff --git a/Ccore/shp.py b/Ccore/shp.py
index ed42af8..e9eaaca 100644
--- a/Ccore/shp.py
+++ b/Ccore/shp.py
@@ -8,6 +8,7 @@ import conveks
 import onelab
 from scipy.optimize import minimize
 
+modelName = 'ccore'
 
 c = onelab.client(__file__)
 
@@ -27,8 +28,6 @@ if(not len(mygetdp)):
 
 c.sendInfo('Will use gmsh={0} and getdp={1}'.format(mygmsh, mygetdp))
 
-modelName = 'ccore'
-
 # append correct path to model file names
 file_geo = c.getPath(modelName+'.geo')
 file_msh = c.getPath(modelName+'.msh')
@@ -46,10 +45,6 @@ c.sendCommand('Solver.AutoMesh = -1;')
 # dry getdp run (without -solve or -pos option) to get model parameters in the GUI
 c.runSubClient('myGetDP', mygetdp)
 
-# define optimization parameters as Onelab parameter (editable in the GUI)
-maxIter = c.defineNumber('Optimization/00Max iterations', value=100)
-maxChange = c.defineNumber('Optimization/01Max change', value=1e-0)
-
 # end of check phase. We're done if we do not start the optimization
 if c.action == 'check' :
    exit(0)
@@ -145,13 +140,22 @@ def test(n):
 # List of design parameters with format
 # [ "Onelab name", value, Min, Max, "variable name" ]
 designParameters={
-    0 : ['Optimization/Parameters/D', 0.01, 0.001, 0.050, 'D'],
+    0 : ['Optimization/Parameters/D', 0.02, 0.001, 0.050, 'D'],
     1 : ['Optimization/Parameters/E', 0.01 , 0.01, 0.15 , 'E']
 }
-
+Optimizers={
+    0 : 'scipy',
+    1 : 'mma'
+}
 
 x = {}
 index = int(c.getNumber('Optimization/Sweep on parameter'))-1
+Optimizer = Optimizers[int(c.getNumber('Optimization/Optimizer'))]
+
+# define optimization parameters as Onelab parameter (editable in the GUI)
+maxIter = c.defineNumber('Optimization/00Max iterations', value=100, visible=(Optimizer=='mma'))
+maxChange = c.defineNumber('Optimization/01Max change', value=1e-1,  visible=(Optimizer=='mma'))
+
 
 if index >= 0: # parameter sweep and optimisation over one selected design parameter
     x[0] = designParameters[index];
@@ -162,85 +166,68 @@ if index >= 0: # parameter sweep and optimisation over one selected design param
     print res
 
 else: # optimisation over full parameter space
-    for label, var in designParameters.iteritems():
-        x[label] = var
-    values = [designParameters[i][1] for i in range(len(designParameters))]
-    bounds = [(designParameters[i][2], designParameters[i][3]) for i in range(len(designParameters))]
-    ## print x
-    ## print values, bounds
-    ## print fobj(values)
-    ## print grad(values)
-    ## exit(1)
-
-    res = minimize(fobj, values, jac=grad, bounds=bounds, callback=callbackF,\
-                method = 'L-BFGS-B', tol=None, options={'ftol': 1e-5, 'gtol': 1e-3, 'disp':True} )
-    print res
-    
-    
-    
-exit(1)
-
-# number of design variables and number of constraints
-numVariables = len(x.keys())
-
-# initial value of the variables in the design space,
-# the lower bound for design variables,
-# as well as the upper bound for design variables.
-initialPoint = np.zeros(numVariables)
-lowerBound = np.zeros(numVariables)
-upperBound = np.zeros(numVariables)
-for label, var in x.iteritems():
-  initialPoint[label] = var[1]
-  lowerBound[label]   = var[2]
-  upperBound[label]   = var[3]
 
+    if Optimizer=='scipy':
+        for label, var in designParameters.iteritems():
+            x[label] = var
+        values = [designParameters[i][1] for i in range(len(designParameters))]
+        bounds = [(designParameters[i][2], designParameters[i][3]) for i in range(len(designParameters))]
 
+        res = minimize(fobj, values, jac=grad, bounds=bounds, callback=callbackF,\
+                method = 'L-BFGS-B', tol=None, options={'ftol': 1e-5, 'gtol': 1e-3, 'disp':True} )
+        print res
+
+    elif Optimizer=='mma':
+        numVariables = len(designParameters.keys())
+        initialPoint = np.zeros(numVariables)
+        lowerBound = np.zeros(numVariables)
+        upperBound = np.zeros(numVariables)
+        for label, var in designParameters.iteritems():
+            x[label] = var
+            initialPoint[label] = var[1]
+            lowerBound[label]   = var[2]
+            upperBound[label]   = var[3]
   
-# Initialize the MMA optimizer
-conveks.mma.initialize(initialPoint, lowerBound, upperBound)
-
-# Set some options for MMA
-conveks.mma.option.setNumber('General.Verbosity', 0) # def 4
-conveks.mma.option.setNumber('General.SaveHistory', 0)
-conveks.mma.option.setNumber('SubProblem.isRobustMoveLimit', 1) # def 1
-conveks.mma.option.setNumber('SubProblem.isRobustAsymptotes', 1) # def 1  
-conveks.mma.option.setNumber('SubProblem.type', 0)
-conveks.mma.option.setNumber('SubProblem.addConvexity', 0)
-conveks.mma.option.setNumber('SubProblem.asymptotesRmax', 100.0)
-conveks.mma.option.setNumber('SubProblem.asymptotesRmin', 0.001) # def 0.001
-
-# Get iteration count (here it will be 1 - could be different in case of restart)
-it = conveks.mma.getOuterIteration()
-change = 1.0
-
-while it <= maxIter and c.getString(c.name+'/Action') != 'stop': #and change > 1e-6 :
-    c.setNumber('Optimization/01Current iteration', value=it, readOnly=1,
-            attributes={'Highlight':'LightYellow'})
-
-    # get (copy of) current point
-    xFromMMA = np.array(conveks.mma.getCurrentPoint())
-    print 'xFromMMA = ', xFromMMA
+        # Initialize the MMA optimizer
+        conveks.mma.initialize(initialPoint, lowerBound, upperBound)
+
+        # Set some options for MMA
+        conveks.mma.option.setNumber('General.Verbosity', 0) # def 4
+        conveks.mma.option.setNumber('General.SaveHistory', 0)
+        conveks.mma.option.setNumber('SubProblem.isRobustMoveLimit', 1) # def 1
+        conveks.mma.option.setNumber('SubProblem.isRobustAsymptotes', 1) # def 1  
+        conveks.mma.option.setNumber('SubProblem.type', 0)
+        conveks.mma.option.setNumber('SubProblem.addConvexity', 0)
+        conveks.mma.option.setNumber('SubProblem.asymptotesRmax', 100.0)
+        conveks.mma.option.setNumber('SubProblem.asymptotesRmin', 0.001) # def 0.001
+
+        # Get iteration count (here it will be 1 - could be different in case of restart)
+        it = conveks.mma.getOuterIteration()
+        change = 1.0
+        while it <= maxIter and c.getString(c.name+'/Action') != 'stop' and change > 1e-4 :
+            c.setNumber('Optimization/01Current iteration', value=it, readOnly=1,
+                         attributes={'Highlight':'LightYellow'})
+
+            xFromMMA = np.array(conveks.mma.getCurrentPoint())
+            objective = fobj(xFromMMA)
+            grad_objective = grad(xFromMMA)
+            constraints = np.array([np.sum(xFromMMA)/100.0-1.0])
+            grad_constraints = np.ones(numVariables)/100.0
     
-    # get the value of the objective function and of the constraints
-    # as well as their sensitivity with respect to design variables at `x'
-    objective = fobj(xFromMMA)
-    grad_objective = grad(xFromMMA)
-    constraints = np.array([np.sum(xFromMMA)/100.0-1.0])
-    grad_constraints = np.ones(numVariables)/100.0
-    
-    if it == 1: fscale = 1.0 / np.abs(objective)
-    objective *= fscale
-    grad_objective *= fscale
+            if it == 1: fscale = 1.0 / np.abs(objective)
+            objective *= fscale
+            grad_objective *= fscale
+
+            print it, change, xFromMMA, objective, grad_objective
+            c.sendInfo('Optimization: it. {}, obj. {}, constr. {}'.format(it,objective,constraints[0]))
 
-    print it, xFromMMA, objective, grad_objective
-    c.sendInfo('Optimization: it. {}, obj. {}, constr. {}'.format(it,objective,constraints[0]))
+            # call MMA and update the current point
+            conveks.mma.setMoveLimits(lowerBound, upperBound, 1e-0) # parametre interessant
+            conveks.mma.updateCurrentPoint(constraints,grad_objective,grad_constraints)
+            change = conveks.mma.getDesignChange()
+            it = conveks.mma.getOuterIteration()
 
-    # call MMA and update the current point
-    conveks.mma.setMoveLimits(lowerBound, upperBound, 1e-0) # parametre interessant
-    conveks.mma.updateCurrentPoint(constraints,grad_objective,grad_constraints)
-    change = conveks.mma.getDesignChange()
-    it = conveks.mma.getOuterIteration()
+        # This should be called at the end
+        conveks.mma.finalize()
 
-# This should be called at the end
-conveks.mma.finalize()
 
-- 
GitLab