source: issm/oecreview/Archive/18296-19100/ISSM-18615-18616.diff

Last change on this file was 19102, checked in by Mathieu Morlighem, 10 years ago

NEW: added 18296-19100
File size: 5.8 KB
RevLine 
[19102]1Index: ../trunk-jpl/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp
2===================================================================
3--- ../trunk-jpl/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp (revision 18615)
4+++ ../trunk-jpl/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp (revision 18616)
5@@ -12,13 +12,14 @@
6 bool control_analysis;
7 int inversiontype;
8 int nsteps;
9- int num_control_type;
10- int num_cm_responses;
11+ int num_controls;
12+ int num_costfunc;
13 int *control_type = NULL;
14 int *maxiter = NULL;
15 int *cm_responses = NULL;
16 IssmDouble *cm_jump = NULL;
17 IssmDouble *optscal = NULL;
18+ IssmDouble *control_scaling_factors = NULL;
19
20 /*retrieve some parameters: */
21 iomodel->Constant(&control_analysis,InversionIscontrolEnum);
22@@ -40,10 +41,10 @@
23 }
24
25 /*Now, recover fit, optscal and maxiter as vectors: */
26- iomodel->FetchData(&control_type,NULL,&num_control_type,InversionControlParametersEnum);
27- iomodel->FetchData(&cm_responses,NULL,&num_cm_responses,InversionCostFunctionsEnum);
28- parameters->AddObject(new IntVecParam(InversionControlParametersEnum,control_type,num_control_type));
29- parameters->AddObject(new IntVecParam(InversionCostFunctionsEnum,cm_responses,num_cm_responses));
30+ iomodel->FetchData(&control_type,NULL,&num_controls,InversionControlParametersEnum);
31+ iomodel->FetchData(&cm_responses,NULL,&num_costfunc,InversionCostFunctionsEnum);
32+ parameters->AddObject(new IntVecParam(InversionControlParametersEnum,control_type,num_controls));
33+ parameters->AddObject(new IntVecParam(InversionCostFunctionsEnum,cm_responses,num_costfunc));
34
35 /*Inversion type specifics*/
36 switch(inversiontype){
37@@ -53,7 +54,7 @@
38 iomodel->FetchData(&cm_jump,&nsteps,NULL,InversionStepThresholdEnum);
39 iomodel->FetchData(&optscal,NULL,NULL,InversionGradientScalingEnum);
40 iomodel->FetchData(&maxiter,NULL,NULL,InversionMaxiterPerStepEnum);
41- parameters->AddObject(new DoubleMatParam(InversionGradientScalingEnum,optscal,nsteps,num_control_type));
42+ parameters->AddObject(new DoubleMatParam(InversionGradientScalingEnum,optscal,nsteps,num_controls));
43 parameters->AddObject(new DoubleVecParam(InversionStepThresholdEnum,cm_jump,nsteps));
44 parameters->AddObject(new IntVecParam(InversionMaxiterPerStepEnum,maxiter,nsteps));
45 break;
46@@ -72,6 +73,8 @@
47 parameters->AddObject(iomodel->CopyConstantObject(InversionGttolEnum));
48 parameters->AddObject(iomodel->CopyConstantObject(InversionMaxstepsEnum));
49 parameters->AddObject(iomodel->CopyConstantObject(InversionMaxiterEnum));
50+ iomodel->FetchData(&control_scaling_factors,NULL,NULL,InversionControlScalingFactorsEnum);
51+ parameters->AddObject(new DoubleVecParam(InversionControlScalingFactorsEnum,control_scaling_factors,num_controls));
52 break;
53 case 3:/*Validation*/
54 break;
55Index: ../trunk-jpl/src/c/cores/controlm1qn3_core.cpp
56===================================================================
57--- ../trunk-jpl/src/c/cores/controlm1qn3_core.cpp (revision 18615)
58+++ ../trunk-jpl/src/c/cores/controlm1qn3_core.cpp (revision 18616)
59@@ -32,6 +32,7 @@
60 double f,dxmin,gttol;
61 int maxsteps,maxiter;
62 int intn,num_controls,solution_type;
63+ IssmDouble *scaling_factors = NULL;
64 IssmDouble *X = NULL;
65 IssmDouble *G = NULL;
66
67@@ -42,6 +43,7 @@
68 femmodel->parameters->FindParam(&maxiter,InversionMaxiterEnum);
69 femmodel->parameters->FindParam(&dxmin,InversionDxminEnum);
70 femmodel->parameters->FindParam(&gttol,InversionGttolEnum);
71+ femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum);
72 femmodel->parameters->SetParam(false,SaveResultsEnum);
73
74 /*Initialize M1QN3 parameters*/
75@@ -73,6 +75,12 @@
76 long n = long(intn);
77 G = xNew<double>(n);
78
79+ /*Scale control for M1QN3*/
80+ if(num_controls!=1) _error_("not supported yet...");
81+ for(long i=0;i<n;i++){
82+ X[i] = X[i]/scaling_factors[0];
83+ }
84+
85 /*Allocate m1qn3 working arrays (see doc)*/
86 long m = 100;
87 long ndz = 4*n+m*(2*n+1);
88@@ -131,15 +139,21 @@
89 FemModel *femmodel = (FemModel*)dzs;
90
91 /*Recover number of cost functions responses*/
92- int num_responses;
93+ int num_responses;
94+ int num_controls;
95+ IssmDouble* scaling_factors = NULL;
96 femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
97+ femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
98+ femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum);
99
100 /*Constrain input vector*/
101 IssmDouble *XL = NULL;
102 IssmDouble *XU = NULL;
103 GetVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
104 GetVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
105+ if(num_controls!=1) _error_("not supported yet");
106 for(long i=0;i<*n;i++){
107+ X[i] = X[i]*scaling_factors[0];
108 if(X[i]>XU[i]) X[i]=XU[i];
109 if(X[i]<XL[i]) X[i]=XL[i];
110 }
111@@ -161,8 +175,6 @@
112 femmodel->CostFunctionx(pf,&Jlist,NULL);
113 _printf0_("f(x) = "<<setw(12)<<setprecision(7)<<*pf<<" | ");
114
115-
116-
117 if(indic==0){
118 /*dry run, no gradient required*/
119
120@@ -189,9 +201,12 @@
121
122 /*Constrain Gradient*/
123 IssmDouble Gnorm = 0.;
124+ if(num_controls!=1) _error_("not supported yet");
125 for(long i=0;i<*n;i++){
126 if(X[i]>=XU[i]) G[i]=0.;
127 if(X[i]<=XL[i]) G[i]=0.;
128+ G[i] = G[i]*scaling_factors[0];
129+ X[i] = X[i]/scaling_factors[0];
130 Gnorm += G[i]*G[i];
131 }
132 Gnorm = sqrt(Gnorm);
133@@ -202,7 +217,7 @@
134 _printf0_("\n");
135
136 /*Clean-up and return*/
137- xDelete<IssmDouble>(Jlist);
138+ xDelete<IssmDouble>(Jlist);
139 xDelete<IssmDouble>(XU);
140 xDelete<IssmDouble>(XL);
141 }
Note: See TracBrowser for help on using the repository browser.