Ice Sheet System Model  4.18
Code documentation
controlm1qn3_core.cpp
Go to the documentation of this file.
1 
5 #include <config.h>
6 #include "./cores.h"
7 #include "../toolkits/toolkits.h"
8 #include "../classes/classes.h"
9 #include "../shared/shared.h"
10 #include "../modules/modules.h"
11 #include "../solutionsequences/solutionsequences.h"
12 
13 #if defined (_HAVE_M1QN3_)
14 /*m1qn3 prototypes {{{*/
15 extern "C" void *ctonbe_; // DIS mode : Conversion
16 extern "C" void *ctcabe_; // DIS mode : Conversion
17 extern "C" void *euclid_; // Scalar product
18 typedef void (*SimulFunc) (long* indic,long* n, double* x, double* pf,double* g,long [],float [],void* dzs);
19 extern "C" void m1qn3_ (void f(long* indic,long* n, double* x, double* pf,double* g,long [],float [],void* dzs),
20  void **, void **, void **,
21  long *, double [], double *, double [], double*, double *,
22  double *, char [], long *, long *, long *, long *, long *, long *, long [], double [], long *,
23  long *, long *, long [], float [],void* );
24 
25 /*Cost function prototype*/
26 void simul(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs);
27 
28 /*Use struct to provide arguments*/
29 typedef struct{
30  FemModel * femmodel;
31  IssmPDouble* Jlist;
32  int M;
33  int N;
34  int* i;
35 } m1qn3_struct;
36 /*}}}*/
37 
38 void controlm1qn3_core(FemModel* femmodel){/*{{{*/
39 
40  /*Intermediaries*/
41  long omode;
42  double f,dxmin,gttol;
43  int maxsteps,maxiter;
44  int intn,numberofvertices,num_controls,num_cost_functions,solution_type;
45  double *scaling_factors = NULL;
46  double *X = NULL;
47  double *G = NULL;
48 
49  /*Recover some parameters*/
50  femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
51  femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
52  femmodel->parameters->FindParam(&num_cost_functions,InversionNumCostFunctionsEnum);
53  femmodel->parameters->FindParam(&maxsteps,InversionMaxstepsEnum);
54  femmodel->parameters->FindParam(&maxiter,InversionMaxiterEnum);
55  femmodel->parameters->FindParamAndMakePassive(&dxmin,InversionDxminEnum);
56  femmodel->parameters->FindParamAndMakePassive(&gttol,InversionGttolEnum);
57  femmodel->parameters->FindParamAndMakePassive(&scaling_factors,NULL,InversionControlScalingFactorsEnum);
58  femmodel->parameters->SetParam(false,SaveResultsEnum);
59  numberofvertices=femmodel->vertices->NumberOfVertices();
60 
61  /*Initialize M1QN3 parameters*/
62  if(VerboseControl())_printf0_(" Initialize M1QN3 parameters\n");
63  SimulFunc costfuncion = &simul; /*Cost function address*/
64  void** prosca = &euclid_; /*Dot product function (euclid is the default)*/
65  char normtype[] = "dfn"; /*Norm type: dfn = scalar product defined by prosca*/
66  long izs[5]; /*Arrays used by m1qn3 subroutines*/
67  long iz[5]; /*Integer m1qn3 working array of size 5*/
68  float rzs[1]; /*Arrays used by m1qn3 subroutines*/
69  long impres = 0; /*verbosity level*/
70  long imode[3] = {0}; /*scaling and starting mode, 0 by default*/
71  long indic = 4; /*compute f and g*/
72  long reverse = 0; /*reverse or direct mode*/
73  long io = 6; /*Channel number for the output*/
74 
75  /*Optimization criterions*/
76  long niter = long(maxsteps); /*Maximum number of iterations*/
77  long nsim = long(maxiter);/*Maximum number of function calls*/
78 
79  /*Get initial guess*/
80  GetPassiveVectorFromControlInputsx(&X,&intn,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value");
81  _assert_(intn==numberofvertices*num_controls);
82 
83  /*Get problem dimension and initialize gradient and initial guess*/
84  long n = long(intn);
85  G = xNew<double>(n);
86 
87  /*Scale control for M1QN3*/
88  for(int c=0;c<num_controls;c++){
89  for(int i=0;i<numberofvertices;i++){
90  int index = numberofvertices*c+i;
91  X[index] = X[index]/scaling_factors[c];
92  }
93  }
94 
95  /*Allocate m1qn3 working arrays (see documentation)*/
96  long m = 100;
97  long ndz = 4*n+m*(2*n+1);
98  double* dz = xNew<double>(ndz);
99 
100  if(VerboseControl())_printf0_(" Computing initial solution\n");
101  _printf0_("\n");
102  _printf0_("Cost function f(x) | Gradient norm |g(x)| | List of contributions\n");
103  _printf0_("____________________________________________________________________\n");
104 
105  /*Prepare structure for m1qn3*/
106  m1qn3_struct mystruct;
107  mystruct.femmodel = femmodel;
108  mystruct.M = maxiter;
109  mystruct.N = num_cost_functions+1;
110  mystruct.Jlist = xNewZeroInit<IssmPDouble>(mystruct.M*mystruct.N);
111  mystruct.i = xNewZeroInit<int>(1);
112 
113  /*Initialize Gradient and cost function of M1QN3*/
114  indic = 4; /*gradient required*/
115  simul(&indic,&n,X,&f,G,izs,rzs,(void*)&mystruct);
116 
117  /*Estimation of the expected decrease in f during the first iteration*/
118  double df1=f;
119 
120  /*Call M1QN3 solver*/
121  m1qn3_(costfuncion,prosca,&ctonbe_,&ctcabe_,
122  &n,X,&f,G,&dxmin,&df1,
123  &gttol,normtype,&impres,&io,imode,&omode,&niter,&nsim,iz,dz,&ndz,
124  &reverse,&indic,izs,rzs,(void*)&mystruct);
125 
126  /*Print exit flag*/
127  switch(int(omode)){
128  case 0: _printf0_(" Stop requested (indic = 0)\n"); break;
129  case 1: _printf0_(" Convergence reached (gradient satisfies stopping criterion)\n"); break;
130  case 2: _printf0_(" Bad initialization\n"); break;
131  case 3: _printf0_(" Line search failure\n"); break;
132  case 4: _printf0_(" Maximum number of iterations exceeded\n");break;
133  case 5: _printf0_(" Maximum number of function calls exceeded\n"); break;
134  case 6: _printf0_(" stopped on dxmin during line search\n"); break;
135  case 7: _printf0_(" <g,d> > 0 or <y,s> <0\n"); break;
136  default: _printf0_(" Unknown end condition\n");
137  }
138 
139  /*Constrain solution vector*/
140  double *XL = NULL;
141  double *XU = NULL;
142  GetPassiveVectorFromControlInputsx(&XL,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
143  GetPassiveVectorFromControlInputsx(&XU,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
144  for(int c=0;c<num_controls;c++){
145  for(int i=0;i<numberofvertices;i++){
146  int index = numberofvertices*c+i;
147  X[index] = X[index]*scaling_factors[c];
148  if(X[index]>XU[index]) X[index]=XU[index];
149  if(X[index]<XL[index]) X[index]=XL[index];
150  }
151  }
152 
153  /*Set X as our new control (need to recast)*/
154  #ifdef _HAVE_AD_
155  IssmDouble* aX=xNew<IssmDouble>(intn);
156  IssmDouble* aG=xNew<IssmDouble>(intn);
157  for(int i=0;i<intn;i++) {
158  aX[i] = reCast<IssmDouble>(X[i]);
159  aG[i] = reCast<IssmDouble>(G[i]);
160  }
161  ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,aG);
162  SetControlInputsFromVectorx(femmodel,aX);
163  xDelete(aX);
164  xDelete(aG);
165  #else
166  SetControlInputsFromVectorx(femmodel,X);
167  ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,G);
168  #endif
169 
170  femmodel->OutputControlsx(&femmodel->results);
171  femmodel->results->AddObject(new GenericExternalResult<double*>(femmodel->results->Size()+1,JEnum,mystruct.Jlist,(*mystruct.i),mystruct.N,0,0));
172 
173  /*Finalize*/
174  if(VerboseControl()) _printf0_(" preparing final solution\n");
175  femmodel->parameters->SetParam(true,SaveResultsEnum);
176  void (*solutioncore)(FemModel*)=NULL;
177  CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type);
178  solutioncore(femmodel);
179 
180  /*Clean-up and return*/
181  xDelete<double>(G);
182  xDelete<double>(X);
183  xDelete<double>(dz);
184  xDelete<double>(XU);
185  xDelete<double>(XL);
186  xDelete<double>(scaling_factors);
187  xDelete<double>(mystruct.Jlist);
188  xDelete<int>(mystruct.i);
189 }/*}}}*/
190 void simul(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs){/*{{{*/
191 
192  /*Recover Arguments*/
193  m1qn3_struct *input_struct = (m1qn3_struct*)dzs;
194  FemModel *femmodel = input_struct->femmodel;
195  IssmPDouble *Jlist = input_struct->Jlist;
196  int JlistM = input_struct->M;
197  int JlistN = input_struct->N;
198  int *Jlisti = input_struct->i;
199 
200  /*Recover some parameters*/
201  int num_responses,num_controls,numberofvertices,solution_type;
202  double* scaling_factors = NULL;
203  femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
204  femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
205  femmodel->parameters->FindParamAndMakePassive(&scaling_factors,NULL,InversionControlScalingFactorsEnum);
206  femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
207  numberofvertices=femmodel->vertices->NumberOfVertices();
208 
209  /*Constrain input vector and update controls*/
210  double *XL = NULL;
211  double *XU = NULL;
212  GetPassiveVectorFromControlInputsx(&XL,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
213  GetPassiveVectorFromControlInputsx(&XU,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
214  for(int c=0;c<num_controls;c++){
215  for(int i=0;i<numberofvertices;i++){
216  int index = numberofvertices*c+i;
217  X[index] = X[index]*scaling_factors[c];
218  if(X[index]>XU[index]) X[index]=XU[index];
219  if(X[index]<XL[index]) X[index]=XL[index];
220  }
221  }
222  #ifdef _HAVE_AD_
223  IssmDouble* aX=xNew<IssmDouble>(*n);
224  for(int i=0;i<*n;i++) aX[i] = reCast<IssmDouble>(X[i]);
225  SetControlInputsFromVectorx(femmodel,aX);
226  xDelete(aX);
227  #else
228  SetControlInputsFromVectorx(femmodel,X);
229  #endif
230 
231  /*Compute solution and adjoint*/
232  void (*solutioncore)(FemModel*)=NULL;
233  void (*adjointcore)(FemModel*)=NULL;
234  CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type);
235  solutioncore(femmodel);
236 
237  /*Check size of Jlist to avoid crashes*/
238  _assert_((*Jlisti)<JlistM);
239  _assert_(JlistN==num_responses+1);
240 
241  /*Compute objective function*/
242  IssmDouble* Jtemp = NULL;
243  IssmDouble J;
244  femmodel->CostFunctionx(&J,&Jtemp,NULL);
245  *pf = reCast<double>(J);
246  _printf0_("f(x) = "<<setw(12)<<setprecision(7)<<*pf<<" | ");
247 
248  /*Record cost function values and delete Jtemp*/
249  for(int i=0;i<num_responses;i++) Jlist[(*Jlisti)*JlistN+i] = reCast<double>(Jtemp[i]);
250  Jlist[(*Jlisti)*JlistN+num_responses] = *pf;
251  xDelete<IssmDouble>(Jtemp);
252 
253  if(*indic==0){
254  /*dry run, no gradient required*/
255 
256  /*Retrieve objective functions independently*/
257  _printf0_(" N/A |\n");
258  for(int i=0;i<num_responses;i++) _printf0_(" "<<setw(12)<<setprecision(7)<<Jlist[(*Jlisti)*JlistN+i]);
259  _printf0_("\n");
260 
261  *Jlisti = (*Jlisti) +1;
262  xDelete<double>(XU);
263  xDelete<double>(XL);
264  return;
265  }
266 
267  /*Compute Adjoint*/
268  AdjointCorePointerFromSolutionEnum(&adjointcore,solution_type);
269  adjointcore(femmodel);
270 
271  /*Compute gradient*/
272  IssmDouble* G2 = NULL;
273  Gradjx(&G2,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
274  for(long i=0;i<*n;i++) G[i] = -reCast<double>(G2[i]);
275  xDelete<IssmDouble>(G2);
276 
277  /*Constrain Gradient*/
278  IssmDouble Gnorm = 0.;
279  for(int c=0;c<num_controls;c++){
280  for(int i=0;i<numberofvertices;i++){
281  int index = numberofvertices*c+i;
282  if(X[index]>=XU[index]) G[index]=0.;
283  if(X[index]<=XL[index]) G[index]=0.;
284  G[index] = G[index]*scaling_factors[c];
285  X[index] = X[index]/scaling_factors[c];
286  Gnorm += G[index]*G[index];
287  }
288  }
289  Gnorm = sqrt(Gnorm);
290 
291  /*Print info*/
292  _printf0_(" "<<setw(12)<<setprecision(7)<<Gnorm<<" |");
293  for(int i=0;i<num_responses;i++) _printf0_(" "<<setw(12)<<setprecision(7)<<Jlist[(*Jlisti)*JlistN+i]);
294  _printf0_("\n");
295 
296  /*Clean-up and return*/
297  *Jlisti = (*Jlisti) +1;
298  xDelete<double>(XU);
299  xDelete<double>(XL);
300  xDelete<double>(scaling_factors);
301 }/*}}}*/
302 
303 #else
304 void controlm1qn3_core(FemModel* femmodel){_error_("M1QN3 not installed");}
305 #endif //_HAVE_M1QN3_
DataSet::Size
int Size()
Definition: DataSet.cpp:399
SaveResultsEnum
@ SaveResultsEnum
Definition: EnumDefinitions.h:302
ControlInputSetGradientx
void ControlInputSetGradientx(Elements *elements, Nodes *nodes, Vertices *vertices, Loads *loads, Materials *materials, Parameters *parameters, IssmDouble *gradient)
Definition: ControlInputSetGradientx.cpp:9
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
FemModel::vertices
Vertices * vertices
Definition: FemModel.h:49
GetPassiveVectorFromControlInputsx
void GetPassiveVectorFromControlInputsx(IssmPDouble **pvector, int *pN, Elements *elements, Nodes *nodes, Vertices *vertices, Loads *loads, Materials *materials, Parameters *parameters, const char *data)
Definition: GetVectorFromControlInputsx.cpp:127
IssmDouble
double IssmDouble
Definition: types.h:37
InversionNumControlParametersEnum
@ InversionNumControlParametersEnum
Definition: EnumDefinitions.h:223
_printf0_
#define _printf0_(StreamArgs)
Definition: Print.h:29
DataSet::AddObject
int AddObject(Object *object)
Definition: DataSet.cpp:252
AdjointCorePointerFromSolutionEnum
void AdjointCorePointerFromSolutionEnum(void(**padjointcore)(FemModel *), int solutiontype)
Definition: AdjointCorePointerFromSolutionEnum.cpp:18
cores.h
FemModel::parameters
Parameters * parameters
Definition: FemModel.h:46
CorePointerFromSolutionEnum
void CorePointerFromSolutionEnum(void(**psolutioncore)(FemModel *), Parameters *parameters, int solutiontype)
Definition: CorePointerFromSolutionEnum.cpp:18
FemModel::CostFunctionx
void CostFunctionx(IssmDouble *pJ, IssmDouble **pJlist, int *pn)
Definition: FemModel.cpp:1062
JEnum
@ JEnum
Definition: EnumDefinitions.h:1134
FemModel::results
Results * results
Definition: FemModel.h:48
Vertices::NumberOfVertices
int NumberOfVertices(void)
Definition: Vertices.cpp:255
InversionMaxstepsEnum
@ InversionMaxstepsEnum
Definition: EnumDefinitions.h:221
FemModel::OutputControlsx
void OutputControlsx(Results **presults)
Definition: FemModel.cpp:2171
Parameters::SetParam
void SetParam(bool boolean, int enum_type)
Definition: Parameters.cpp:441
FemModel::nodes
Nodes * nodes
Definition: FemModel.h:56
xDelete
void xDelete(T **&aT_pp)
Definition: MemOps.h:97
FemModel::materials
Materials * materials
Definition: FemModel.h:45
Gradjx
void Gradjx(Vector< IssmDouble > **pgradient, IssmDouble **pnorm_list, Elements *elements, Nodes *nodes, Vertices *vertices, Loads *loads, Materials *materials, Parameters *parameters)
Definition: Gradjx.cpp:9
SolutionTypeEnum
@ SolutionTypeEnum
Definition: EnumDefinitions.h:398
FemModel::loads
Loads * loads
Definition: FemModel.h:54
Parameters::FindParamAndMakePassive
void FindParamAndMakePassive(IssmPDouble *pscalar, int enum_type)
Definition: Parameters.cpp:379
controlm1qn3_core
void controlm1qn3_core(FemModel *femmodel)
Definition: controlm1qn3_core.cpp:304
FemModel::elements
Elements * elements
Definition: FemModel.h:44
FemModel
Definition: FemModel.h:31
GenericExternalResult
Definition: GenericExternalResult.h:21
InversionControlScalingFactorsEnum
@ InversionControlScalingFactorsEnum
Definition: EnumDefinitions.h:210
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
InversionNumCostFunctionsEnum
@ InversionNumCostFunctionsEnum
Definition: EnumDefinitions.h:224
InversionMaxiterEnum
@ InversionMaxiterEnum
Definition: EnumDefinitions.h:219
Parameters::FindParam
void FindParam(bool *pinteger, int enum_type)
Definition: Parameters.cpp:262
InversionDxminEnum
@ InversionDxminEnum
Definition: EnumDefinitions.h:212
SetControlInputsFromVectorx
void SetControlInputsFromVectorx(FemModel *femmodel, IssmDouble *vector)
Definition: SetControlInputsFromVectorx.cpp:9
InversionGttolEnum
@ InversionGttolEnum
Definition: EnumDefinitions.h:216
IssmPDouble
IssmDouble IssmPDouble
Definition: types.h:38
VerboseControl
bool VerboseControl(void)
Definition: Verbosity.cpp:27
femmodel
FemModel * femmodel
Definition: esmfbinders.cpp:16
Generated on Thu Jul 2 2020 08:09:18 for Ice Sheet System Model by   doxygen 1.8.19