Changeset 21449

Timestamp:

12/22/16 10:39:50 (8 years ago)

Author:

Mathieu Morlighem

Message:

CHG: keep track of all cost functions

File:

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

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

@@ -26 +26 @@
 void simul(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs);
 
+/*Use struct to provide arguments*/
+typedef struct{
+        FemModel  * femmodel;
+        IssmDouble* Jlist;
+        int         M;
+        int         N;
+        int*        i;
+} m1qn3_struct;
+
 void controlm1qn3_core(FemModel* femmodel){
 
@@ -32 +41 @@
         double       f,dxmin,gttol; 
         int          maxsteps,maxiter;
-        int          intn,numberofvertices,num_controls,solution_type;
+        int          intn,numberofvertices,num_controls,num_cost_functions,solution_type;
         IssmDouble  *scaling_factors = NULL;
         IssmDouble  *X  = NULL;
@@ -40 +49 @@
         femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
         femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
+        femmodel->parameters->FindParam(&num_cost_functions,InversionNumCostFunctionsEnum);
         femmodel->parameters->FindParam(&maxsteps,InversionMaxstepsEnum);
         femmodel->parameters->FindParam(&maxiter,InversionMaxiterEnum);
@@ -86 +96 @@
         }
 
-        /*Allocate m1qn3 working arrays (see doc)*/
+        /*Allocate m1qn3 working arrays (see documentation)*/
         long      m   = 100;
         long      ndz = 4*n+m*(2*n+1);
@@ -96 +106 @@
         _printf0_("____________________________________________________________________\n");
 
+        /*Prepare structure for m1qn3*/
+        m1qn3_struct mystruct;
+        mystruct.femmodel = femmodel;
+        mystruct.M        = maxiter;
+        mystruct.N        = num_cost_functions+1;
+        mystruct.Jlist    = xNewZeroInit<IssmDouble>(mystruct.M*mystruct.N);
+        mystruct.i        = xNewZeroInit<int>(1);
+
         /*Initialize Gradient and cost function of M1QN3*/
-        indic = 4; //adjoint and gradient required
-        simul(&indic,&n,X,&f,G,izs,rzs,(void*)femmodel);
+        indic = 4; /*gradient required*/
+        simul(&indic,&n,X,&f,G,izs,rzs,(void*)&mystruct);
 
         /*Estimation of the expected decrease in f during the first iteration*/
@@ -107 +125 @@
                                 &n,X,&f,G,&dxmin,&df1,
                                 &gttol,normtype,&impres,&io,imode,&omode,&niter,&nsim,iz,dz,&ndz,
-                                &reverse,&indic,izs,rzs,(void*)femmodel);
+                                &reverse,&indic,izs,rzs,(void*)&mystruct);
 
         switch(int(omode)){
@@ -137 +155 @@
         ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,G);
         femmodel->OutputControlsx(&femmodel->results);
-        femmodel->results->AddObject(new GenericExternalResult<double>(femmodel->results->Size()+1,JEnum,f));
+        femmodel->results->AddObject(new GenericExternalResult<double*>(femmodel->results->Size()+1,JEnum,mystruct.Jlist,mystruct.M,mystruct.N,0,0));
 
         /*Finalize*/
@@ -153 +171 @@
         xDelete<double>(XL);
         xDelete<double>(scaling_factors);
+        xDelete<double>(mystruct.Jlist);
+        xDelete<int>(mystruct.i);
 }
 
@@ -158 +178 @@
 void simul(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs){
 
-        /*Recover Femmodel*/
-        int         solution_type;
-        FemModel   *femmodel  = (FemModel*)dzs;
-
-        /*Recover number of cost functions responses*/
-        int num_responses,num_controls,numberofvertices;
+        /*Recover Arguments*/
+        m1qn3_struct *input_struct = (m1qn3_struct*)dzs;
+        FemModel     *femmodel     = input_struct->femmodel;
+        IssmDouble   *Jlist        = input_struct->Jlist;
+        int           JlistM       = input_struct->M;
+        int           JlistN       = input_struct->N;
+        int          *Jlisti       = input_struct->i;
+
+        /*Recover some parameters*/
+        int num_responses,num_controls,numberofvertices,solution_type;
         IssmDouble* scaling_factors = NULL;
         femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
         femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
         femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum);
+        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
         numberofvertices=femmodel->vertices->NumberOfVertices();
 
-        /*Constrain input vector*/
+        /*Constrain input vector and update controls*/
         IssmDouble  *XL = NULL;
         IssmDouble  *XU = NULL;
@@ -183 +208 @@
                 }
         }
-
-        /*Update control input*/
         SetControlInputsFromVectorx(femmodel,X);
-
-        /*Recover some parameters*/
-        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
 
         /*Compute solution and adjoint*/
@@ -196 +216 @@
         solutioncore(femmodel);
 
+        /*Check size of Jlist to avoid crashes*/
+        _assert_((*Jlisti)<JlistM);
+        _assert_(JlistN==num_responses+1);
+
         /*Compute objective function*/
-        IssmDouble* Jlist = NULL;
-        femmodel->CostFunctionx(pf,&Jlist,NULL);
+        IssmDouble* Jtemp = NULL;
+        femmodel->CostFunctionx(pf,&Jtemp,NULL);
         _printf0_("f(x) = "<<setw(12)<<setprecision(7)<<*pf<<"  |  ");
+
+        /*Record cost function values and delete Jtemp*/
+        for(int i=0;i<num_responses;i++) Jlist[(*Jlisti)*JlistN+i] = Jtemp[i];
+        Jlist[(*Jlisti)*JlistN+num_responses] = *pf;
+        xDelete<IssmDouble>(Jtemp);
 
         if(*indic==0){
@@ -206 +235 @@
                 /*Retrieve objective functions independently*/
                 _printf0_("            N/A |\n");
-                for(int i=0;i<num_responses;i++) _printf0_(" "<<setw(12)<<setprecision(7)<<Jlist[i]);
+                for(int i=0;i<num_responses;i++) _printf0_(" "<<setw(12)<<setprecision(7)<<Jlist[(*Jlisti)*JlistN+i]);
                 _printf0_("\n");
 
-                xDelete<IssmDouble>(Jlist);
+                *Jlisti = (*Jlisti) +1;
                 xDelete<IssmDouble>(XU);
                 xDelete<IssmDouble>(XL);
@@ -241 +270 @@
         /*Print info*/
         _printf0_("       "<<setw(12)<<setprecision(7)<<Gnorm<<" |");
-        for(int i=0;i<num_responses;i++) _printf0_(" "<<setw(12)<<setprecision(7)<<Jlist[i]);
+        for(int i=0;i<num_responses;i++) _printf0_(" "<<setw(12)<<setprecision(7)<<Jlist[(*Jlisti)*JlistN+i]);
         _printf0_("\n");
 
         /*Clean-up and return*/
-        xDelete<IssmDouble>(Jlist);
+        *Jlisti = (*Jlisti) +1;
         xDelete<IssmDouble>(XU);
         xDelete<IssmDouble>(XL);