Changeset 22829

Timestamp:

06/06/18 12:03:23 (7 years ago)

Author:

erobo

Message:

CHG: allowing AD M1QN3 to remember the 'best' paramters

File:

• issm/trunk-jpl/src/c/cores/controladm1qn3_core.cpp (modified) (8 diffs)

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

@@ -30 +30 @@
         int          N;
         int*         i;
+        IssmPDouble* X_best;
+        IssmPDouble* G_best;
+        IssmPDouble* J_best;
 } m1qn3_struct;
 
@@ -77 +80 @@
         }
 
-        IssmPDouble  *Jlist        = input_struct->Jlist;
+        IssmPDouble*  Jlist        = input_struct->Jlist;
         int           JlistM       = input_struct->M;
         int           JlistN       = input_struct->N;
-        int          *Jlisti       = input_struct->i;
+        int*          Jlisti       = input_struct->i;
+        IssmPDouble*  J_best       = input_struct->J_best;
+        IssmPDouble*  X_best                    = input_struct->X_best;
+        IssmPDouble*  G_best                    = input_struct->G_best;
         int           intn         = (int)*n;
 
@@ -292 +298 @@
         Jlist[(*Jlisti)*JlistN+num_responses] = reCast<IssmPDouble>(J);
 
+        
+        if(*J_best<0 || J<*J_best){
+                *J_best = reCast<IssmPDouble>(J);
+                for(int i=0;i<intn;i++){
+                        X_best[i] = reCast<IssmPDouble>(X[i]);
+                        G_best[i] = reCast<IssmPDouble>(G[i]);
+                }
+        }
+
+/*
+        IssmDouble* test = xNew<IssmDouble>(intn);
+        femmodel->parameters->FindParam(&test,NULL,InversionXBestEnum);
+        for(int i=0;i<10;i++)_printf_("X "<<test[i]<<"\n");
+        xDelete<IssmDouble>(intn);
+*/
         if(*indic==0){
                 /*dry run, no gradient required*/
@@ -413 +434 @@
         long      ndz = 4*n+m*(2*n+1);
         double*   dz  = xNew<double>(ndz);
-
+        IssmDouble J_best = -10.;
         if(VerboseControl())_printf0_("   Computing initial solution\n");
         _printf0_("\n");
@@ -426 +447 @@
         mystruct.Jlist    = xNewZeroInit<IssmPDouble>(mystruct.M*mystruct.N);
         mystruct.i        = xNewZeroInit<int>(1);
+        mystruct.J_best   = xNewZeroInit<IssmPDouble>(1);
+        mystruct.X_best = xNewZeroInit<IssmPDouble>(intn);
+        mystruct.G_best   = xNewZeroInit<IssmPDouble>(intn);
         /*Initialize Gradient and cost function of M1QN3*/
         indic = 4; /*gradient required*/
@@ -465 +489 @@
                 N_add +=N[c];
         }
-                
+        
+        IssmDouble* aX_best = NULL;
+        IssmDouble* aG_best = NULL;
+
+        femmodel->parameters->FindParam(&aX_best,NULL,InversionXBestEnum);
+        femmodel->parameters->FindParam(&aG_best,NULL,InversionGBestEnum);
+        
         /*Set X as our new control*/
         IssmDouble* aX=xNew<IssmDouble>(intn);
+        IssmDouble* aG=xNew<IssmDouble>(intn);
+        double* X_best = xNew<double>(intn);
+        double* G_best = xNew<double>(intn);
+
         for(int i=0;i<intn;i++) {
                 aX[i] = reCast<IssmDouble>(X[i]); 
-                }
+                aG[i] = reCast<IssmDouble>(G[i]);
+                X_best[i] = reCast<double>(aX_best[i]); 
+                G_best[i] = reCast<double>(aG_best[i]); 
+                }
+
+        ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,aG);
         SetControlInputsFromVectorx(femmodel,aX);
+        
         xDelete(aX);
 
-        /*Set final gradient in inputs*/
-        IssmDouble* aG=xNew<IssmDouble>(intn);
-        for(int i=0;i<intn;i++) {
-                aG[i] = reCast<IssmDouble>(G[i]);
-        }
-        ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,aG);
 
         if (solution_type == TransientSolutionEnum){
@@ -493 +527 @@
                         GenericExternalResult<IssmPDouble*>* X_output = new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,control_enum[i],&X[offset],N[i],numberofvertices,1,0.);
 
+                        GenericExternalResult<IssmPDouble*>* Gbest_output = new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,Outputdefinition90Enum+i,&G_best[offset],N[i],numberofvertices,1,0.);
+                        GenericExternalResult<IssmPDouble*>* Xbest_output = new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,Outputdefinition80Enum+i,&X_best[offset],N[i],numberofvertices,1,0.);
+
                         /*transpose for consistency with MATLAB's formating*/
                         G_output->Transpose();
                         X_output->Transpose();
+                        Gbest_output->Transpose();
+                        Xbest_output->Transpose();
 
                         /*Add to results*/
                         femmodel->results->AddObject(G_output);
                         femmodel->results->AddObject(X_output);
-
+                        femmodel->results->AddObject(Gbest_output);
+                        femmodel->results->AddObject(Xbest_output);
+                        
                         offset += N[i]*numberofvertices;
                 }
@@ -524 +565 @@
         xDelete<double>(G);
         xDelete<double>(X);
+        xDelete<double>(X_best);
+        xDelete<double>(G_best);
         xDelete<double>(dz);
         xDelete<double>(XU);