source: issm/oecreview/Archive/18296-19100/ISSM-18615-18616.diff

../trunk-jpl/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp

@@ -12 +12 @@
         bool        control_analysis;
         int         inversiontype;
         int         nsteps;
-        int         num_control_type;
-        int         num_cm_responses;
+        int         num_controls;
+        int         num_costfunc;
         int        *control_type     = NULL;
         int        *maxiter          = NULL;
         int        *cm_responses     = NULL;
         IssmDouble *cm_jump          = NULL;
         IssmDouble *optscal          = NULL;
+        IssmDouble *control_scaling_factors = NULL;
 
         /*retrieve some parameters: */
         iomodel->Constant(&control_analysis,InversionIscontrolEnum);
@@ -40 +41 @@
                 }
 
                 /*Now, recover fit, optscal and maxiter as vectors: */
-                iomodel->FetchData(&control_type,NULL,&num_control_type,InversionControlParametersEnum);
-                iomodel->FetchData(&cm_responses,NULL,&num_cm_responses,InversionCostFunctionsEnum);
-                parameters->AddObject(new IntVecParam(InversionControlParametersEnum,control_type,num_control_type));
-                parameters->AddObject(new IntVecParam(InversionCostFunctionsEnum,cm_responses,num_cm_responses));
+                iomodel->FetchData(&control_type,NULL,&num_controls,InversionControlParametersEnum);
+                iomodel->FetchData(&cm_responses,NULL,&num_costfunc,InversionCostFunctionsEnum);
+                parameters->AddObject(new IntVecParam(InversionControlParametersEnum,control_type,num_controls));
+                parameters->AddObject(new IntVecParam(InversionCostFunctionsEnum,cm_responses,num_costfunc));
 
                 /*Inversion type specifics*/
                 switch(inversiontype){
@@ -53 +54 @@
                                 iomodel->FetchData(&cm_jump,&nsteps,NULL,InversionStepThresholdEnum);
                                 iomodel->FetchData(&optscal,NULL,NULL,InversionGradientScalingEnum);
                                 iomodel->FetchData(&maxiter,NULL,NULL,InversionMaxiterPerStepEnum);
-                                parameters->AddObject(new DoubleMatParam(InversionGradientScalingEnum,optscal,nsteps,num_control_type));
+                                parameters->AddObject(new DoubleMatParam(InversionGradientScalingEnum,optscal,nsteps,num_controls));
                                 parameters->AddObject(new DoubleVecParam(InversionStepThresholdEnum,cm_jump,nsteps));
                                 parameters->AddObject(new IntVecParam(InversionMaxiterPerStepEnum,maxiter,nsteps));
                                 break;
@@ -72 +73 @@
                                 parameters->AddObject(iomodel->CopyConstantObject(InversionGttolEnum));
                                 parameters->AddObject(iomodel->CopyConstantObject(InversionMaxstepsEnum));
                                 parameters->AddObject(iomodel->CopyConstantObject(InversionMaxiterEnum));
+                                iomodel->FetchData(&control_scaling_factors,NULL,NULL,InversionControlScalingFactorsEnum);
+                                parameters->AddObject(new DoubleVecParam(InversionControlScalingFactorsEnum,control_scaling_factors,num_controls));
                                 break;
                         case 3:/*Validation*/
                                 break;

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

@@ -32 +32 @@
         double       f,dxmin,gttol; 
         int          maxsteps,maxiter;
         int          intn,num_controls,solution_type;
+        IssmDouble  *scaling_factors = NULL;
         IssmDouble  *X  = NULL;
         IssmDouble  *G  = NULL;
 
@@ -42 +43 @@
         femmodel->parameters->FindParam(&maxiter,InversionMaxiterEnum);
         femmodel->parameters->FindParam(&dxmin,InversionDxminEnum);
         femmodel->parameters->FindParam(&gttol,InversionGttolEnum);
+        femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum);
         femmodel->parameters->SetParam(false,SaveResultsEnum);
 
         /*Initialize M1QN3 parameters*/
@@ -73 +75 @@
         long n = long(intn);
         G = xNew<double>(n);
 
+        /*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];
+        }
+
         /*Allocate m1qn3 working arrays (see doc)*/
         long      m   = 100;
         long      ndz = 4*n+m*(2*n+1);
@@ -131 +139 @@
         FemModel   *femmodel  = (FemModel*)dzs;
 
         /*Recover number of cost functions responses*/
-        int num_responses;
+        int         num_responses;
+        int         num_controls;
+        IssmDouble* scaling_factors = NULL;
         femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
+        femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
+        femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum);
 
         /*Constrain input vector*/
         IssmDouble  *XL = NULL;
         IssmDouble  *XU = NULL;
         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];
         }
@@ -161 +175 @@
         femmodel->CostFunctionx(pf,&Jlist,NULL);
         _printf0_("f(x) = "<<setw(12)<<setprecision(7)<<*pf<<"  |  ");
 
-
-
         if(indic==0){
                 /*dry run, no gradient required*/
 
@@ -189 +201 @@
 
         /*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];
         }
         Gnorm = sqrt(Gnorm);
@@ -202 +217 @@
         _printf0_("\n");
 
         /*Clean-up and return*/
-                xDelete<IssmDouble>(Jlist);
+        xDelete<IssmDouble>(Jlist);
         xDelete<IssmDouble>(XU);
         xDelete<IssmDouble>(XL);
 }