Changeset 18128

Timestamp:

06/10/14 10:18:16 (11 years ago)

Author:

Mathieu Morlighem

Message:

NEW: working on BrentSearch interface to make it tao like

Location:

issm/trunk-jpl/src

Files:

• c/classes/Elements/Penta.cpp (modified) (4 diffs)
• c/classes/Elements/Tria.cpp (modified) (1 diff)
• c/cores/control_core.cpp (modified) (8 diffs)
• c/cores/controlm1qn3_core.cpp (modified) (5 diffs)
• c/cores/controltao_core.cpp (modified) (3 diffs)
• c/cores/controlvalidation_core.cpp (modified) (1 diff)
• c/modules/Gradjx/Gradjx.cpp (modified) (1 diff)
• c/modules/ModelProcessorx/Control/CreateParametersControl.cpp (modified) (3 diffs)
• c/modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.cpp (modified) (3 diffs)
• c/modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.h (modified) (1 diff)
• c/shared/Numerics/BrentSearch.cpp (modified) (2 diffs)
• c/shared/Numerics/OptPars.h (modified) (1 diff)
• c/shared/Numerics/numerics.h (modified) (1 diff)
• m/classes/inversion.m (modified) (2 diffs)

issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

@@ -2851 +2851 @@
 
         /*Get input (either in element or material)*/
+        if(control_enum==MaterialsRheologyBbarEnum) control_enum=MaterialsRheologyBEnum;
         Input* input=inputs->GetInput(control_enum);
         if(!input) _error_("Input " << EnumToStringx(control_enum) << " not found in element");
@@ -2865 +2866 @@
 
         IssmDouble  values[NUMVERTICES];
-        int     vertexpidlist[NUMVERTICES];
-        Input  *input     = NULL;
-        Input  *new_input = NULL;
+        int         vertexpidlist[NUMVERTICES],control_init;
+        Input      *input     = NULL;
+        Input      *new_input = NULL;
+
+        /*Specific case for depth averaged quantities*/
+        control_init=control_enum;
+        if(control_enum==MaterialsRheologyBbarEnum){
+                control_enum=MaterialsRheologyBEnum;
+                if(!IsOnBase()) return;
+        }
+        if(control_enum==DamageDbarEnum){
+                control_enum=DamageDEnum;
+                if(!IsOnBase()) return;
+        }
 
         /*Get out if this is not an element input*/
@@ -2876 +2888 @@
 
         /*Get values on vertices*/
-        for (int i=0;i<NUMVERTICES;i++){
+        for(int i=0;i<NUMVERTICES;i++){
                 values[i]=vector[vertexpidlist[i]];
         }
@@ -2887 +2899 @@
 
         ((ControlInput*)input)->SetInput(new_input);
+
+        if(control_init==MaterialsRheologyBbarEnum){
+                this->InputExtrude(control_enum);
+        }
+        if(control_init==DamageDbarEnum){
+                this->InputExtrude(control_enum);
+        }
 }
 /*}}}*/

issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

@@ -983 +983 @@
                 this->inputs->AddInput(new IntInput(ApproximationEnum,reCast<int>(iomodel->Data(FlowequationElementEquationEnum)[index])));
         }
-
 
         /*Control Inputs*/

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

@@ -11 +11 @@
 
 /*Local prototypes*/
-bool controlconvergence(IssmDouble J, IssmDouble tol_cm);
-IssmDouble objectivefunction(IssmDouble search_scalar,void* optargs);
+IssmDouble FormFunction(IssmDouble* X,void* usr);
+IssmDouble FormFunctionGradient(IssmDouble** pG,IssmDouble* X,void* usr);
+typedef struct {
+        FemModel* femmodel;
+        int       nsize;
+} AppCtx;
 
 void control_core(FemModel* femmodel){
@@ -19 +23 @@
 
         /*parameters: */
-        int        num_controls;
+        int        num_controls,nsize;
         int        nsteps;
         IssmDouble tol_cm;
@@ -27 +31 @@
 
         int        *control_type   = NULL;
-        IssmDouble *maxiter        = NULL;
+        int*        maxiter        = NULL;
         IssmDouble *cm_jump        = NULL;
 
         /*intermediary: */
-        IssmDouble search_scalar = 1.;
         OptPars    optpars;
 
@@ -61 +64 @@
         if(isFS) solutioncore(femmodel);
 
-        /*Initialize cost function: */
-        J=xNew<IssmDouble>(nsteps);
+        /*Get initial guess*/
+        Vector<IssmDouble> *Xpetsc = NULL;
+        GetVectorFromControlInputsx(&Xpetsc,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value");
+        IssmDouble* X0 = Xpetsc->ToMPISerial();
+        Xpetsc->GetSize(&nsize);
+        delete Xpetsc;
 
         /*Initialize some of the BrentSearch arguments: */
-        optpars.xmin=0; optpars.xmax=1;
-
-        /*Start looping: */
-        for(int n=0;n<nsteps;n++){
-
-                /*Display info*/
-                if(VerboseControl()) _printf0_("\n" << "   control method step " << n+1 << "/" << nsteps << "\n");
-
-
-                /*In steady state inversion, compute new temperature field now*/
-                if(solution_type==SteadystateSolutionEnum) solutioncore(femmodel);
-
-                if(VerboseControl()) _printf0_("   compute adjoint state:\n");
-                adjointcore(femmodel);
-                gradient_core(femmodel,n,search_scalar==0.);
-
-                if(VerboseControl()) _printf0_("   optimizing along gradient direction\n");
-                optpars.maxiter=reCast<int,IssmDouble>(maxiter[n]); optpars.cm_jump=cm_jump[n];
-                BrentSearch(&search_scalar,J+n,&optpars,&objectivefunction,(void*)femmodel);
-
-                if(VerboseControl()) _printf0_("   updating parameter using optimized search scalar\n"); //true means update save controls
-                InputControlUpdatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,search_scalar,true);
-
-                if(controlconvergence(J[n],tol_cm)) break;
-        }
+        optpars.xmin    = 0; 
+        optpars.xmax    = 1;
+        optpars.nsteps  = nsteps;
+        optpars.nsize   = nsize;
+        optpars.maxiter = maxiter;
+        optpars.cm_jump = cm_jump;
+
+        /*Initialize function argument*/
+        AppCtx usr;
+        usr.femmodel = femmodel;
+        usr.nsize    = nsize;
+
+        /*Call Brent optimization*/
+        BrentSearch(&J,optpars,X0,&FormFunction,&FormFunctionGradient,(void*)&usr);
 
         if(VerboseControl()) _printf0_("   preparing final solution\n");
+        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");
+        for(long i=0;i<nsize;i++){
+                if(X0[i]>XU[i]) X0[i]=XU[i];
+                if(X0[i]<XL[i]) X0[i]=XL[i];
+        }
+        xDelete<IssmDouble>(XU);
+        xDelete<IssmDouble>(XL);
+        SetControlInputsFromVectorx(femmodel,X0);
         femmodel->parameters->SetParam(true,SaveResultsEnum);
         solutioncore(femmodel);
@@ -111 +118 @@
         /*Free ressources: */
         xDelete<int>(control_type);
-        xDelete<IssmDouble>(maxiter);
+        xDelete<int>(maxiter);
         xDelete<IssmDouble>(cm_jump);
         xDelete<IssmDouble>(J);
+        xDelete<IssmDouble>(X0);
 }
-bool controlconvergence(IssmDouble J, IssmDouble tol_cm){
-
-        bool converged=false;
-
-        /*Has convergence been reached?*/
-        if (!xIsNan<IssmDouble>(tol_cm) && J<tol_cm){
-                converged=true;
-                if(VerboseConvergence()) _printf0_("      Convergence criterion reached: J = " << J << " < " << tol_cm);
-        }
-
-        return converged;
-}
-
-IssmDouble objectivefunction(IssmDouble search_scalar,void* optargs){
+
+IssmDouble FormFunction(IssmDouble* X,void* usrvoid){
 
         /*output: */
@@ -134 +130 @@
 
         /*parameters: */
-        int        solution_type,analysis_type;
-        bool       isFS       = false;
+        int        solution_type,analysis_type,num_responses;
+        bool       isFS           = false;
         bool       conserve_loads = true;
-        FemModel  *femmodel       = (FemModel*)optargs;
+        AppCtx*    usr = (AppCtx*)usrvoid;
+        FemModel  *femmodel  = usr->femmodel;
+        int        nsize     = usr->nsize;
 
         /*Recover parameters: */
@@ -143 +141 @@
         femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
         femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
-
-        /*set analysis type to compute velocity: */
+        femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
+
+        /*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");
+        for(long i=0;i<nsize;i++){
+                if(X[i]>XU[i]) X[i]=XU[i];
+                if(X[i]<XL[i]) X[i]=XL[i];
+        }
+
+        /*Update control input*/
+        SetControlInputsFromVectorx(femmodel,X);
+
+        /*solve forward: */
         switch(solution_type){
                 case SteadystateSolutionEnum:
+                        femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
+                        stressbalance_core(femmodel);   //We need a 3D velocity!! (vz is required for the next thermal run)
+                        break;
                 case StressbalanceSolutionEnum:
                         femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
+                        solutionsequence_nonlinear(femmodel,conserve_loads); 
                         break;
                 case BalancethicknessSolutionEnum:
                         femmodel->SetCurrentConfiguration(BalancethicknessAnalysisEnum);
+                        solutionsequence_linear(femmodel); 
                         break;
                 case BalancethicknessSoftSolutionEnum:
-                        femmodel->SetCurrentConfiguration(BalancethicknessAnalysisEnum);
+                        /*NOTHING*/
                         break;
                 case Balancethickness2SolutionEnum:
                         femmodel->SetCurrentConfiguration(Balancethickness2AnalysisEnum);
+                        solutionsequence_linear(femmodel); 
                         break;
                 default:
@@ -163 +181 @@
         }
 
-        /*update parameter according to scalar: */ //false means: do not save control
-        InputControlUpdatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,search_scalar,false);
-
-        /*Run stressbalance with updated inputs: */
-        if (solution_type==SteadystateSolutionEnum){
-                stressbalance_core(femmodel);   //We need a 3D velocity!! (vz is required for the next thermal run)
-        }
-        else if (solution_type==StressbalanceSolutionEnum){
-                solutionsequence_nonlinear(femmodel,conserve_loads); 
-        }
-        else if (solution_type==BalancethicknessSolutionEnum){
-                solutionsequence_linear(femmodel); 
-        }
-        else if (solution_type==Balancethickness2SolutionEnum){
-                solutionsequence_linear(femmodel); 
-        }
-        else if (solution_type==BalancethicknessSoftSolutionEnum){
-                /*Don't do anything*/
-        }
-        else{
-                _error_("Solution " << EnumToStringx(solution_type) << " not implemented yet");
-        }
-
         /*Compute misfit for this velocity field.*/
-        femmodel->CostFunctionx(&J,NULL,NULL);
+        IssmDouble* Jlist = NULL;
+        femmodel->CostFunctionx(&J,&Jlist,NULL);
+        //_printf0_("f(x) = "<<setw(12)<<setprecision(7)<<*pf<<"  |  ");
+
+        /*Retrieve objective functions independently*/
+        //for(int i=0;i<num_responses;i++) _printf0_(" "<<setw(12)<<setprecision(7)<<Jlist[i]);
+        //_printf0_("\n");
 
         /*Free ressources:*/
+        xDelete<IssmDouble>(XU);
+        xDelete<IssmDouble>(XL);
+        xDelete<IssmDouble>(Jlist);
         return J;
 }
+IssmDouble FormFunctionGradient(IssmDouble** pG,IssmDouble* X,void* usrvoid){
+
+        /*output: */
+        IssmDouble J;
+
+        /*parameters: */
+        void (*adjointcore)(FemModel*)=NULL;
+        int         solution_type,analysis_type,num_responses,num_controls,numvertices;
+        bool        isFS           = false;
+        bool        conserve_loads = true;
+        IssmDouble *scalar_list    = NULL;
+        IssmDouble *Jlist          = NULL;
+        IssmDouble *G              = NULL;
+        IssmDouble *norm_list      = NULL;
+        AppCtx     *usr            = (AppCtx*)usrvoid;
+        FemModel   *femmodel       = usr->femmodel;
+        int         nsize          = usr->nsize;
+
+        /*Recover parameters: */
+        femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum);
+        femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
+        femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
+        femmodel->parameters->FindParam(&scalar_list,NULL,NULL,InversionGradientScalingEnum);
+        femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);   _assert_(num_controls);
+        numvertices=femmodel->vertices->NumberOfVertices();
+
+        /*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");
+        for(long i=0;i<nsize;i++){
+                if(X[i]>XU[i]) X[i]=XU[i];
+                if(X[i]<XL[i]) X[i]=XL[i];
+        }
+
+        /*Update control input*/
+        SetControlInputsFromVectorx(femmodel,X);
+
+        /*Compute Adjoint*/
+        AdjointCorePointerFromSolutionEnum(&adjointcore,solution_type);
+        adjointcore(femmodel);
+
+        /*Compute gradient*/
+        Gradjx(&G,&norm_list,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+
+        /*Compute scaling factor*/
+        IssmDouble scalar = scalar_list[0]/norm_list[0];
+        for(int i=1;i<num_controls;i++) scalar=min(scalar,scalar_list[i]/norm_list[i]);
+
+        /*Constrain Gradient*/
+        for(int i=0;i<num_controls;i++){
+                for(int j=0;j<numvertices;j++){
+                        G[i*numvertices+j] = scalar*G[i*numvertices+j];
+                }
+        }
+
+        /*Needed for output results (FIXME: should be placed 6 lines below)*/
+        ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,G);
+
+        for(long i=0;i<nsize;i++){
+                if(X[i]>=XU[i]) G[i]=0.;
+                if(X[i]<=XL[i]) G[i]=0.;
+        }
+
+        /*solve forward: (FIXME: not needed actually...)*/
+        switch(solution_type){
+                case SteadystateSolutionEnum:
+                        femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
+                        stressbalance_core(femmodel);   //We need a 3D velocity!! (vz is required for the next thermal run)
+                        break;
+                case StressbalanceSolutionEnum:
+                        femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
+                        solutionsequence_nonlinear(femmodel,conserve_loads); 
+                        break;
+                case BalancethicknessSolutionEnum:
+                        femmodel->SetCurrentConfiguration(BalancethicknessAnalysisEnum);
+                        solutionsequence_linear(femmodel); 
+                        break;
+                case BalancethicknessSoftSolutionEnum:
+                        /*NOTHING*/
+                        break;
+                case Balancethickness2SolutionEnum:
+                        femmodel->SetCurrentConfiguration(Balancethickness2AnalysisEnum);
+                        solutionsequence_linear(femmodel); 
+                        break;
+                default:
+                        _error_("Solution " << EnumToStringx(solution_type) << " not implemented yet");
+        }
+
+        /*Compute misfit for this velocity field.*/
+        femmodel->CostFunctionx(&J,&Jlist,NULL);
+        //_printf0_("f(x) = "<<setw(12)<<setprecision(7)<<*pf<<"  |  ");
+        //for(int i=0;i<num_responses;i++) _printf0_(" "<<setw(12)<<setprecision(7)<<Jlist[i]);
+        //_printf0_("\n");
+
+        /*Clean-up and return*/
+        xDelete<IssmDouble>(XU);
+        xDelete<IssmDouble>(XL);
+        xDelete<IssmDouble>(norm_list);
+        xDelete<IssmDouble>(scalar_list);
+        xDelete<IssmDouble>(Jlist);
+        *pG = G;
+        return J;
+}

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

@@ -106 +106 @@
 
         /*Get solution*/
-        SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X);
+        SetControlInputsFromVectorx(femmodel,X);
         ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,G);
         femmodel->OutputControlsx(&femmodel->results);
@@ -131 +131 @@
         FemModel   *femmodel  = (FemModel*)dzs;
 
-        /*Recover responses*/
-        int         num_responses;
-        int        *responses = NULL;
-        femmodel->parameters->FindParam(&responses,&num_responses,InversionCostFunctionsEnum);
+        /*Recover number of cost functions responses*/
+        int num_responses;
+        femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
 
         /*Constrain input vector*/
@@ -147 +146 @@
 
         /*Update control input*/
-        SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X);
+        SetControlInputsFromVectorx(femmodel,X);
 
         /*Recover some parameters*/
@@ -170 +169 @@
         if(indic==0){
                 /*dry run, no gradient required*/
-                xDelete<int>(responses);
                 xDelete<IssmDouble>(XU);
                 xDelete<IssmDouble>(XL);
@@ -193 +191 @@
 
         /*Clean-up and return*/
-        xDelete<int>(responses);
         xDelete<IssmDouble>(XU);
         xDelete<IssmDouble>(XL);

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

@@ -93 +93 @@
         G=new Vector<IssmDouble>(0); VecFree(&G->pvector->vector);
         TaoGetGradientVector(tao,&G->pvector->vector);
-        SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X);
+        SetControlInputsFromVectorx(femmodel,X);
         ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,G);
         femmodel->OutputControlsx(&femmodel->results);
@@ -113 +113 @@
         TaoFinalize();
 }
-int FormFunctionGradient(TaoSolver tao, Vec Xpetsc, IssmDouble *fcn,Vec G,void *userCtx){
+int FormFunctionGradient(TaoSolver tao, Vec Xpetsc, IssmDouble *fcn,Vec G,void *uservoid){
 
         /*Retreive arguments*/
         int                  solution_type;
-        AppCtx              *user            = (AppCtx *)userCtx;
+        AppCtx              *user            = (AppCtx *)uservoid;
         FemModel            *femmodel        = user->femmodel;
         Vector<IssmDouble>  *gradient        = NULL;
@@ -127 +127 @@
         /*Set new variable*/
         //VecView(X,PETSC_VIEWER_STDOUT_WORLD);
-        SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X);
+        SetControlInputsFromVectorx(femmodel,X);
         delete X;
 

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

@@ -73 +73 @@
 
                 /*Calculate j(k+alpha delta k) */
-                SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X);
+                SetControlInputsFromVectorx(femmodel,X);
                 solutioncore(femmodel);
                 femmodel->CostFunctionx(&j,NULL,NULL);

issm/trunk-jpl/src/c/modules/Gradjx/Gradjx.cpp

@@ -50 +50 @@
                 delete gradient_list[i];
         }
+        //gradient->Echo();
+        //_error_("S");
 
         /*Check that gradient is clean*/

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

@@ -16 +16 @@
         int         num_cm_responses;
         int        *control_type     = NULL;
+        int        *maxiter          = NULL;
         int        *cm_responses     = NULL;
         IssmDouble *cm_jump          = NULL;
         IssmDouble *optscal          = NULL;
-        IssmDouble *maxiter          = NULL;
 
         /*retrieve some parameters: */
@@ -56 +56 @@
                                 parameters->AddObject(new DoubleMatParam(InversionGradientScalingEnum,optscal,nsteps,num_control_type));
                                 parameters->AddObject(new DoubleVecParam(InversionStepThresholdEnum,cm_jump,nsteps));
-                                parameters->AddObject(new DoubleVecParam(InversionMaxiterPerStepEnum,maxiter,nsteps));
+                                parameters->AddObject(new IntVecParam(InversionMaxiterPerStepEnum,maxiter,nsteps));
                                 break;
                         case 1:/*TAO*/
@@ -83 +83 @@
                 xDelete<int>(control_type);
                 xDelete<int>(cm_responses);
+                xDelete<int>(maxiter);
                 iomodel->DeleteData(cm_jump,InversionStepThresholdEnum);
                 iomodel->DeleteData(optscal,InversionGradientScalingEnum);
-                iomodel->DeleteData(maxiter,InversionMaxiterPerStepEnum);
         }
 }

issm/trunk-jpl/src/c/modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.cpp

@@ -7 +7 @@
 #include "../../toolkits/toolkits.h"
 
-void SetControlInputsFromVectorx(Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,IssmDouble* vector){
+void SetControlInputsFromVectorx(FemModel* femmodel,IssmDouble* vector){
 
         int  num_controls;
@@ -13 +13 @@
 
         /*Retrieve some parameters*/
-        parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
-        parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
+        femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
+        femmodel->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
 
         for(int i=0;i<num_controls;i++){
-                for(int j=0;j<elements->Size();j++){
-                        Element* element=(Element*)elements->GetObjectByOffset(j);
+                for(int j=0;j<femmodel->elements->Size();j++){
+                        Element* element=(Element*)femmodel->elements->GetObjectByOffset(j);
                         element->SetControlInputsFromVector(vector,control_type[i],i);
                 }
@@ -26 +26 @@
 }
 
-void SetControlInputsFromVectorx(Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,Vector<IssmDouble>* vector){
+void SetControlInputsFromVectorx(FemModel* femmodel,Vector<IssmDouble>* vector){
 
-        IssmDouble* serial_vector=NULL;
-
-        serial_vector=vector->ToMPISerial();
-
-        SetControlInputsFromVectorx(elements,nodes, vertices, loads, materials, parameters,serial_vector);
-
-        /*Free ressources:*/
+        IssmDouble* serial_vector=vector->ToMPISerial();
+        SetControlInputsFromVectorx(femmodel,serial_vector);
         xDelete<IssmDouble>(serial_vector);
 }

issm/trunk-jpl/src/c/modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.h

@@ -8 +8 @@
 
 /* local prototypes: */
-void SetControlInputsFromVectorx(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads, Materials* materials,  Parameters* parameters,Vector<IssmDouble>* vector);
-void SetControlInputsFromVectorx(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads, Materials* materials,  Parameters* parameters,IssmDouble* vector);
+void SetControlInputsFromVectorx(FemModel* femmodel,Vector<IssmDouble>* vector);
+void SetControlInputsFromVectorx(FemModel* femmodel,IssmDouble* vector);
 
 #endif 

issm/trunk-jpl/src/c/shared/Numerics/BrentSearch.cpp

@@ -20 +20 @@
 #include "./isnan.h"
 
-void BrentSearch(IssmDouble* psearch_scalar,IssmDouble* pJ,OptPars* optpars,IssmDouble (*f)(IssmDouble,void*),void* optargs){
+void BrentSearch(IssmDouble** pJ,OptPars optpars,IssmDouble* X0,IssmDouble (*f)(IssmDouble*,void*),IssmDouble (*g)(IssmDouble**,IssmDouble*,void*),void* usr){
 
         /* This routine is optimizing a given function using Brent's method
@@ -26 +26 @@
 
         /*Intermediary*/
+        int        iter;
         IssmDouble si,gold,intervalgold,oldintervalgold;
-        IssmDouble parab_num,parab_den;
-        IssmDouble distance,cm_jump;
+        IssmDouble parab_num,parab_den,distance;
         IssmDouble fxmax,fxmin,fxbest;
         IssmDouble fx,fx1,fx2;
-        IssmDouble xmax,xmin,xbest;
-        IssmDouble x,x1,x2,xm;
+        IssmDouble x,x1,x2,xm,xbest;
         IssmDouble tol1,tol2,seps;
         IssmDouble tolerance = 1.e-4;
-        int        maxiter ,iter;
-        bool       loop= true,goldenflag;
 
         /*Recover parameters:*/
-        xmin=optpars->xmin;
-        xmax=optpars->xmax;
-        maxiter=optpars->maxiter;
-        cm_jump=optpars->cm_jump;
-
-        /*initialize counter and get response at the boundaries*/
-        iter=0;
-        fxmin = (*f)(xmin,optargs);
-        if (xIsNan<IssmDouble>(fxmin)) _error_("Function evaluation returned NaN");
-        cout<<setprecision(5);
-        if(VerboseControl()) _printf0_("\n");
-        if(VerboseControl()) _printf0_("       Iteration         x           f(x)       Tolerance         Procedure\n");
-        if(VerboseControl()) _printf0_("\n");
-        if(VerboseControl()) _printf0_("           N/A    "<<setw(12)<<xmin<<"  "<<setw(12)<<fxmin<<"           N/A         boundary\n");
-        fxmax = (*f)(xmax,optargs);
-        if (xIsNan<IssmDouble>(fxmax)) _error_("Function evaluation returned NaN");
-        if(VerboseControl()) _printf0_("           N/A    "<<setw(12)<<xmax<<"  "<<setw(12)<<fxmax<<"           N/A         boundary\n");
-
-        /*test if jump option activated and xmin==0*/
-        if (!xIsNan<IssmDouble>(cm_jump) && (xmin==0) && (fxmax/fxmin)<cm_jump){
-                *psearch_scalar=xmax;
-                *pJ=fxmax;
-                return;
+        int         nsteps  = optpars.nsteps;
+        int         nsize   = optpars.nsize;
+        IssmDouble  xmin    = optpars.xmin;
+        IssmDouble  xmax    = optpars.xmax;
+        int        *maxiter = optpars.maxiter;
+        IssmDouble *cm_jump = optpars.cm_jump;
+
+        /*Initialize gradient and controls*/
+        IssmDouble* G = NULL;
+        IssmDouble* J = xNew<IssmDouble>(nsteps);
+        IssmDouble* X = xNew<IssmDouble>(nsize);
+
+        /*start iterations*/
+        for(int n=0;n<nsteps;n++){
+
+                /*Reset some variables*/
+                iter = 0;
+                xmin = 0.;
+                xmax = 1.;
+                bool loop = true;
+                cout<<setprecision(5);
+
+                /*Get current Gradient at xmin=0*/
+                if(VerboseControl()) _printf0_("\n" << "   step " << n+1 << "/" << nsteps << "\n");
+                fxmin = (*g)(&G,X0,usr); if(xIsNan<IssmDouble>(fxmin)) _error_("Function evaluation returned NaN");
+                if(VerboseControl()) _printf0_("\n");
+                if(VerboseControl()) _printf0_("       Iteration         x           f(x)       Tolerance         Procedure\n");
+                if(VerboseControl()) _printf0_("\n");
+                if(VerboseControl()) _printf0_("           N/A    "<<setw(12)<<xmin<<"  "<<setw(12)<<fxmin<<"           N/A         boundary\n");
+                
+                /*Get f(xmax)*/
+                for(int i=0;i<nsize;i++) X[i]=X0[i]+xmax*G[i];
+                fxmax = (*f)(X,usr); if (xIsNan<IssmDouble>(fxmax)) _error_("Function evaluation returned NaN");
+                if(VerboseControl()) _printf0_("           N/A    "<<setw(12)<<xmax<<"  "<<setw(12)<<fxmax<<"           N/A         boundary\n");
+
+                /*test if jump option activated and xmin==0*/
+                if(!xIsNan<IssmDouble>(cm_jump[n]) && (xmin==0) && (fxmax/fxmin)<cm_jump[n]){
+                        for(int i=0;i<nsize;i++) X0[i]=X0[i]+xmax*G[i];
+                        xDelete<IssmDouble>(G);
+                        J[n]=fxmax;
+                        continue;
+                }
+
+                /*initialize optimization variables*/
+                seps=sqrt(DBL_EPSILON);    //precision of a IssmDouble
+                distance=0.0;              //new_x=old_x + distance
+                gold=0.5*(3.0-sqrt(5.0));  //gold = 1 - golden ratio
+                intervalgold=0.0;          //distance used by Golden procedure
+
+                /*1: initialize the values of the 4 x needed (x1,x2,x,xbest)*/
+                x1=xmin+gold*(xmax-xmin);
+                x2=x1;
+                xbest=x1;
+                x=xbest;
+
+                /*2: call the function to be evaluated*/
+                for(int i=0;i<nsize;i++) X[i]=X0[i]+x*G[i];
+                fxbest = (*f)(X,usr); if(xIsNan<IssmDouble>(fxbest)) _error_("Function evaluation returned NaN");
+                iter++;
+
+                /*3: update the other variables*/
+                fx1=fxbest;
+                fx2=fxbest;
+                /*xm is always in the middle of a and b*/
+                xm=0.5*(xmin+xmax);                           
+                /*update tolerances*/
+                tol1=seps*sqrt(pow(xbest,2))+tolerance/3.0;
+                tol2=2.0*tol1;
+
+                /*4: print result*/
+                if(VerboseControl())
+                 _printf0_("         "<<setw(5)<<iter<<"    "<<setw(12)<<xbest<<"  "<<setw(12)<<fxbest<<"  "<<setw(12)<<pow(pow(xbest-xm,2),0.5)<<"         initial\n");
+                if (!xIsNan<IssmDouble>(cm_jump[n]) && (xmin==0) && ((fxbest/fxmin)<cm_jump[n])){
+                        if(VerboseControl()) _printf0_("      optimization terminated: current x satisfies criteria 'cm_jump'=" << cm_jump[n] << "\n");
+                        loop=false;
+                }
+
+                while(loop){
+
+                        bool goldenflag=true;
+
+                        // Is a parabolic fit possible ?
+                        if (sqrt(pow(intervalgold,2))>tol1){
+
+                                // Yes, so fit parabola
+                                goldenflag=false;
+                                parab_num=(xbest-x1)*(xbest-x1)*(fxbest-fx2)-(xbest-x2)*(xbest-x2)*(fxbest-fx1);;
+                                parab_den=2.0*(xbest-x1)*(fxbest-fx2)-2.0*(xbest-x2)*(fxbest-fx1);
+
+                                //reverse p if necessary
+                                if(parab_den>0.0){ 
+                                        parab_num=-parab_num;
+                                }
+                                parab_den=sqrt(pow(parab_den,2));
+                                oldintervalgold=intervalgold;
+                                intervalgold=distance;
+
+                                // Is the parabola acceptable (we use seps here because in some configuration parab_num==parab_den*(xmax-xbest)
+                                // and the result is not repeatable anymore
+                                if (( sqrt(pow(parab_num,2)) < sqrt(pow(0.5*parab_den*oldintervalgold,2))) &&
+                                                        (parab_num>parab_den*(xmin-xbest)+seps) && 
+                                                        (parab_num<parab_den*(xmax-xbest)-seps)){
+
+                                        // Yes, parabolic interpolation step
+                                        distance=parab_num/parab_den;
+                                        x=xbest+distance;
+
+                                        // f must not be evaluated too close to min_x or max_x
+                                        if (((x-xmin)<tol2) || ((xmax-x)<tol2)){
+                                                if ((xm-xbest)<0.0) si=-1;
+                                                else                si=1;
+                                                //compute new distance
+                                                distance=tol1*si;
+                                        }
+                                }
+                                else{
+                                        // Not acceptable, must do a golden section step
+                                        goldenflag=true;
+                                }
+                        }
+
+                        //Golden procedure
+                        if(goldenflag){
+                                // compute the new distance d
+                                if(xbest>=xm){
+                                        intervalgold=xmin-xbest;    
+                                }
+                                else{ 
+                                        intervalgold=xmax-xbest;  
+                                }
+                                distance=gold*intervalgold;
+                        }
+
+                        // The function must not be evaluated too close to xbest
+                        if(distance<0) si=-1;
+                        else           si=1;
+                        if(sqrt(pow(distance,2))>tol1) x=xbest+si*sqrt(pow(distance,2));
+                        else                           x=xbest+si*tol1;
+
+                        //evaluate function on x
+                        for(int i=0;i<nsize;i++) X[i]=X0[i]+x*G[i];
+                        fx = (*f)(X,usr); if(xIsNan<IssmDouble>(fx)) _error_("Function evaluation returned NaN");
+                        iter++;
+
+                        // Update a, b, xm, x1, x2, tol1, tol2
+                        if (fx<=fxbest){
+                                if (x>=xbest) xmin=xbest;
+                                else          xmax=xbest;
+                                x1=x2;    fx1=fx2;
+                                x2=xbest; fx2=fxbest;
+                                xbest=x;  fxbest=fx;
+                        }
+                        else{ // fx > fxbest
+                                if (x<xbest) xmin=x;
+                                else         xmax=x;
+                                if ((fx<=fx2) || (x2==xbest)){
+                                        x1=x2; fx1=fx2;
+                                        x2=x;  fx2=fx;
+                                }
+                                else if ( (fx <= fx1) || (x1 == xbest) || (x1 == x2) ){
+                                        x1=x;  fx1=fx;
+                                }
+                        }
+                        xm = 0.5*(xmin+xmax);
+                        tol1=seps*pow(pow(xbest,2),0.5)+tolerance/3.0;
+                        tol2=2.0*tol1;
+                        if(VerboseControl())
+                         _printf0_("         "<<setw(5)<<iter<<"    "<<setw(12)<<x<<"  "<<setw(12)<<fx<<"  "<<setw(12)<<pow(pow(xbest-xm,2),0.5)<<
+                                                 "         "<<(goldenflag?"golden":"parabolic")<<"\n");
+
+                        /*Stop the optimization?*/
+                        if (sqrt(pow(xbest-xm,2)) < (tol2-0.5*(xmax-xmin))){
+                                if(VerboseControl()) _printf0_("      optimization terminated: current x satisfies criteria 'tolx'=" << tolerance << "\n");
+                                loop=false;
+                        }
+                        else if (iter>=maxiter[n]){
+                                if(VerboseControl()) _printf0_("      exiting: Maximum number of iterations has been exceeded  ('maxiter'=" << maxiter[n] << ")\n");
+                                loop=false;
+                        }
+                        else if (!xIsNan<IssmDouble>(cm_jump[n]) && (xmin==0) && ((fxbest/fxmin)<cm_jump[n])){
+                                if(VerboseControl()) _printf0_("      optimization terminated: current x satisfies criteria 'cm_jump'=" << cm_jump[n] << "\n");
+                                loop=false;
+                        }
+                        else{
+                                //continue
+                                loop=true;
+                        }
+                }//end while
+
+                //Now, check that the value on the boundaries are not better than current fxbest
+                if (fxbest>fxmin){
+                        xbest=optpars.xmin; fxbest=fxmin;
+                }
+                if (fxbest>fxmax){
+                        xbest=optpars.xmax; fxbest=fxmax;
+                }
+
+                /*Assign output pointers: */
+                for(int i=0;i<nsize;i++) X0[i]=X0[i]+xbest*G[i];
+                xDelete<IssmDouble>(G);
+                J[n]=fxbest;
         }
-
-        /*initialize optimization variables*/
-        seps=sqrt(DBL_EPSILON);    //precision of a IssmDouble
-        distance=0.0;              //new_x=old_x + distance
-        gold=0.5*(3.0-sqrt(5.0));  //gold = 1 - golden ratio
-        intervalgold=0.0;          //distance used by Golden procedure
-
-        /*1: initialize the values of the 4 x needed (x1,x2,x,xbest)*/
-        x1=xmin+gold*(xmax-xmin);
-        x2=x1;
-        xbest=x1;
-        x=xbest;
-
-        /*2: call the function to be evaluated*/
-        fxbest = (*f)(x,optargs);
-        if(xIsNan<IssmDouble>(fxbest)) _error_("Function evaluation returned NaN");
-        iter=iter+1;
-
-        /*3: update the other variables*/
-        fx1=fxbest;
-        fx2=fxbest;
-        /*xm is always in the middle of a and b*/
-        xm=0.5*(xmin+xmax);                           
-        /*update tolerances*/
-        tol1=seps*sqrt(pow(xbest,2))+tolerance/3.0;
-        tol2=2.0*tol1;
-
-        /*4: print result*/
-        if(VerboseControl())
-         _printf0_("         "<<setw(5)<<iter<<"    "<<setw(12)<<xbest<<"  "<<setw(12)<<fxbest<<"  "<<setw(12)<<pow(pow(xbest-xm,2),0.5)<<"         initial\n");
-        if (!xIsNan<IssmDouble>(cm_jump) && (xmin==0) && ((fxbest/fxmin)<cm_jump)){
-                if(VerboseControl()) _printf0_("      optimization terminated: current x satisfies criteria 'cm_jump'=" << cm_jump << "\n");
-                loop=false;
-        }
-
-        while(loop){
-
-                goldenflag=true;
-
-                // Is a parabolic fit possible ?
-                if (sqrt(pow(intervalgold,2))>tol1){
-
-                        // Yes, so fit parabola
-                        goldenflag=false;
-                        parab_num=(xbest-x1)*(xbest-x1)*(fxbest-fx2)-(xbest-x2)*(xbest-x2)*(fxbest-fx1);;
-                        parab_den=2.0*(xbest-x1)*(fxbest-fx2)-2.0*(xbest-x2)*(fxbest-fx1);
-
-                        //reverse p if necessary
-                        if(parab_den>0.0){ 
-                                parab_num=-parab_num;
-                        }
-                        parab_den=sqrt(pow(parab_den,2));
-                        oldintervalgold=intervalgold;
-                        intervalgold=distance;
-
-                        // Is the parabola acceptable (we use seps here because in some configuration parab_num==parab_den*(xmax-xbest)
-                        // and the result is not repeatable anymore
-                        if (( sqrt(pow(parab_num,2)) < sqrt(pow(0.5*parab_den*oldintervalgold,2))) &&
-                                                (parab_num>parab_den*(xmin-xbest)+seps) && 
-                                                (parab_num<parab_den*(xmax-xbest)-seps)){
-
-                                // Yes, parabolic interpolation step
-                                distance=parab_num/parab_den;
-                                x=xbest+distance;
-
-                                // f must not be evaluated too close to min_x or max_x
-                                if (((x-xmin)<tol2) || ((xmax-x)<tol2)){
-                                        if ((xm-xbest)<0.0) si=-1;
-                                        else                si=1;
-                                        //compute new distance
-                                        distance=tol1*si;
-                                }
-                        }
-                        else{
-                                // Not acceptable, must do a golden section step
-                                goldenflag=true;
-                        }
-                }
-
-                //Golden procedure
-                if(goldenflag){
-                        // compute the new distance d
-                        if(xbest>=xm){
-                                intervalgold=xmin-xbest;    
-                        }
-                        else{ 
-                                intervalgold=xmax-xbest;  
-                        }
-                        distance=gold*intervalgold;
-                }
-
-                // The function must not be evaluated too close to xbest
-                if(distance<0) si=-1;
-                else           si=1;
-                if(sqrt(pow(distance,2))>tol1) x=xbest+si*sqrt(pow(distance,2));
-                else                           x=xbest+si*tol1;
-
-                //evaluate function on x
-                fx = (*f)(x,optargs);
-                if(xIsNan<IssmDouble>(fx)) _error_("Function evaluation returned NaN");
-                iter=iter+1;
-
-                // Update a, b, xm, x1, x2, tol1, tol2
-                if (fx<=fxbest){
-                        if (x>=xbest) xmin=xbest;
-                        else          xmax=xbest;
-                        x1=x2;    fx1=fx2;
-                        x2=xbest; fx2=fxbest;
-                        xbest=x;  fxbest=fx;
-                }
-                else{ // fx > fxbest
-                        if (x<xbest) xmin=x;
-                        else         xmax=x;
-                        if ((fx<=fx2) || (x2==xbest)){
-                                x1=x2; fx1=fx2;
-                                x2=x;  fx2=fx;
-                        }
-                        else if ( (fx <= fx1) || (x1 == xbest) || (x1 == x2) ){
-                                x1=x;  fx1=fx;
-                        }
-                }
-                xm = 0.5*(xmin+xmax);
-                tol1=seps*pow(pow(xbest,2),0.5)+tolerance/3.0;
-                tol2=2.0*tol1;
-                if(VerboseControl())
-                 _printf0_("         "<<setw(5)<<iter<<"    "<<setw(12)<<x<<"  "<<setw(12)<<fx<<"  "<<setw(12)<<pow(pow(xbest-xm,2),0.5)<<
-                                         "         "<<(goldenflag?"golden":"parabolic")<<"\n");
-
-                /*Stop the optimization?*/
-                if (sqrt(pow(xbest-xm,2)) < (tol2-0.5*(xmax-xmin))){
-                        if(VerboseControl()) _printf0_("      optimization terminated: current x satisfies criteria 'tolx'=" << tolerance << "\n");
-                        loop=false;
-                }
-                else if (iter>=maxiter){
-                        if(VerboseControl()) _printf0_("      exiting: Maximum number of iterations has been exceeded  ('maxiter'=" << maxiter << ")\n");
-                        loop=false;
-                }
-                else if (!xIsNan<IssmDouble>(cm_jump) && (xmin==0) && ((fxbest/fxmin)<cm_jump)){
-                        if(VerboseControl()) _printf0_("      optimization terminated: current x satisfies criteria 'cm_jump'=" << cm_jump << "\n");
-                        loop=false;
-                }
-                else{
-                        //continue
-                        loop=true;
-                }
-        }//end while
-
-        //Now, check that the value on the boundaries are not better than current fxbest
-        if (fxbest>fxmin){
-                xbest=optpars->xmin; fxbest=fxmin;
-        }
-        if (fxbest>fxmax){
-                xbest=optpars->xmax; fxbest=fxmax;
-        }
-
-        /*Assign output pointers: */
-        *psearch_scalar=xbest;
-        *pJ=fxbest;
+        
+        /*return*/
+        xDelete<IssmDouble>(X);
+        *pJ=J;
 }

issm/trunk-jpl/src/c/shared/Numerics/OptPars.h

@@ -10 +10 @@
 struct OptPars{
 
-        IssmDouble xmin;
-        IssmDouble xmax;
-        IssmDouble cm_jump;
-        int maxiter;
+        IssmDouble  xmin;
+        IssmDouble  xmax;
+        IssmDouble *cm_jump;
+        int* maxiter;
+        int  nsteps;
+        int  nsize;
 
 };

issm/trunk-jpl/src/c/shared/Numerics/numerics.h

@@ -30 +30 @@
 int         min(int a,int b);
 int         max(int a,int b);
-void        BrentSearch(IssmDouble *psearch_scalar,IssmDouble*pJ,OptPars*optpars,IssmDouble (*f)(IssmDouble,void*),void* optargs);
+void        BrentSearch(IssmDouble** pJ,OptPars optpars,IssmDouble* X0,IssmDouble (*f)(IssmDouble*,void*),IssmDouble (*g)(IssmDouble**,IssmDouble*,void*),void* usr);
 void        cross(IssmDouble *result,IssmDouble*vector1,IssmDouble*vector2);
 void        XZvectorsToCoordinateSystem(IssmDouble *T,IssmDouble*xzvectors);

issm/trunk-jpl/src/m/classes/inversion.m

@@ -25 +25 @@
         end
         methods
-         function createxml(obj,fid) % {{{
-            fprintf(fid, '<!-- inversion -->\n');            
-                    
-            % inversion parameters
-            fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="inversion parameters">','<section name="inversion" />');                    
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="iscontrol" type="',class(obj.iscontrol),'" default="',convert2str(obj.iscontrol),'">','     <section name="inversion" />','     <help> is inversion activated? </help>','  </parameter>');
-            
-            % incompleteadjoing drop-down (0 or 1)
-            fprintf(fid,'%s\n%s\n%s\n%s\n','  <parameter key ="incomplete_adjoint" type="alternative" optional="false">','     <section name="inversion" />','     <help> 1: linear viscosity, 0: non-linear viscosity </help>');
-            fprintf(fid,'%s\n','       <option value="0" type="string" default="true"> </option>');
-            fprintf(fid,'%s\n%s\n','       <option value="1" type="string" default="false"> </option>','</parameter>');
-            
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="control_parameters" type="',class(obj.control_parameters),'" default="',convert2str(obj.control_parameters),'">','     <section name="inversion" />','     <help> ex: {''FrictionCoefficient''}, or {''MaterialsRheologyBbar''} </help>','  </parameter>');
-                
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="nsteps" type="',class(obj.nsteps),'" default="',convert2str(obj.nsteps),'">','     <section name="inversion" />','     <help> number of optimization searches </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="cost_functions" type="',class(obj.cost_functions),'" default="',convert2str(obj.cost_functions),'">','     <section name="inversion" />','     <help> indicate the type of response for each optimization step  </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="cost_functions_coefficients" type="',class(obj.cost_functions_coefficients),'" default="',convert2str(obj.cost_functions_coefficients),'">','     <section name="inversion" />','     <help> cost_functions_coefficients applied to the misfit of each vertex and for each control_parameter </help>','  </parameter>');
-                
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="cost_function_threshold" type="',class(obj.cost_function_threshold),'" default="',convert2str(obj.cost_function_threshold),'">','     <section name="inversion" />','     <help> misfit convergence criterion. Default is 1%, NaN if not applied </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="maxiter_per_step" type="',class(obj.maxiter_per_step),'" default="',convert2str(obj.maxiter_per_step),'">','     <section name="inversion" />','     <help> maximum iterations during each optimization step  </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="gradient_scaling" type="',class(obj.gradient_scaling),'" default="',convert2str(obj.gradient_scaling),'">','     <section name="inversion" />','     <help> scaling factor on gradient direction during optimization, for each optimization step </help>','  </parameter>');
-               
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="step_threshold" type="',class(obj.step_threshold),'" default="',convert2str(obj.step_threshold),'">','     <section name="inversion" />','     <help> decrease threshold for misfit, default is 30% </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="min_parameters" type="',class(obj.min_parameters),'" default="',convert2str(obj.min_parameters),'">','     <section name="inversion" />','     <help> absolute minimum acceptable value of the inversed parameter on each vertex </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="max_parameters" type="',class(obj.max_parameters),'" default="',convert2str(obj.max_parameters),'">','     <section name="inversion" />','     <help> absolute maximum acceptable value of the inversed parameter on each vertex </help>','  </parameter>');
-               
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="vx_obs" type="',class(obj.vx_obs),'" default="',convert2str(obj.vx_obs),'">','     <section name="inversion" />','     <help> observed velocity x component [m/yr] </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="vy_obs" type="',class(obj.vy_obs),'" default="',convert2str(obj.vy_obs),'">','     <section name="inversion" />','     <help> observed velocity y component [m/yr]  </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="vel_obs" type="',class(obj.vel_obs),'" default="',convert2str(obj.vel_obs),'">','     <section name="inversion" />','     <help> observed velocity magnitude [m/yr] </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="thickness_obs" type="',class(obj.thickness_obs),'" default="',convert2str(obj.thickness_obs),'">','     <section name="inversion" />','     <help> observed thickness [m]) </help>','  </parameter>');
-               
-            fprintf(fid,'%s\n%s\n','</frame>');    
-            
-            fprintf(fid,'%s\n%s\n%s\n','<frame key="2" label="Available cost functions">','<section name="inversion" />');                    
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="SurfaceAbsVelMisfit" type="','string','" default="','101','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="SurfaceRelVelMisfit" type="','string','" default="','102','">','     <section name="inversion" />','     <help>   </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="SurfaceLogVelMisfit" type="','string','" default="','103','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
-                
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="SurfaceLogVxVyMisfit" type="','string','" default="','104','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="SurfaceAverageVelMisfit" type="','string','" default="','105','">','     <section name="inversion" />','     <help>   </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="ThicknessAbsMisfit" type="','string','" default="','106','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
-               
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="DragCoefficientAbsGradient" type="','string','" default="','107','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="RheologyBbarAbsGradient" type="','string','" default="','108','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="ThicknessAbsGradient" type="','string','" default="','109','">','     <section name="inversion" />','     <help> </help>','  </parameter>');
-               
-            fprintf(fid,'%s\n%s\n','</frame>');    
-        
-        end % }}}       
+                function createxml(obj,fid) % {{{
+                        fprintf(fid, '<!-- inversion -->\n');            
+
+                        % inversion parameters
+                        fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="inversion parameters">','<section name="inversion" />');                    
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="iscontrol" type="',class(obj.iscontrol),'" default="',convert2str(obj.iscontrol),'">','     <section name="inversion" />','     <help> is inversion activated? </help>','  </parameter>');
+
+                        % incompleteadjoing drop-down (0 or 1)
+                        fprintf(fid,'%s\n%s\n%s\n%s\n','  <parameter key ="incomplete_adjoint" type="alternative" optional="false">','     <section name="inversion" />','     <help> 1: linear viscosity, 0: non-linear viscosity </help>');
+                        fprintf(fid,'%s\n','       <option value="0" type="string" default="true"> </option>');
+                        fprintf(fid,'%s\n%s\n','       <option value="1" type="string" default="false"> </option>','</parameter>');
+
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="control_parameters" type="',class(obj.control_parameters),'" default="',convert2str(obj.control_parameters),'">','     <section name="inversion" />','     <help> ex: {''FrictionCoefficient''}, or {''MaterialsRheologyBbar''} </help>','  </parameter>');
+
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="nsteps" type="',class(obj.nsteps),'" default="',convert2str(obj.nsteps),'">','     <section name="inversion" />','     <help> number of optimization searches </help>','  </parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="cost_functions" type="',class(obj.cost_functions),'" default="',convert2str(obj.cost_functions),'">','     <section name="inversion" />','     <help> indicate the type of response for each optimization step  </help>','  </parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="cost_functions_coefficients" type="',class(obj.cost_functions_coefficients),'" default="',convert2str(obj.cost_functions_coefficients),'">','     <section name="inversion" />','     <help> cost_functions_coefficients applied to the misfit of each vertex and for each control_parameter </help>','  </parameter>');
+
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="cost_function_threshold" type="',class(obj.cost_function_threshold),'" default="',convert2str(obj.cost_function_threshold),'">','     <section name="inversion" />','     <help> misfit convergence criterion. Default is 1%, NaN if not applied </help>','  </parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="maxiter_per_step" type="',class(obj.maxiter_per_step),'" default="',convert2str(obj.maxiter_per_step),'">','     <section name="inversion" />','     <help> maximum iterations during each optimization step  </help>','  </parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="gradient_scaling" type="',class(obj.gradient_scaling),'" default="',convert2str(obj.gradient_scaling),'">','     <section name="inversion" />','     <help> scaling factor on gradient direction during optimization, for each optimization step </help>','  </parameter>');
+
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="step_threshold" type="',class(obj.step_threshold),'" default="',convert2str(obj.step_threshold),'">','     <section name="inversion" />','     <help> decrease threshold for misfit, default is 30% </help>','  </parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="min_parameters" type="',class(obj.min_parameters),'" default="',convert2str(obj.min_parameters),'">','     <section name="inversion" />','     <help> absolute minimum acceptable value of the inversed parameter on each vertex </help>','  </parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="max_parameters" type="',class(obj.max_parameters),'" default="',convert2str(obj.max_parameters),'">','     <section name="inversion" />','     <help> absolute maximum acceptable value of the inversed parameter on each vertex </help>','  </parameter>');
+
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="vx_obs" type="',class(obj.vx_obs),'" default="',convert2str(obj.vx_obs),'">','     <section name="inversion" />','     <help> observed velocity x component [m/yr] </help>','  </parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="vy_obs" type="',class(obj.vy_obs),'" default="',convert2str(obj.vy_obs),'">','     <section name="inversion" />','     <help> observed velocity y component [m/yr]  </help>','  </parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="vel_obs" type="',class(obj.vel_obs),'" default="',convert2str(obj.vel_obs),'">','     <section name="inversion" />','     <help> observed velocity magnitude [m/yr] </help>','  </parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="thickness_obs" type="',class(obj.thickness_obs),'" default="',convert2str(obj.thickness_obs),'">','     <section name="inversion" />','     <help> observed thickness [m]) </help>','  </parameter>');
+
+                        fprintf(fid,'%s\n%s\n','</frame>');    
+
+                        fprintf(fid,'%s\n%s\n%s\n','<frame key="2" label="Available cost functions">','<section name="inversion" />');                    
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="SurfaceAbsVelMisfit" type="','string','" default="','101','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="SurfaceRelVelMisfit" type="','string','" default="','102','">','     <section name="inversion" />','     <help>   </help>','  </parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="SurfaceLogVelMisfit" type="','string','" default="','103','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
+
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="SurfaceLogVxVyMisfit" type="','string','" default="','104','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="SurfaceAverageVelMisfit" type="','string','" default="','105','">','     <section name="inversion" />','     <help>   </help>','  </parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="ThicknessAbsMisfit" type="','string','" default="','106','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
+
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="DragCoefficientAbsGradient" type="','string','" default="','107','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="RheologyBbarAbsGradient" type="','string','" default="','108','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="ThicknessAbsGradient" type="','string','" default="','109','">','     <section name="inversion" />','     <help> </help>','  </parameter>');
+
+                        fprintf(fid,'%s\n%s\n','</frame>');    
+
+                end % }}}       
                 function obj = inversion(varargin) % {{{
                         switch nargin
@@ -196 +196 @@
                         if ~obj.iscontrol, return; end
                         WriteData(fid,'object',obj,'fieldname','nsteps','format','Integer');
-                        WriteData(fid,'object',obj,'fieldname','maxiter_per_step','format','DoubleMat','mattype',3);
+                        WriteData(fid,'object',obj,'fieldname','maxiter_per_step','format','IntMat','mattype',3);
                         WriteData(fid,'object',obj,'fieldname','cost_functions_coefficients','format','DoubleMat','mattype',1);
                         WriteData(fid,'object',obj,'fieldname','gradient_scaling','format','DoubleMat','mattype',3);