Changeset 18714

Timestamp:

10/30/14 15:04:40 (10 years ago)

Author:

Mathieu Morlighem

Message:

NEW: enable multi-parametric optimization with m1qn3

File:

• issm/trunk-jpl/src/c/cores/controlm1qn3_core.cpp (modified) (9 diffs)

issm/trunk-jpl/src/c/cores/controlm1qn3_core.cpp

@@ -32 +32 @@
         double       f,dxmin,gttol; 
         int          maxsteps,maxiter;
-        int          intn,num_controls,solution_type;
+        int          intn,numberofvertices,num_controls,solution_type;
         IssmDouble  *scaling_factors = NULL;
         IssmDouble  *X  = NULL;
@@ -46 +46 @@
         femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum);
         femmodel->parameters->SetParam(false,SaveResultsEnum);
+        numberofvertices=femmodel->vertices->NumberOfVertices();
 
         /*Initialize M1QN3 parameters*/
@@ -62 +63 @@
 
         /*Optimization criterions*/
-        long      niter        = long(maxsteps); /*Maximum number of iterations*/
-        long      nsim         = long(maxiter);/*Maximum number of function calls*/
+        long niter = long(maxsteps); /*Maximum number of iterations*/
+        long nsim  = long(maxiter);/*Maximum number of function calls*/
 
         /*Get initial guess*/
@@ -71 +72 @@
         Xpetsc->GetSize(&intn);
         delete Xpetsc;
+        _assert_(intn==numberofvertices*num_controls);
 
         /*Get problem dimension and initialize gradient and initial guess*/
@@ -77 +79 @@
 
         /*Scale control for M1QN3*/
-        if(num_controls!=1) _error_("not supported yet...");
-        for(long i=0;i<n;i++){
-                X[i] = X[i]/scaling_factors[0];
+        for(int i=0;i<numberofvertices;i++){
+                for(int c=0;c<num_controls;c++){
+                        int index = num_controls*i+c;
+                        X[index] = X[index]/scaling_factors[c];
+                }
         }
 
@@ -118 +122 @@
         GetVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
         GetVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
-        if(num_controls!=1) _error_("not supported yet");
-        for(long i=0;i<n;i++){
-                X[i] = X[i]*scaling_factors[0];
-                if(X[i]>XU[i]) X[i]=XU[i];
-                if(X[i]<XL[i]) X[i]=XL[i];
+        for(int i=0;i<numberofvertices;i++){
+                for(int c=0;c<num_controls;c++){
+                        int index = num_controls*i+c;
+                        X[index] = X[index]*scaling_factors[c];
+                        if(X[index]>XU[index]) X[index]=XU[index];
+                        if(X[index]<XL[index]) X[index]=XL[index];
+                }
         }
         SetControlInputsFromVectorx(femmodel,X);
@@ -152 +158 @@
 
         /*Recover number of cost functions responses*/
-        int         num_responses;
-        int         num_controls;
+        int num_responses,num_controls,numberofvertices;
         IssmDouble* scaling_factors = NULL;
         femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
         femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
         femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum);
+        numberofvertices=femmodel->vertices->NumberOfVertices();
 
         /*Constrain input vector*/
@@ -164 +170 @@
         GetVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
         GetVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
-        if(num_controls!=1) _error_("not supported yet");
-        for(long i=0;i<*n;i++){
-                X[i] = X[i]*scaling_factors[0];
-                if(X[i]>XU[i]) X[i]=XU[i];
-                if(X[i]<XL[i]) X[i]=XL[i];
+        for(int i=0;i<numberofvertices;i++){
+                for(int c=0;c<num_controls;c++){
+                        int index = num_controls*i+c;
+                        X[index] = X[index]*scaling_factors[c];
+                        if(X[index]>XU[index]) X[index]=XU[index];
+                        if(X[index]<XL[index]) X[index]=XL[index];
+                }
         }
 
@@ -214 +222 @@
         /*Constrain Gradient*/
         IssmDouble  Gnorm = 0.;
-        if(num_controls!=1) _error_("not supported yet");
-        for(long i=0;i<*n;i++){
-                if(X[i]>=XU[i]) G[i]=0.;
-                if(X[i]<=XL[i]) G[i]=0.;
-                G[i] = G[i]*scaling_factors[0];
-                X[i] = X[i]/scaling_factors[0];
-                Gnorm += G[i]*G[i];
+        for(int i=0;i<numberofvertices;i++){
+                for(int c=0;c<num_controls;c++){
+                        int index = num_controls*i+c;
+                        if(X[index]>=XU[index]) G[index]=0.;
+                        if(X[index]<=XL[index]) G[index]=0.;
+                        G[index] = G[index]*scaling_factors[c];
+                        X[index] = X[index]/scaling_factors[c];
+                        Gnorm += G[index]*G[index];
+                }
         }
         Gnorm = sqrt(Gnorm);