Changeset 18942

Timestamp:

12/04/14 15:43:24 (10 years ago)

Author:

Eric.Larour

Message:

CHG: rehooking ad gradient core

File:

• issm/trunk-jpl/src/c/cores/controlad_core.cpp (modified) (10 diffs)

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

@@ -25 +25 @@
 /*Cost function prototype*/
 void simulad(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs);
+FemModel* presimulad(int* pintn, double** pX, FemModel* femmodel);
+void postsimulad(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs);
 /*}}}*/
 void controlad_core(FemModel* femmodel){ /*{{{*/
 
         /*Intermediaries*/
+        FemModel*    femmodelad=NULL;
         int          i;
         long         omode;
@@ -34 +37 @@
         IssmDouble   dxmind,gttold; 
         int          maxsteps,maxiter;
-        int          intn,numberofvertices,num_controls,solution_type;
-        IssmDouble  *scaling_factors = NULL;
+        int          intn,solution_type;
         IssmPDouble  *X  = NULL;
         IssmDouble   *Xd  = NULL;
-        IssmDouble   *Gd  = NULL;
         IssmPDouble  *G  = NULL;
-        bool onsid=true;
 
         /*Recover some parameters*/
         femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
-        femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
         femmodel->parameters->FindParam(&maxsteps,InversionMaxstepsEnum);
         femmodel->parameters->FindParam(&maxiter,InversionMaxiterEnum);
         femmodel->parameters->FindParam(&dxmind,InversionDxminEnum); dxmin=reCast<IssmPDouble>(dxmind);
         femmodel->parameters->FindParam(&gttold,InversionGttolEnum); gttol=reCast<IssmPDouble>(gttold);
-        femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum);
         femmodel->parameters->SetParam(false,SaveResultsEnum);
-        numberofvertices=femmodel->vertices->NumberOfVertices();
 
         /*Initialize M1QN3 parameters*/
@@ -71 +68 @@
         long nsim  = long(maxiter);/*Maximum number of function calls*/
 
-        /*Get initial guess*/
-        Vector<IssmDouble> *Xad = NULL;
-        GetVectorFromControlInputsx(&Xad,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value",onsid);
-        Xd = Xad->ToMPISerial();
-        Xad->GetSize(&intn);
-        X=xNew<IssmPDouble>(intn);
-        for(i=0;i<intn;i++) X[i]=reCast<IssmPDouble>(Xd[i]);
-        delete Xad;
-        _assert_(intn==numberofvertices*num_controls);
-        
-        /*Get problem dimension and initialize gradient and initial guess*/
+        /*Run the first part of simulad, in order to get things started!:*/
+        femmodelad=presimulad(&intn,&X,femmodel);
+
+        /*Get problem dimension and initialize gradient: */
         long n = long(intn);
         G = xNew<IssmPDouble>(n);
-        Gd = xNew<IssmDouble>(n);
-
-        /*Scale control for M1QN3*/
-        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]/reCast<IssmPDouble>(scaling_factors[c]);
-                }
-        }
 
         /*Allocate m1qn3 working arrays (see doc)*/
@@ -104 +85 @@
         _printf0_("____________________________________________________________________\n");
 
-        //first run before firing up the control optimization
-        simulad(&indic,&n,X,&f,G,izs,rzs,(void*)femmodel);
+        //run post simular phase, to fire up the control optimization
+        postsimulad(&indic,&n,X,&f,G,izs,rzs,(void*)femmodelad); 
         double f1=f;
 
@@ -125 +106 @@
         }
         
-        /*Constrain solution vector*/
-        IssmDouble  *XL = NULL;
-        IssmDouble  *XU = NULL;
-        GetVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound",onsid);
-        GetVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound",onsid);
-        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]*reCast<IssmPDouble>(scaling_factors[c]);
-                        if(X[index]>XU[index]) X[index]=reCast<IssmPDouble>(XU[index]);
-                        if(X[index]<XL[index]) X[index]=reCast<IssmPDouble>(XL[index]);
-                }
-        }
-
         /*Save results:*/
         femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,AutodiffJacobianEnum,G,n,1,1,0.0));
@@ -147 +114 @@
         xDelete<double>(X);
         xDelete<double>(dz);
-        xDelete<IssmDouble>(scaling_factors);
 }/*}}}*/
-void simulad(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs){ /*{{{*/
+FemModel* presimulad(int* pintn, double** pX, FemModel* femmodel){ /*{{{*/
 
         /*Intermediaries:*/
@@ -157 +123 @@
         char* toolkitsfilename=NULL;
         char* lockfilename=NULL;
-        IssmPDouble* G2=NULL;
-        bool onsid=true;
-        IssmDouble  *XL = NULL;
-        IssmDouble  *XU = NULL;
         int         solution_type;
-        FemModel   *femmodel  =  NULL;
-        FemModel   *femmodelad  = NULL;
         IssmDouble    pfd;
+        IssmDouble*   Xd=NULL;
+        int           intn;
+        IssmPDouble*   X=NULL;
         int            i;
         
-        /*Recover Femmodel*/
-        femmodel  = (FemModel*)dzs;
-
-        /*Recover number of cost functions responses*/
-        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*/
-        GetVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound",onsid);
-        GetVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound",onsid);
-        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]*reCast<IssmPDouble>(scaling_factors[c]);
-                        if(X[index]>reCast<IssmPDouble>(XU[index])) X[index]=reCast<IssmPDouble>(XU[index]);
-                        if(X[index]<reCast<IssmPDouble>(XL[index])) X[index]=reCast<IssmPDouble>(XL[index]);
-                }
-        }
-
         /*Now things get complicated. The femmodel we recovered did not initialize an AD trace, so we can't compute gradients with it. We are going to recreate 
          *a new femmodel, identical in all aspects to the first one, with trace on though, which will allow us to run the forward mode and get the gradient 
@@ -199 +139 @@
         femmodel->parameters->FindParam(&lockfilename,LockFileNameEnum);
 
-        femmodelad=new FemModel(rootpath, inputfilename, outputfilename, toolkitsfilename, lockfilename, femmodel->comm, femmodel->solution_type,X);
-        femmodel=femmodelad; //We can do this, because femmodel is being called from outside, not by reference, so we won't erase it
-        
+        femmodel=new FemModel(rootpath, inputfilename, outputfilename, toolkitsfilename, lockfilename, femmodel->comm, femmodel->solution_type,NULL);
+
+        
+        /*Get initial guess:*/
+        femmodel->parameters->FindParam(&Xd,&intn,AutodiffXpEnum);
+        X=xNew<IssmPDouble>(intn); for(i=0;i<intn;i++) X[i]=reCast<IssmPDouble>(Xd[i]); 
+
+        xDelete<char>(rootpath);
+        xDelete<char>(inputfilename);
+        xDelete<char>(outputfilename);
+        xDelete<char>(toolkitsfilename);
+        xDelete<char>(lockfilename);
+        xDelete<IssmDouble>(Xd);
+
+        *pintn=intn;
+        *pX=X;
+
+        return femmodel;
+
+} /*}}}*/
+void postsimulad(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs){ /*{{{*/
+
+        /*Intermediaries:*/
+        char* rootpath=NULL;
+        char* inputfilename=NULL;
+        char* outputfilename=NULL;
+        char* toolkitsfilename=NULL;
+        char* lockfilename=NULL;
+        IssmPDouble* G2=NULL;
+        int         solution_type;
+        FemModel   *femmodel  =  NULL;
+        IssmDouble    pfd;
+        int            i;
+        
+        /*Recover Femmodel*/
+        femmodel  = (FemModel*)dzs;
+
+        /*Recover number of cost functions responses*/
+        int num_responses;
+        femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
+
         /*Recover some parameters*/
         femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
@@ -236 +214 @@
         /*Constrain X and G*/
         IssmDouble  Gnorm = 0.;
-        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]*reCast<IssmPDouble>(scaling_factors[c]);
-                        X[index] = X[index]/reCast<IssmPDouble>(scaling_factors[c]);
-                        Gnorm += G[index]*G[index];
-                }
-        }
+        for(long i=0;i<*n;i++) Gnorm += G[i]*G[i];
         Gnorm = sqrt(Gnorm);
 
@@ -255 +224 @@
         /*Clean-up and return*/
         xDelete<IssmDouble>(Jlist);
-        xDelete<IssmDouble>(XU);
-        xDelete<IssmDouble>(XL);
         xDelete<IssmPDouble>(G2);
-        xDelete<IssmDouble>(scaling_factors);
+        
+        xDelete<char>(rootpath);
+        xDelete<char>(inputfilename);
+        xDelete<char>(outputfilename);
+        xDelete<char>(toolkitsfilename);
+        xDelete<char>(lockfilename);
+
+} /*}}}*/
+void simulad(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs){ /*{{{*/
+
+        /*Intermediaries:*/
+        char* rootpath=NULL;
+        char* inputfilename=NULL;
+        char* outputfilename=NULL;
+        char* toolkitsfilename=NULL;
+        char* lockfilename=NULL;
+        IssmPDouble* G2=NULL;
+        int         solution_type;
+        FemModel   *femmodel  =  NULL;
+        FemModel   *femmodelad  = NULL;
+        IssmDouble    pfd;
+        int            i;
+        
+        /*Recover Femmodel*/
+        femmodel  = (FemModel*)dzs;
+
+        /*Recover number of cost functions responses*/
+        int num_responses;
+        femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
+
+        /*Now things get complicated. The femmodel we recovered did not initialize an AD trace, so we can't compute gradients with it. We are going to recreate 
+         *a new femmodel, identical in all aspects to the first one, with trace on though, which will allow us to run the forward mode and get the gradient 
+         in one run of the solution core. So first recover the filenames required for the FemModel constructor, then call a new ad tailored constructor:*/
+        femmodel->parameters->FindParam(&rootpath,RootPathEnum);
+        femmodel->parameters->FindParam(&inputfilename,InputFileNameEnum);
+        femmodel->parameters->FindParam(&outputfilename,OutputFileNameEnum);
+        femmodel->parameters->FindParam(&toolkitsfilename,ToolkitsFileNameEnum);
+        femmodel->parameters->FindParam(&lockfilename,LockFileNameEnum);
+
+        femmodelad=new FemModel(rootpath, inputfilename, outputfilename, toolkitsfilename, lockfilename, femmodel->comm, femmodel->solution_type,X);
+        femmodel=femmodelad; //We can do this, because femmodel is being called from outside, not by reference, so we won't erase it
+        
+        /*Recover some parameters*/
+        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
+
+        /*Compute solution:*/
+        void (*solutioncore)(FemModel*)=NULL;
+        CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type);
+        solutioncore(femmodel);
+
+        /*Compute objective function*/
+        IssmDouble* Jlist = NULL;
+        femmodel->CostFunctionx(&pfd,&Jlist,NULL); *pf=reCast<IssmPDouble>(pfd);
+        _printf0_("f(x) = "<<setw(12)<<setprecision(7)<<*pf<<"  |  ");
+        
+        /*Compute gradient using AD. Gradient is in the results after the ad_core is called*/
+        adgradient_core(femmodel); 
+
+        if(IssmComm::GetRank()==0){
+                IssmPDouble* G_temp=NULL;
+                GenericExternalResult<IssmPDouble*>* gradient=(GenericExternalResult<IssmPDouble*>*)femmodel->results->FindResult(AutodiffJacobianEnum); _assert_(gradient);
+                G_temp=gradient->GetValues();
+                /*copy onto G2, to avoid a leak: */
+                G2=xNew<IssmPDouble>(*n); 
+                xMemCpy<IssmPDouble>(G2,G_temp,*n);
+        }
+        else G2=xNew<IssmPDouble>(*n);
+
+        /*MPI broadcast results:*/
+        ISSM_MPI_Bcast(G2,*n,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
+        
+        /*Send gradient to m1qn3 core: */
+        for(long i=0;i<*n;i++) G[i] = G2[i];
+        
+        /*Recover Gnorm: */
+        IssmDouble  Gnorm = 0.;
+        for(int i=0;i<*n;i++) Gnorm += G[i]*G[i];
+        Gnorm = sqrt(Gnorm);
+
+        /*Print info*/
+        _printf0_("       "<<setw(12)<<setprecision(7)<<Gnorm<<" |");
+        for(int i=0;i<num_responses;i++) _printf0_(" "<<setw(12)<<setprecision(7)<<Jlist[i]);
+        _printf0_("\n");
+
+        /*Clean-up and return*/
+        xDelete<IssmDouble>(Jlist);
+        xDelete<IssmPDouble>(G2);
         
         xDelete<char>(rootpath);