Index: /issm/trunk-jpl/src/c/classes/Cfsurfacesquare.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Cfsurfacesquare.cpp (revision 22514) +++ /issm/trunk-jpl/src/c/classes/Cfsurfacesquare.cpp (revision 22515) @@ -119,6 +119,8 @@ IssmDouble J=0.; IssmDouble J_sum=0.; - + printf("-------------- file: Cfsurfacesquare.cpp line: %i\n",__LINE__); + _printf_("current time: "<elements->Size();i++){ Element* element=(Element*)femmodel->elements->GetObjectByOffset(i); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 22514) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 22515) @@ -1602,4 +1602,43 @@ } else _error_("Cannot add input for vector type " << vector_type << " (not supported)"); +} +/*}}}*/ +void Element::ControlInputCreate(IssmDouble* vector,IssmDouble* min_vector,IssmDouble* max_vector,IoModel* iomodel,int M,int N,int input_enum,int id){/*{{{*/ + + /*Intermediaries*/ + int numvertices = this->GetNumberOfVertices(); + int *vertexids = xNew(numvertices); + IssmDouble *values = xNew(numvertices); + IssmDouble *values_min = xNew(numvertices); + IssmDouble *values_max = xNew(numvertices); + + /*Some sanity checks*/ + _assert_(vector); + _assert_(min_vector); + _assert_(max_vector); + + /*For now we only support nodal vectors*/ + if(M!=iomodel->numberofvertices) _error_("not supported"); + if(N!=1) _error_("not supported"); + + /*Recover vertices ids needed to initialize inputs*/ + _assert_(iomodel->elements); + for(int i=0;i(iomodel->elements[numvertices*this->Sid()+i]); //ids for vertices are in the elements array from Matlab + } + + /*Are we in transient or static? */ + if(M==iomodel->numberofvertices){ + for(int i=0;iAddControlInput(input_enum,values,min_vector,max_vector,P1Enum,id); + } + + /*clean up*/ + xDelete(values); + xDelete(vertexids); } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 22514) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 22515) @@ -122,4 +122,5 @@ void InputChangeName(int enum_type,int enum_type_old); void InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code); + void ControlInputCreate(IssmDouble* doublearray,IssmDouble* independents_min,IssmDouble* independents_max,IoModel* iomodel,int M,int N,int input_enum,int id); void DatasetInputAdd(int enum_type,IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code,int input_enum); void InputDuplicate(int original_enum,int new_enum); @@ -188,4 +189,6 @@ virtual void AddBasalInput(int input_enum, IssmDouble* values, int interpolation_enum)=0; virtual void AddInput(int input_enum, IssmDouble* values, int interpolation_enum)=0; + virtual void AddControlInput(int input_enum, IssmDouble* values,IssmDouble* values_min,IssmDouble* values_max, int interpolation_enum,int id){_error_("not supported yet");}; + virtual void DatasetInputCreate(IssmDouble* array,int M,int N,int* individual_enums,int num_inputs,IoModel* iomodel,int input_enum){_error_("not supported");}; virtual void AverageOntoPartition(Vector* partition_contributions,Vector* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part)=0; virtual void BasalNodeIndices(int* pnumindices,int** pindices,int finiteelement){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 22514) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 22515) @@ -1400,5 +1400,5 @@ for(j=0;jData("md.inversion.min_parameters")[(penta_vertex_ids[j]-1)*num_control_type+i]/yts; for(j=0;jData("md.inversion.max_parameters")[(penta_vertex_ids[j]-1)*num_control_type+i]/yts; - this->inputs->AddInput(new ControlInput(BalancethicknessThickeningRateEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + this->inputs->AddInput(new ControlInput(BalancethicknessThickeningRateEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,P1Enum,i+1)); } break; @@ -1408,5 +1408,5 @@ for(j=0;jData("md.inversion.min_parameters")[(penta_vertex_ids[j]-1)*num_control_type+i]/yts; for(j=0;jData("md.inversion.max_parameters")[(penta_vertex_ids[j]-1)*num_control_type+i]/yts; - this->inputs->AddInput(new ControlInput(VxEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + this->inputs->AddInput(new ControlInput(VxEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,P1Enum,i+1)); } break; @@ -1416,5 +1416,5 @@ for(j=0;jData("md.inversion.min_parameters")[(penta_vertex_ids[j]-1)*num_control_type+i]/yts; for(j=0;jData("md.inversion.max_parameters")[(penta_vertex_ids[j]-1)*num_control_type+i]/yts; - this->inputs->AddInput(new ControlInput(VyEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + this->inputs->AddInput(new ControlInput(VyEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,P1Enum,i+1)); } break; @@ -1424,5 +1424,5 @@ for(j=0;jData("md.inversion.min_parameters")[(penta_vertex_ids[j]-1)*num_control_type+i]; for(j=0;jData("md.inversion.max_parameters")[(penta_vertex_ids[j]-1)*num_control_type+i]; - this->inputs->AddInput(new ControlInput(FrictionCoefficientEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + this->inputs->AddInput(new ControlInput(FrictionCoefficientEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,P1Enum,i+1)); } break; @@ -1432,5 +1432,5 @@ for(j=0;jData("md.inversion.min_parameters")[(penta_vertex_ids[j]-1)*num_control_type+i]; for(j=0;jData("md.inversion.max_parameters")[(penta_vertex_ids[j]-1)*num_control_type+i]; - this->inputs->AddInput(new ControlInput(MaterialsRheologyBEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + this->inputs->AddInput(new ControlInput(MaterialsRheologyBEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,P1Enum,i+1)); } break; @@ -1440,5 +1440,5 @@ for(j=0;jData("md.inversion.min_parameters")[(penta_vertex_ids[j]-1)*num_control_type+i]; for(j=0;jData("md.inversion.max_parameters")[(penta_vertex_ids[j]-1)*num_control_type+i]; - this->inputs->AddInput(new ControlInput(DamageDEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + this->inputs->AddInput(new ControlInput(DamageDEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,P1Enum,i+1)); } break; Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.cpp (revision 22514) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.cpp (revision 22515) @@ -394,5 +394,5 @@ for(j=0;jData("md.inversion.min_parameters")[(tetra_vertex_ids[j]-1)*num_control_type+i]/yts; for(j=0;jData("md.inversion.max_parameters")[(tetra_vertex_ids[j]-1)*num_control_type+i]/yts; - this->inputs->AddInput(new ControlInput(BalancethicknessThickeningRateEnum,TetraInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + this->inputs->AddInput(new ControlInput(BalancethicknessThickeningRateEnum,TetraInputEnum,nodeinputs,cmmininputs,cmmaxinputs,P1Enum,i+1)); } break; @@ -402,5 +402,5 @@ for(j=0;jData("md.inversion.min_parameters")[(tetra_vertex_ids[j]-1)*num_control_type+i]/yts; for(j=0;jData("md.inversion.max_parameters")[(tetra_vertex_ids[j]-1)*num_control_type+i]/yts; - this->inputs->AddInput(new ControlInput(VxEnum,TetraInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + this->inputs->AddInput(new ControlInput(VxEnum,TetraInputEnum,nodeinputs,cmmininputs,cmmaxinputs,P1Enum,i+1)); } break; @@ -410,5 +410,5 @@ for(j=0;jData("md.inversion.min_parameters")[(tetra_vertex_ids[j]-1)*num_control_type+i]/yts; for(j=0;jData("md.inversion.max_parameters")[(tetra_vertex_ids[j]-1)*num_control_type+i]/yts; - this->inputs->AddInput(new ControlInput(VyEnum,TetraInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + this->inputs->AddInput(new ControlInput(VyEnum,TetraInputEnum,nodeinputs,cmmininputs,cmmaxinputs,P1Enum,i+1)); } break; @@ -418,5 +418,5 @@ for(j=0;jData("md.inversion.min_parameters")[(tetra_vertex_ids[j]-1)*num_control_type+i]; for(j=0;jData("md.inversion.max_parameters")[(tetra_vertex_ids[j]-1)*num_control_type+i]; - this->inputs->AddInput(new ControlInput(FrictionCoefficientEnum,TetraInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + this->inputs->AddInput(new ControlInput(FrictionCoefficientEnum,TetraInputEnum,nodeinputs,cmmininputs,cmmaxinputs,P1Enum,i+1)); } break; @@ -426,5 +426,5 @@ for(j=0;jData("md.inversion.min_parameters")[(tetra_vertex_ids[j]-1)*num_control_type+i]; for(j=0;jData("md.inversion.max_parameters")[(tetra_vertex_ids[j]-1)*num_control_type+i]; - this->inputs->AddInput(new ControlInput(MaterialsRheologyBEnum,TetraInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + this->inputs->AddInput(new ControlInput(MaterialsRheologyBEnum,TetraInputEnum,nodeinputs,cmmininputs,cmmaxinputs,P1Enum,i+1)); } break; @@ -434,5 +434,5 @@ for(j=0;jData("md.inversion.min_parameters")[(tetra_vertex_ids[j]-1)*num_control_type+i]; for(j=0;jData("md.inversion.max_parameters")[(tetra_vertex_ids[j]-1)*num_control_type+i]; - this->inputs->AddInput(new ControlInput(DamageDEnum,TetraInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + this->inputs->AddInput(new ControlInput(DamageDEnum,TetraInputEnum,nodeinputs,cmmininputs,cmmaxinputs,P1Enum,i+1)); } break; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 22514) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 22515) @@ -168,4 +168,34 @@ _assert_(this->inputs); this->inputs->AddInput(new TriaInput(input_enum,values,interpolation_enum)); +} +/*}}}*/ +void Tria::AddControlInput(int input_enum,IssmDouble* values,IssmDouble* values_min,IssmDouble* values_max, int interpolation_enum,int id){/*{{{*/ + + /*Call inputs method*/ + _assert_(this->inputs); + this->inputs->AddInput(new ControlInput(input_enum,TriaInputEnum,values,values_min,values_max,interpolation_enum,id)); +} +/*}}}*/ +void Tria::DatasetInputCreate(IssmDouble* array,int M,int N,int* individual_enums,int num_inputs,IoModel* iomodel,int input_enum){/*{{{*/ + + IssmDouble nodeinputs[NUMVERTICES]; + if(num_inputs<1) _error_("Cannot create a DatasetInput of size <1"); + if(M!=iomodel->numberofvertices) _error_("not supported yet"); + if(N!=num_inputs) _error_("sizes are not consistent"); + + int tria_vertex_ids[3]; + + for(int k=0;k<3;k++){ + tria_vertex_ids[k]=reCast(iomodel->elements[3*this->Sid()+k]); //ids for vertices are in the elements array from Matlab + } + /*Create inputs and add to DataSetInput*/ + DatasetInput* datasetinput=new DatasetInput(input_enum); + for(int i=0;iAddInput(new TriaInput(input_enum,nodeinputs,P1Enum),individual_enums[i]); + } + + /*Add datasetinput to element inputs*/ + this->inputs->AddInput(datasetinput); } /*}}}*/ @@ -1911,5 +1941,5 @@ IssmDouble cmmininputs[3]; IssmDouble cmmaxinputs[3]; - bool control_analysis = false; + bool control_analysis,ad_analysis = false; int num_control_type,num_responses; char** controls = NULL; @@ -1919,6 +1949,11 @@ iomodel->FindConstant(&yts,"md.constants.yts"); iomodel->FindConstant(&control_analysis,"md.inversion.iscontrol"); - if(control_analysis) iomodel->FindConstant(&num_control_type,"md.inversion.num_control_parameters"); - if(control_analysis) iomodel->FindConstant(&num_responses,"md.inversion.num_cost_functions"); + iomodel->FindConstant(&ad_analysis, "md.autodiff.isautodiff"); + if(control_analysis && !ad_analysis) iomodel->FindConstant(&num_control_type,"md.inversion.num_control_parameters"); + if(control_analysis && !ad_analysis) iomodel->FindConstant(&num_responses,"md.inversion.num_cost_functions"); + if(control_analysis && ad_analysis) iomodel->FindConstant(&num_control_type,"md.autodiff.num_independent_objects"); + if(control_analysis && ad_analysis) iomodel->FindConstant(&num_responses,"md.autodiff.num_dependent_objects"); + + /*Recover vertices ids needed to initialize inputs*/ @@ -1934,5 +1969,7 @@ /*Control Inputs*/ if (control_analysis){ - iomodel->FindConstant(&controls,NULL,"md.inversion.control_parameters"); + if(!ad_analysis)iomodel->FindConstant(&controls,NULL,"md.inversion.control_parameters"); + if(ad_analysis)iomodel->FindConstant(&controls,NULL,"md.autodiff.independent_object_names"); + for(i=0;iData("md.inversion.min_parameters")[(tria_vertex_ids[j]-1)*num_control_type+i]/yts; for(j=0;jData("md.inversion.max_parameters")[(tria_vertex_ids[j]-1)*num_control_type+i]/yts; - this->inputs->AddInput(new ControlInput(BalancethicknessThickeningRateEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + this->inputs->AddInput(new ControlInput(BalancethicknessThickeningRateEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,P1Enum,i+1)); } break; @@ -1952,5 +1989,5 @@ for(j=0;jData("md.inversion.min_parameters")[(tria_vertex_ids[j]-1)*num_control_type+i]/yts; for(j=0;jData("md.inversion.max_parameters")[(tria_vertex_ids[j]-1)*num_control_type+i]/yts; - this->inputs->AddInput(new ControlInput(VxEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + this->inputs->AddInput(new ControlInput(VxEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,P1Enum,i+1)); } break; @@ -1960,5 +1997,5 @@ for(j=0;jData("md.inversion.min_parameters")[(tria_vertex_ids[j]-1)*num_control_type+i]/yts; for(j=0;jData("md.inversion.max_parameters")[(tria_vertex_ids[j]-1)*num_control_type+i]/yts; - this->inputs->AddInput(new ControlInput(VyEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + this->inputs->AddInput(new ControlInput(VyEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,P1Enum,i+1)); } break; @@ -1968,5 +2005,5 @@ for(j=0;jData("md.inversion.min_parameters")[(tria_vertex_ids[j]-1)*num_control_type+i]; for(j=0;jData("md.inversion.max_parameters")[(tria_vertex_ids[j]-1)*num_control_type+i]; - this->inputs->AddInput(new ControlInput(ThicknessEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + this->inputs->AddInput(new ControlInput(ThicknessEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,P1Enum,i+1)); } break; @@ -1976,5 +2013,5 @@ for(j=0;jData("md.inversion.min_parameters")[(tria_vertex_ids[j]-1)*num_control_type+i]; for(j=0;jData("md.inversion.max_parameters")[(tria_vertex_ids[j]-1)*num_control_type+i]; - this->inputs->AddInput(new ControlInput(BalancethicknessSpcthicknessEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + this->inputs->AddInput(new ControlInput(BalancethicknessSpcthicknessEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,P1Enum,i+1)); } break; @@ -1984,5 +2021,5 @@ for(j=0;jData("md.inversion.min_parameters")[(tria_vertex_ids[j]-1)*num_control_type+i]; for(j=0;jData("md.inversion.max_parameters")[(tria_vertex_ids[j]-1)*num_control_type+i]; - this->inputs->AddInput(new ControlInput(BalancethicknessOmegaEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + this->inputs->AddInput(new ControlInput(BalancethicknessOmegaEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,P1Enum,i+1)); } break; @@ -1992,5 +2029,5 @@ for(j=0;jData("md.inversion.min_parameters")[(tria_vertex_ids[j]-1)*num_control_type+i]; for(j=0;jData("md.inversion.max_parameters")[(tria_vertex_ids[j]-1)*num_control_type+i]; - this->inputs->AddInput(new ControlInput(FrictionCoefficientEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + this->inputs->AddInput(new ControlInput(FrictionCoefficientEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,P1Enum,i+1)); } break; @@ -2000,5 +2037,5 @@ for(j=0;jData("md.inversion.min_parameters")[(tria_vertex_ids[j]-1)*num_control_type+i]; for(j=0;jData("md.inversion.max_parameters")[(tria_vertex_ids[j]-1)*num_control_type+i]; - this->inputs->AddInput(new ControlInput(MaterialsRheologyBbarEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + this->inputs->AddInput(new ControlInput(MaterialsRheologyBbarEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,P1Enum,i+1)); } break; @@ -2008,5 +2045,5 @@ for(j=0;jData("md.inversion.min_parameters")[(tria_vertex_ids[j]-1)*num_control_type+i]; for(j=0;jData("md.inversion.max_parameters")[(tria_vertex_ids[j]-1)*num_control_type+i]; - this->inputs->AddInput(new ControlInput(DamageDbarEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + this->inputs->AddInput(new ControlInput(DamageDbarEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,P1Enum,i+1)); } break; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 22514) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 22515) @@ -162,4 +162,6 @@ void AddBasalInput(int input_enum, IssmDouble* values, int interpolation_enum); void AddInput(int input_enum, IssmDouble* values, int interpolation_enum); + void AddControlInput(int input_enum, IssmDouble* values,IssmDouble* values_min,IssmDouble* values_max, int interpolation_enum,int id); + void DatasetInputCreate(IssmDouble* array,int M,int N,int* individual_enums,int num_inputs,IoModel* iomodel,int input_enum); IssmDouble GetArea(void); IssmDouble GetArea3D(void); Index: /issm/trunk-jpl/src/c/classes/Inputs/ControlInput.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/ControlInput.cpp (revision 22514) +++ /issm/trunk-jpl/src/c/classes/Inputs/ControlInput.cpp (revision 22515) @@ -22,24 +22,26 @@ } /*}}}*/ -ControlInput::ControlInput(int in_enum_type,int enum_input,IssmDouble* pvalues,IssmDouble* pmin,IssmDouble* pmax,int id){/*{{{*/ +ControlInput::ControlInput(int in_enum_type,int input_layout_enum,IssmDouble* pvalues,IssmDouble* pmin,IssmDouble* pmax,int interp,int id){/*{{{*/ control_id=id; enum_type=in_enum_type; - switch(enum_input){ + _assert_(interp==P1Enum); + + switch(input_layout_enum){ case TriaInputEnum: - values =new TriaInput(enum_type,pvalues,P1Enum); - savedvalues=new TriaInput(enum_type,pvalues,P1Enum); - minvalues =new TriaInput(enum_type,pmin,P1Enum); - maxvalues =new TriaInput(enum_type,pmax,P1Enum); + values =new TriaInput(enum_type,pvalues,interp); + savedvalues=new TriaInput(enum_type,pvalues,interp); + minvalues =new TriaInput(enum_type,pmin,interp); + maxvalues =new TriaInput(enum_type,pmax,interp); break; case PentaInputEnum: - values =new PentaInput(enum_type,pvalues,P1Enum); - savedvalues=new PentaInput(enum_type,pvalues,P1Enum); - minvalues =new PentaInput(enum_type,pmin,P1Enum); - maxvalues =new PentaInput(enum_type,pmax,P1Enum); + values =new PentaInput(enum_type,pvalues,interp); + savedvalues=new PentaInput(enum_type,pvalues,interp); + minvalues =new PentaInput(enum_type,pmin,interp); + maxvalues =new PentaInput(enum_type,pmax,interp); break; default: - _error_("Input of Enum " << EnumToStringx(enum_input) << " not supported yet by ControlInput"); + _error_("Input of Enum " << EnumToStringx(input_layout_enum) << " not supported yet by ControlInput"); } gradient =NULL; Index: /issm/trunk-jpl/src/c/classes/Inputs/ControlInput.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/ControlInput.h (revision 22514) +++ /issm/trunk-jpl/src/c/classes/Inputs/ControlInput.h (revision 22515) @@ -26,5 +26,5 @@ /*ControlInput constructors, destructors: {{{*/ ControlInput(); - ControlInput(int enum_type,int enum_input,IssmDouble* pvalues,IssmDouble* pmin,IssmDouble* pmax,int id); + ControlInput(int enum_type,int enum_input,IssmDouble* pvalues,IssmDouble* pmin,IssmDouble* pmax,int interp,int id); ~ControlInput(); /*}}}*/ Index: /issm/trunk-jpl/src/c/cores/controladm1qn3_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/controladm1qn3_core.cpp (revision 22514) +++ /issm/trunk-jpl/src/c/cores/controladm1qn3_core.cpp (revision 22515) @@ -355,7 +355,9 @@ long niter = long(maxsteps); /*Maximum number of iterations*/ long nsim = long(maxiter);/*Maximum number of function calls*/ - + /*Get initial guess*/ Vector *Xpetsc = NULL; + + /*THIS IS WHERE IT FAILS*/ GetPassiveVectorFromControlInputsx(&Xpetsc,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value"); X = Xpetsc->ToMPISerial(); @@ -363,5 +365,5 @@ delete Xpetsc; _assert_(intn==numberofvertices*num_controls); - + /*Get problem dimension and initialize gradient and initial guess*/ long n = long(intn); @@ -395,4 +397,5 @@ /*Initialize Gradient and cost function of M1QN3*/ indic = 4; /*gradient required*/ +printf("-------------- file: controladm1qn3_core.cpp line: %i\n",__LINE__); simul_ad(&indic,&n,X,&f,G,izs,rzs,(void*)&mystruct); /*Estimation of the expected decrease in f during the first iteration*/ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp (revision 22514) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp (revision 22515) @@ -28,8 +28,4 @@ if(control_analysis){ - /*How many controls and how many responses?*/ - parameters->AddObject(iomodel->CopyConstantObject("md.inversion.num_control_parameters",InversionNumControlParametersEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.inversion.num_cost_functions",InversionNumCostFunctionsEnum)); - /*What solution type?*/ if(solution_type==SteadystateSolutionEnum){ @@ -40,27 +36,90 @@ } - /*recover controls and convert to Enums*/ - iomodel->FindConstant(&controls,&num_controls,"md.inversion.control_parameters"); - if(num_controls<1) _error_("no controls found"); - int* control_enums=xNew(num_controls); - for(int i=0;i(controls[i]); + switch(inversiontype){ + _printf_("Inversiontype: "<AddObject(iomodel->CopyConstantObject("md.inversion.num_control_parameters",InversionNumControlParametersEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.inversion.num_cost_functions",InversionNumCostFunctionsEnum)); + + /*recover controls and convert to Enums*/ + iomodel->FindConstant(&controls,&num_controls,"md.inversion.control_parameters"); + if(num_controls<1) _error_("no controls found"); + int* control_enums=xNew(num_controls); + for(int i=0;i(controls[i]); + } + xDelete(controls); + parameters->AddObject(new IntVecParam(InversionControlParametersEnum,control_enums,num_controls)); + + iomodel->FindConstant(&cm_responses,&num_costfunc,"md.inversion.cost_functions"); + if(num_costfunc<1) _error_ ("no cost functions found"); + int* costfunc_enums=xNew(num_costfunc); + for(int i=0;i(cm_responses[i]); + } + xDelete(cm_responses); + parameters->AddObject(new IntVecParam(InversionCostFunctionsEnum,costfunc_enums,num_costfunc)); + + xDelete(control_enums); + xDelete(costfunc_enums); + + break; } - xDelete(controls); - parameters->AddObject(new IntVecParam(InversionControlParametersEnum,control_enums,num_controls)); + case 4:/*AD M1QN3*/ + { + /*Intermediaries*/ + int num_independent_objects,M; + char** names = NULL; + + /*this is done somewhere else*/ + parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.num_independent_objects",InversionNumControlParametersEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.num_dependent_objects",InversionNumCostFunctionsEnum)); + + /*Step1: create controls (independents)*/ + iomodel->FetchData(&num_independent_objects,"md.autodiff.num_independent_objects"); + _assert_(num_independent_objects>0); + iomodel->FetchData(&names,&M,"md.autodiff.independent_object_names"); + _assert_(M==num_independent_objects); + int* ind_enums=xNew(num_independent_objects); + for(int i=0;i(names[i]); + } - iomodel->FindConstant(&cm_responses,&num_costfunc,"md.inversion.cost_functions"); - if(num_costfunc<1) _error_ ("no cost functions found"); - int* costfunc_enums=xNew(num_costfunc); - for(int i=0;i(cm_responses[i]); + parameters->AddObject(new IntVecParam(InversionControlParametersEnum,ind_enums,num_independent_objects)); + //iomodel->FetchData(&num_costfunc,"md.numberedcostfunction.num_cost_functions"); + _assert_(num_costfunc>0); + iomodel->FindConstant(&cm_responses,&num_costfunc,"md.autodiff.dependent_object_names"); + if(num_costfunc<1) _error_ ("no cost functions found"); + int* costfunc_enums=xNew(num_costfunc); + for(int i=0;i(cm_responses[i]); + } + xDelete(cm_responses); + parameters->AddObject(new IntVecParam(InversionCostFunctionsEnum,costfunc_enums,num_costfunc)); + + iomodel->FetchData(&control_scaling_factors,NULL,NULL,"md.autodiff.independent_scaling_factors"); + parameters->AddObject(new DoubleVecParam(InversionControlScalingFactorsEnum,control_scaling_factors,num_independent_objects)); + + /*cleanup*/ + for(int i=0;i(names[i]); + } + xDelete(names); + xDelete(ind_enums); + xDelete(costfunc_enums); + break; + } + default: + _error_("not supported"); } - xDelete(cm_responses); - parameters->AddObject(new IntVecParam(InversionCostFunctionsEnum,costfunc_enums,num_costfunc)); - - xDelete(control_enums); - xDelete(costfunc_enums); /*Inversion type specifics*/ @@ -107,12 +166,9 @@ parameters->AddObject(iomodel->CopyConstantObject("md.inversion.maxsteps",InversionMaxstepsEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.inversion.maxiter",InversionMaxiterEnum)); - iomodel->FetchData(&control_scaling_factors,NULL,NULL,"md.inversion.control_scaling_factors"); - parameters->AddObject(new DoubleVecParam(InversionControlScalingFactorsEnum,control_scaling_factors,num_controls)); - break; + break; default: _error_("not supported"); } - xDelete(control_enums); xDelete(maxiter); xDelete(control_scaling_factors); Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp (revision 22514) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp (revision 22515) @@ -9,6 +9,8 @@ #include "../ModelProcessorx.h" + +#if !defined(_HAVE_ADOLC_) void UpdateElementsAndMaterialsControl(Elements* elements,Materials* materials, IoModel* iomodel){ - + /*Intermediary*/ bool control_analysis; @@ -21,4 +23,6 @@ char **cost_functions = NULL; + if(_HAVE_ADOLC_) return; + /*Fetch parameters: */ iomodel->FindConstant(&control_analysis,"md.inversion.iscontrol"); @@ -27,5 +31,5 @@ /*Now, return if no control*/ if(!control_analysis) return; - + /*Process controls and convert from string to enums*/ iomodel->FindConstant(&controls,&num_controls,"md.inversion.control_parameters"); @@ -125,2 +129,72 @@ xDelete(controls); } +#else +void UpdateElementsAndMaterialsControl(Elements* elements,Materials* materials, IoModel* iomodel){ + + /*Intermediaries*/ + int num_independent_objects,M,N; + char** names = NULL; + int* types = NULL; + IssmDouble* independent = NULL; + IssmDouble** independents_min = NULL; + IssmDouble** independents_max = NULL; + + /*Step1: create controls (independents)*/ + iomodel->FetchData(&num_independent_objects,"md.autodiff.num_independent_objects"); + _assert_(num_independent_objects>0); + iomodel->FetchData(&names,&M,"md.autodiff.independent_object_names"); + _assert_(M==num_independent_objects); + iomodel->FetchData(&types,NULL,NULL,"md.autodiff.independent_object_types"); + + + /*TODO: fetch min and max*/ + independents_min = xNew(num_independent_objects); + independents_max = xNew(num_independent_objects); + for(int i=0;iFetchData(&independent,&M,&N,iofieldname); + iomodel->FetchData(&independents_min[i],&M,&N,"md.autodiff.independent_min_parameters"); + iomodel->FetchData(&independents_max[i],&M,&N,"md.autodiff.independent_max_parameters"); + + _assert_(independent); + + for(int j=0;jSize();j++){ + Element* element=xDynamicCast(elements->GetObjectByOffset(j)); + element->ControlInputCreate(independent,independents_min[i],independents_max[i],iomodel,M,N,input_enum,i+1); + } + xDelete(independent); + } + else{ + _error_("Independent cannot be of size " << types[i]); + } + } + + /*cleanup*/ + for(int i=0;i(names[i]); + xDelete(independents_min[i]); + xDelete(independents_max[i]); + } + xDelete(independents_min); + xDelete(independents_max); + xDelete(names); + xDelete(types); + + /*Step2: create cost functions (dependents)*/ + + return; +} +#endif Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (revision 22514) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (revision 22515) @@ -15,4 +15,5 @@ bool control_analysis; bool dakota_analysis; + bool adolc_analysis; int nnat,dummy; int* nature=NULL; @@ -22,4 +23,5 @@ iomodel->FindConstant(&dakota_analysis,"md.qmu.isdakota"); iomodel->FindConstant(&materials_type,"md.materials.type"); + iomodel->FindConstant(&adolc_analysis,"md.autodiff.isautodiff"); /*Did we already create the elements? : */ @@ -30,5 +32,6 @@ /*Create elements*/ - if(control_analysis)iomodel->FetchData(2,"md.inversion.min_parameters","md.inversion.max_parameters"); + if(control_analysis && !adolc_analysis)iomodel->FetchData(2,"md.inversion.min_parameters","md.inversion.max_parameters"); + switch(iomodel->meshelementtype){ case TriaEnum: @@ -218,6 +221,5 @@ /*Free data: */ - iomodel->DeleteData(7,"md.mesh.upperelements","md.mesh.lowerelements","md.material.rheology_B", - "md.material.rheology_n","md.damage.D","md.inversion.min_parameters","md.inversion.max_parameters"); + iomodel->DeleteData(7,"md.mesh.upperelements","md.mesh.lowerelements","md.material.rheology_B","md.material.rheology_n","md.damage.D","md.inversion.min_parameters","md.inversion.max_parameters"); /*Add new constant material property to materials, at the end: */ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateOutputDefinitions.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateOutputDefinitions.cpp (revision 22514) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateOutputDefinitions.cpp (revision 22515) @@ -507,10 +507,13 @@ /*Intermediary*/ - int numout; - char **ncf_name_s = NULL; - char **ncf_definitionstring_s = NULL; - char **cost_functions = NULL; - int cost_function ,domaintype; - int num_cost_functions; + int numout,numout2; + char **ncf_name_s = NULL; + char **ncf_definitionstring_s = NULL; + char **cost_functions = NULL; + IssmDouble **cost_functions_weights = NULL; + int* cost_functions_weights_M = NULL; + int* cost_functions_weights_N = NULL; + int cost_function,domaintype; + int num_cost_functions; /*Process cost functions and convert from string to enums*/ @@ -523,9 +526,8 @@ } - //iomodel->FetchData(1,"md.numberedcostfunction.cost_functions_coefficients"); - iomodel->FetchMultipleData(&ncf_name_s,&numout,"md.numberedcostfunction.name"); - iomodel->FetchMultipleData(&ncf_definitionstring_s,&numout,"md.numberedcostfunction.definitionstring"); - if(numout>1) _error_("not implemented yet, check code here"); + iomodel->FetchMultipleData(&ncf_definitionstring_s,&numout2,"md.numberedcostfunction.definitionstring"); _assert_(numout2 == numout); + iomodel->FetchMultipleData(&cost_functions_weights,&cost_functions_weights_M,&cost_functions_weights_N,&numout2,"md.numberedcostfunction.cost_functions_coefficients"); _assert_(numout2 == numout); + if(numout!=1) _error_("not implemented yet, check code here"); /*Fetch Observations */ @@ -543,16 +545,17 @@ if(domaintype!=Domain2DverticalEnum) iomodel->FetchDataToInput(elements,"md.numberedcostfunction.vy_obs",InversionVyObsEnum); } - } for(j=0;jSize();k++){ + Element* element=xDynamicCast(elements->GetObjectByOffset(k)); + element->DatasetInputCreate(cost_functions_weights[i],cost_functions_weights_M[i],cost_functions_weights_N[i],cost_function_enums,num_cost_functions,iomodel,InversionCostFunctionsCoefficientsEnum); + } output_definitions->AddObject(new Numberedcostfunction(ncf_name_s[j],StringToEnumx(ncf_definitionstring_s[j]),num_cost_functions,cost_function_enums)); - - } - + } /*Free data: */ - //iomodel->DeleteData(1,"md.inversion.cost_functions_coefficients"); iomodel->DeleteData(2,"md.numberedcostfunction.name","md.numberedcostfunction.definitionstring"); xDelete(cost_function_enums); @@ -564,7 +567,11 @@ xDelete(ncf_name_s[j]); xDelete(ncf_definitionstring_s[j]); + xDelete(cost_functions_weights[j]); } xDelete(ncf_name_s); xDelete(ncf_definitionstring_s); + xDelete(cost_functions_weights_M); + xDelete(cost_functions_weights_N); + xDelete(cost_functions_weights); } /*}}}*/ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 22514) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 22515) @@ -347,4 +347,5 @@ #ifdef _HAVE_ADOLC_ if(VerboseMProcessor()) _printf0_(" starting autodiff parameters \n"); + printf("-------------- file: CreateParameters.cpp line: %i\n",__LINE__); CreateParametersAutodiff(parameters,iomodel); if(VerboseMProcessor()) _printf0_(" ending autodiff parameters \n"); Index: /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp (revision 22514) +++ /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp (revision 22515) @@ -11,124 +11,124 @@ /*output*/ char* fieldname = NULL; - int param_enum = -1; + int input_enum = -1; if(strcmp(string_in,"Thickness")==0){ const char* field = "md.geometry.thickness"; - param_enum = ThicknessEnum; + input_enum = ThicknessEnum; fieldname=xNew((strlen(field)+1)); xMemCpy(fieldname,field,(strlen(field)+1)); } else if(strcmp(string_in,"MaterialsRheologyB")==0){ const char* field = "md.materials.rheology_B"; - param_enum = MaterialsRheologyBEnum; + input_enum = MaterialsRheologyBEnum; fieldname=xNew((strlen(field)+1)); xMemCpy(fieldname,field,(strlen(field)+1)); } else if(strcmp(string_in,"SmbMassBalance")==0){ const char* field = "md.smb.mass_balance"; - param_enum = SmbMassBalanceEnum; + input_enum = SmbMassBalanceEnum; fieldname=xNew((strlen(field)+1)); xMemCpy(fieldname,field,(strlen(field)+1)); } else if(strcmp(string_in,"SmbAccumulation")==0){ const char* field = "md.smb.accumulation"; - param_enum = SmbAccumulationEnum; + input_enum = SmbAccumulationEnum; fieldname=xNew((strlen(field)+1)); xMemCpy(fieldname,field,(strlen(field)+1)); } else if(strcmp(string_in,"SmbMelt")==0){ const char* field = "md.smb.melt"; - param_enum = SmbMeltEnum; + input_enum = SmbMeltEnum; fieldname=xNew((strlen(field)+1)); xMemCpy(fieldname,field,(strlen(field)+1)); } else if(strcmp(string_in,"SmbRefreeze")==0){ const char* field = "md.smb.refreeze"; - param_enum = SmbRefreezeEnum; + input_enum = SmbRefreezeEnum; fieldname=xNew((strlen(field)+1)); xMemCpy(fieldname,field,(strlen(field)+1)); } else if(strcmp(string_in,"SmbRunoff")==0){ const char* field = "md.smb.runoff"; - param_enum = SmbRunoffEnum; + input_enum = SmbRunoffEnum; fieldname=xNew((strlen(field)+1)); xMemCpy(fieldname,field,(strlen(field)+1)); } else if(strcmp(string_in,"SmbEvaporation")==0){ const char* field = "md.smb.evaporation"; - param_enum = SmbEvaporationEnum; + input_enum = SmbEvaporationEnum; fieldname=xNew((strlen(field)+1)); xMemCpy(fieldname,field,(strlen(field)+1)); } else if(strcmp(string_in,"SmbTa")==0){ const char* field = "md.smb.Ta"; - param_enum = SmbTaEnum; + input_enum = SmbTaEnum; fieldname=xNew((strlen(field)+1)); xMemCpy(fieldname,field,(strlen(field)+1)); } else if(strcmp(string_in,"SmbV")==0){ const char* field = "md.smb.V"; - param_enum = SmbVEnum; + input_enum = SmbVEnum; fieldname=xNew((strlen(field)+1)); xMemCpy(fieldname,field,(strlen(field)+1)); } else if(strcmp(string_in,"SmbDswrf")==0){ const char* field = "md.smb.dswrf"; - param_enum = SmbDswrfEnum; + input_enum = SmbDswrfEnum; fieldname=xNew((strlen(field)+1)); xMemCpy(fieldname,field,(strlen(field)+1)); } else if(strcmp(string_in,"SmbDlwrf")==0){ const char* field = "md.smb.dlwrf"; - param_enum = SmbDlwrfEnum; + input_enum = SmbDlwrfEnum; fieldname=xNew((strlen(field)+1)); xMemCpy(fieldname,field,(strlen(field)+1)); } else if(strcmp(string_in,"SmbP")==0){ const char* field = "md.smb.P"; - param_enum = SmbPEnum; + input_enum = SmbPEnum; fieldname=xNew((strlen(field)+1)); xMemCpy(fieldname,field,(strlen(field)+1)); } else if(strcmp(string_in,"SmbEAir")==0){ const char* field = "md.smb.eAir"; - param_enum = SmbEAirEnum; + input_enum = SmbEAirEnum; fieldname=xNew((strlen(field)+1)); xMemCpy(fieldname,field,(strlen(field)+1)); } else if(strcmp(string_in,"SmbPAir")==0){ const char* field = "md.smb.pAir"; - param_enum = SmbPAirEnum; + input_enum = SmbPAirEnum; fieldname=xNew((strlen(field)+1)); xMemCpy(fieldname,field,(strlen(field)+1)); } else if(strcmp(string_in,"SmbVz")==0){ const char* field = "md.smb.Vz"; - param_enum = SmbVzEnum; + input_enum = SmbVzEnum; fieldname=xNew((strlen(field)+1)); xMemCpy(fieldname,field,(strlen(field)+1)); } else if(strcmp(string_in,"SmbTz")==0){ const char* field = "md.smb.Tz"; - param_enum = SmbTzEnum; + input_enum = SmbTzEnum; fieldname=xNew((strlen(field)+1)); xMemCpy(fieldname,field,(strlen(field)+1)); } else if(strcmp(string_in,"SmbC")==0){ const char* field = "md.smb.C"; - param_enum = SmbCEnum; + input_enum = SmbCEnum; fieldname=xNew((strlen(field)+1)); xMemCpy(fieldname,field,(strlen(field)+1)); } else if(strcmp(string_in,"BasalforcingsFloatingiceMeltingRate")==0){ const char* field = "md.basalforcings.floatingice_melting_rate"; - param_enum = BasalforcingsFloatingiceMeltingRateEnum; - fieldname=xNew((strlen(field)+1)); xMemCpy(fieldname,field,(strlen(field)+1)); - } - else if(strcmp(string_in,"FrictionCoefficient")==0){ + input_enum = BasalforcingsFloatingiceMeltingRateEnum; + fieldname=xNew((strlen(field)+1)); xMemCpy(fieldname,field,(strlen(field)+1)); + } + else if(strcmp(string_in,"FrictionCoefficient")==0 || strcmp(string_in,"md.friction.coefficient")==0){ const char* field = "md.friction.coefficient"; - param_enum = FrictionCoefficientEnum; + input_enum = FrictionCoefficientEnum; fieldname=xNew((strlen(field)+1)); xMemCpy(fieldname,field,(strlen(field)+1)); } else if(strcmp(string_in,"Vx")==0){ const char* field = "md.initialization.vx"; - param_enum = VxEnum; + input_enum = VxEnum; fieldname=xNew((strlen(field)+1)); xMemCpy(fieldname,field,(strlen(field)+1)); } else if(strcmp(string_in,"Vy")==0){ const char* field = "md.initialization.vy"; - param_enum = VyEnum; + input_enum = VyEnum; fieldname=xNew((strlen(field)+1)); xMemCpy(fieldname,field,(strlen(field)+1)); } else if(strcmp(string_in,"BalancethicknessThickeningRate")==0){ const char* field = "md.balancethickness.thickening_rate"; - param_enum = BalancethicknessThickeningRateEnum; + input_enum = BalancethicknessThickeningRateEnum; fieldname=xNew((strlen(field)+1)); xMemCpy(fieldname,field,(strlen(field)+1)); } else if(strcmp(string_in,"BalancethicknessSpcthickness")==0){ const char* field = "md.balancethickness.spcthickness"; - param_enum = BalancethicknessSpcthicknessEnum; + input_enum = BalancethicknessSpcthicknessEnum; fieldname=xNew((strlen(field)+1)); xMemCpy(fieldname,field,(strlen(field)+1)); } @@ -138,5 +138,5 @@ /*Assign output pointers*/ - *out_enum = param_enum; + *out_enum = input_enum; *pfield = fieldname; return; Index: /issm/trunk-jpl/src/m/classes/adm1qn3inversion.m =================================================================== --- /issm/trunk-jpl/src/m/classes/adm1qn3inversion.m (revision 22514) +++ /issm/trunk-jpl/src/m/classes/adm1qn3inversion.m (revision 22515) @@ -7,20 +7,8 @@ properties (SetAccess=public) iscontrol = 0 - control_parameters = NaN - control_scaling_factors = NaN maxsteps = 0 maxiter = 0 dxmin = 0 gttol = 0 - cost_functions = NaN - cost_functions_coefficients = NaN - min_parameters = NaN - max_parameters = NaN - vx_obs = NaN - vy_obs = NaN - vz_obs = NaN - vel_obs = NaN - thickness_obs = NaN - surface_obs = NaN end @@ -37,11 +25,4 @@ end % }}} function self = extrude(self,md) % {{{ - self.vx_obs=project3d(md,'vector',self.vx_obs,'type','node'); - self.vy_obs=project3d(md,'vector',self.vy_obs,'type','node'); - self.vel_obs=project3d(md,'vector',self.vel_obs,'type','node'); - self.thickness_obs=project3d(md,'vector',self.thickness_obs,'type','node'); - if numel(self.cost_functions_coefficients)>1,self.cost_functions_coefficients=project3d(md,'vector',self.cost_functions_coefficients,'type','node');end; - if numel(self.min_parameters)>1,self.min_parameters=project3d(md,'vector',self.min_parameters,'type','node');end; - if numel(self.max_parameters)>1,self.max_parameters=project3d(md,'vector',self.max_parameters,'type','node');end; end % }}} function self = setdefaultparameters(self) % {{{ @@ -50,15 +31,9 @@ %parameter to be inferred by control methods (only %drag and B are supported yet) - self.control_parameters={'FrictionCoefficient'}; - - %Scaling factor for each control - self.control_scaling_factors=1; + %self.control_parameters={'FrictionCoefficient'}; %number of iterations self.maxsteps=20; self.maxiter=40; - - %several responses can be used: - self.cost_functions=101; %m1qn3 parameters @@ -75,59 +50,19 @@ md = checkmessage(md,['M1QN3 has not been installed, ISSM needs to be reconfigured and recompiled with M1QN3']); end - num_controls=numel(md.inversion.control_parameters); - num_costfunc=size(md.inversion.cost_functions,2); md = checkfield(md,'fieldname','inversion.iscontrol','values',[0 1]); - md = checkfield(md,'fieldname','inversion.control_parameters','cell',1,'values',supportedcontrols()); - md = checkfield(md,'fieldname','inversion.control_scaling_factors','size',[1 num_controls],'>',0,'NaN',1,'Inf',1); md = checkfield(md,'fieldname','inversion.maxsteps','numel',1,'>=',0); md = checkfield(md,'fieldname','inversion.maxiter','numel',1,'>=',0); md = checkfield(md,'fieldname','inversion.dxmin','numel',1,'>',0); md = checkfield(md,'fieldname','inversion.gttol','numel',1,'>',0); - md = checkfield(md,'fieldname','inversion.cost_functions','size',[1 num_costfunc],'values',supportedcostfunctions()); - md = checkfield(md,'fieldname','inversion.cost_functions_coefficients','size',[md.mesh.numberofvertices num_costfunc],'>=',0); - md = checkfield(md,'fieldname','inversion.min_parameters','size',[md.mesh.numberofvertices num_controls]); - md = checkfield(md,'fieldname','inversion.max_parameters','size',[md.mesh.numberofvertices num_controls]); - if strcmp(solution,'BalancethicknessSolution') - md = checkfield(md,'fieldname','inversion.thickness_obs','size',[md.mesh.numberofvertices 1],'NaN',1,'Inf',1); - md = checkfield(md,'fieldname','inversion.surface_obs','size',[md.mesh.numberofvertices 1],'NaN',1,'Inf',1); - elseif strcmp(solution,'BalancethicknessSoftSolution') - md = checkfield(md,'fieldname','inversion.thickness_obs','size',[md.mesh.numberofvertices 1],'NaN',1,'Inf',1); - else - md = checkfield(md,'fieldname','inversion.vx_obs','size',[md.mesh.numberofvertices 1],'NaN',1,'Inf',1); - if ~strcmp(domaintype(md.mesh),'2Dvertical'), - md = checkfield(md,'fieldname','inversion.vy_obs','size',[md.mesh.numberofvertices 1],'NaN',1,'Inf',1); - end - end - end % }}} + end % }}} function disp(self) % {{{ disp(sprintf(' adm1qn3inversion parameters:')); fielddisplay(self,'iscontrol','is inversion activated?'); - fielddisplay(self,'control_parameters','ex: {''FrictionCoefficient''}, or {''MaterialsRheologyBbar''}'); - fielddisplay(self,'control_scaling_factors','order of magnitude of each control (useful for multi-parameter optimization)'); fielddisplay(self,'maxsteps','maximum number of iterations (gradient computation)'); fielddisplay(self,'maxiter','maximum number of Function evaluation (forward run)'); fielddisplay(self,'dxmin','convergence criterion: two points less than dxmin from eachother (sup-norm) are considered identical'); fielddisplay(self,'gttol','convergence criterion: ||g(X)||/||g(X0)|| (g(X0): gradient at initial guess X0)'); - fielddisplay(self,'cost_functions','indicate the type of response for each optimization step'); - fielddisplay(self,'cost_functions_coefficients','cost_functions_coefficients applied to the misfit of each vertex and for each control_parameter'); - fielddisplay(self,'min_parameters','absolute minimum acceptable value of the inversed parameter on each vertex'); - fielddisplay(self,'max_parameters','absolute maximum acceptable value of the inversed parameter on each vertex'); - fielddisplay(self,'vx_obs','observed velocity x component [m/yr]'); - fielddisplay(self,'vy_obs','observed velocity y component [m/yr]'); - fielddisplay(self,'vel_obs','observed velocity magnitude [m/yr]'); - fielddisplay(self,'thickness_obs','observed thickness [m]'); - fielddisplay(self,'surface_obs','observed surface elevation [m]'); - disp('Available cost functions:'); - disp(' 101: SurfaceAbsVelMisfit'); - disp(' 102: SurfaceRelVelMisfit'); - disp(' 103: SurfaceLogVelMisfit'); - disp(' 104: SurfaceLogVxVyMisfit'); - disp(' 105: SurfaceAverageVelMisfit'); - disp(' 201: ThicknessAbsMisfit'); - disp(' 501: DragCoefficientAbsGradient'); - disp(' 502: RheologyBbarAbsGradient'); - disp(' 503: ThicknessAbsGradient'); end % }}} function marshall(self,prefix,md,fid) % {{{ @@ -138,54 +73,17 @@ WriteData(fid,prefix,'name','md.inversion.type','data',4,'format','Integer'); if ~self.iscontrol, return; end - WriteData(fid,prefix,'object',self,'class','inversion','fieldname','control_scaling_factors','format','DoubleMat','mattype',3); WriteData(fid,prefix,'object',self,'class','inversion','fieldname','maxsteps','format','Integer'); WriteData(fid,prefix,'object',self,'class','inversion','fieldname','maxiter','format','Integer'); WriteData(fid,prefix,'object',self,'class','inversion','fieldname','dxmin','format','Double'); WriteData(fid,prefix,'object',self,'class','inversion','fieldname','gttol','format','Double'); - WriteData(fid,prefix,'object',self,'class','inversion','fieldname','cost_functions_coefficients','format','DoubleMat','mattype',1); - WriteData(fid,prefix,'object',self,'class','inversion','fieldname','min_parameters','format','DoubleMat','mattype',3); - WriteData(fid,prefix,'object',self,'class','inversion','fieldname','max_parameters','format','DoubleMat','mattype',3); - WriteData(fid,prefix,'object',self,'class','inversion','fieldname','vx_obs','format','DoubleMat','mattype',1,'scale',1./yts); - WriteData(fid,prefix,'object',self,'class','inversion','fieldname','vy_obs','format','DoubleMat','mattype',1,'scale',1./yts); - WriteData(fid,prefix,'object',self,'class','inversion','fieldname','vz_obs','format','DoubleMat','mattype',1,'scale',1./yts); - WriteData(fid,prefix,'object',self,'class','inversion','fieldname','vel_obs','format','DoubleMat','mattype',1,'scale',1./yts); - if(numel(self.thickness_obs)==md.mesh.numberofelements), - mattype=2; - else - mattype=1; - end - WriteData(fid,prefix,'object',self,'class','inversion','fieldname','thickness_obs','format','DoubleMat','mattype',mattype); - WriteData(fid,prefix,'object',self,'class','inversion','fieldname','surface_obs','format','DoubleMat','mattype',mattype); - %process control parameters - num_control_parameters=numel(self.control_parameters); - WriteData(fid,prefix,'object',self,'fieldname','control_parameters','format','StringArray'); - WriteData(fid,prefix,'data',num_control_parameters,'name','md.inversion.num_control_parameters','format','Integer'); - - %process cost functions - num_cost_functions=size(self.cost_functions,2); - data=marshallcostfunctions(self.cost_functions); - WriteData(fid,prefix,'data',data,'name','md.inversion.cost_functions','format','StringArray'); - WriteData(fid,prefix,'data',num_cost_functions,'name','md.inversion.num_cost_functions','format','Integer'); end % }}} function savemodeljs(self,fid,modelname) % {{{ writejsdouble(fid,[modelname '.inversion.iscontrol'],self.iscontrol); - writejscellstring(fid,[modelname '.inversion.control_parameters'],self.control_parameters); - writejsdouble(fid,[modelname '.inversion.control_scaling_factors'],self.control_scaling_factors); writejsdouble(fid,[modelname '.inversion.maxsteps'],self.maxsteps); writejsdouble(fid,[modelname '.inversion.maxiter'],self.maxiter); writejsdouble(fid,[modelname '.inversion.dxmin'],self.dxmin); writejsdouble(fid,[modelname '.inversion.gttol'],self.gttol); - writejs2Darray(fid,[modelname '.inversion.cost_functions'],self.cost_functions); - writejs2Darray(fid,[modelname '.inversion.cost_functions_coefficients'],self.cost_functions_coefficients); - writejs1Darray(fid,[modelname '.inversion.min_parameters'],self.min_parameters); - writejs1Darray(fid,[modelname '.inversion.max_parameters'],self.max_parameters); - writejs1Darray(fid,[modelname '.inversion.vx_obs'],self.vx_obs); - writejs1Darray(fid,[modelname '.inversion.vy_obs'],self.vy_obs); - writejs1Darray(fid,[modelname '.inversion.vz_obs'],self.vz_obs); - writejs1Darray(fid,[modelname '.inversion.vel_obs'],self.vel_obs); - writejs1Darray(fid,[modelname '.inversion.thickness_obs'],self.thickness_obs); - writejs1Darray(fid,[modelname '.inversion.surface_obs'],self.surface_obs); end % }}} Index: /issm/trunk-jpl/src/m/classes/autodiff.m =================================================================== --- /issm/trunk-jpl/src/m/classes/autodiff.m (revision 22514) +++ /issm/trunk-jpl/src/m/classes/autodiff.m (revision 22515) @@ -130,7 +130,15 @@ names{i}=indep.name; types(i)=indep.typetoscalar(); + min_parameters(:,i)=indep.min_parameters; + max_parameters(:,i)=indep.max_parameters; + scaling_factors(i)=indep.control_scaling_factor; end WriteData(fid,prefix,'data',names,'name','md.autodiff.independent_object_names','format','StringArray'); WriteData(fid,prefix,'data',types,'name','md.autodiff.independent_object_types','format','IntMat','mattype',3); + WriteData(fid,prefix,'data',min_parameters,'name','md.autodiff.independent_min_parameters','format','DoubleMat','mattype',3); + WriteData(fid,prefix,'data',max_parameters,'name','md.autodiff.independent_max_parameters','format','DoubleMat','mattype',3); + WriteData(fid,prefix,'data',scaling_factors,'name','md.autodiff.independent_scaling_factors','format','DoubleMat','mattype',3); + + end %}}} Index: /issm/trunk-jpl/src/m/classes/independent.m =================================================================== --- /issm/trunk-jpl/src/m/classes/independent.m (revision 22514) +++ /issm/trunk-jpl/src/m/classes/independent.m (revision 22515) @@ -11,4 +11,7 @@ fov_forward_indices = []; nods = 0; + min_parameters = NaN; + max_parameters = NaN; + control_scaling_factor = NaN end methods @@ -45,4 +48,8 @@ end md = checkfield(md,'fieldname',['autodiff.independents{' num2str(i) '}.fov_forward_indices'],'>=',1,'<=',self.nods,'size',[NaN 1]); + md = checkfield(md,'fieldname',['autodiff.independents{' num2str(i) '}.min_parameters'],'size',[md.mesh.numberofvertices 1]); + md = checkfield(md,'fieldname',['autodiff.independents{' num2str(i) '}.max_parameters'],'size',[md.mesh.numberofvertices 1]); + md = checkfield(md,'fieldname',['autodiff.independents{' num2str(i) '}.control_scaling_factors'],'size',[1 1],'>',0,'NaN',1,'Inf',1); + end @@ -54,4 +61,7 @@ fielddisplay(self,'type','type of variable (''vertex'' or ''scalar'')'); fielddisplay(self,'nods','size of dependent variables'); + fielddisplay(self,'min_parameters','absolute minimum acceptable value of the inversed parameter on each vertex'); + fielddisplay(self,'max_parameters','absolute maximum acceptable value of the inversed parameter on each vertex'); + fielddisplay(self,'control_scaling_factor','order of magnitude of each control (useful for multi-parameter optimization)'); if ~isnan(self.fos_forward_index), fielddisplay(self,'fos_forward_index','index for fos_foward driver of ADOLC');