Changeset 18617

Timestamp:

10/10/14 15:05:23 (10 years ago)

Author:

Mathieu Morlighem

Message:

CHG: use scaling factor for inversion validation

Location:

issm/trunk-jpl/src

Files:

• c/cores/controlvalidation_core.cpp (modified) (6 diffs)
• c/modules/ModelProcessorx/Control/CreateParametersControl.cpp (modified) (1 diff)
• m/classes/inversionvalidation.m (modified) (5 diffs)

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

@@ -15 +15 @@
         IssmDouble  j0,j,yts;
         IssmDouble  Ialpha,exponent,alpha;
+        IssmDouble* scaling_factors = NULL;
         IssmDouble* jlist = NULL;
         IssmDouble *G = NULL;
@@ -29 +30 @@
         femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
         femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
+        femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum);
 
         /*Get initial guess*/
@@ -51 +53 @@
         if(VerboseControl()) _printf0_("   Compute Initial cost function\n");
         femmodel->CostFunctionx(&j0,&jlist,NULL);
-        _printf0_("Initial J(x+dk)      |  List of contributions\n");
-        _printf0_("_____________________________________________\n");
-        _printf0_("J(x) = "<<setw(12)<<setprecision(7)<<j0<<"  |  ");
-        for(int i=0;i<num_responses;i++) _printf0_(" "<<setw(12)<<setprecision(7)<<jlist[i]);
-        _printf0_("\n");
+        _printf0_("Initial cost function J(x) = "<<setw(12)<<setprecision(7)<<j0<<"\n");
         xDelete<IssmDouble>(jlist);
 
@@ -62 +60 @@
         for(int i=0;i<n;i++) G[i] = -G[i];
 
-        /*Range of tests*/
-        IssmDouble exp0 = 0.;
-        IssmDouble incr = -0.2;
-        IssmDouble exp1 = -18.;
-        int        num  = reCast<int,IssmDouble>((exp1-exp0)/incr);
-
         /*Allocate output*/
+        int num = 26;
         IssmDouble* output = xNew<IssmDouble>(2*num);
 
@@ -76 +69 @@
         for(int m=0;m<num;m++){
 
-                /*Calculate alpha = 10^-exponent*/
-                exponent = exp0+m*incr;
-                alpha    = pow(10.,exponent);
-
                 /*Create new vector*/
+                alpha    = pow(2.,-m);
                 for(int i=0;i<n;i++) X[i] = X0[i] + alpha;
 
@@ -89 +79 @@
 
                 IssmDouble Den = 0.;
-                for(int i=0;i<n;i++) Den += alpha* G[i] * 1.;
+                for(int i=0;i<n;i++) Den += alpha* G[i] * scaling_factors[0];
                 Ialpha = fabs((j - j0)/Den - 1.);
 

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

@@ -78 +78 @@
                                 break;
                         case 3:/*Validation*/
+                                iomodel->FetchData(&control_scaling_factors,NULL,NULL,InversionControlScalingFactorsEnum);
+                                parameters->AddObject(new DoubleVecParam(InversionControlScalingFactorsEnum,control_scaling_factors,num_controls));
                                 break;
                         default:

issm/trunk-jpl/src/m/classes/inversionvalidation.m ¶

@@ -9 +9 @@
                 incomplete_adjoint          = 0
                 control_parameters          = NaN
+                control_scaling_factors     = NaN
                 cost_functions              = NaN
                 cost_functions_coefficients = NaN
@@ -40 +41 @@
                         self.control_parameters={'FrictionCoefficient'};
 
+                        %Scaling factor for each control
+                        self.control_scaling_factors=1;
+
                         %several responses can be used:
                         self.cost_functions=101;
@@ -56 +60 @@
                                 {'BalancethicknessThickeningRate' 'FrictionCoefficient' 'MaterialsRheologyBbar' 'DamageDbar' 'Vx' 'Vy' 'Thickness',...
                                 'BalancethicknessOmega' 'BalancethicknessApparentMassbalance'});
+                        md = checkfield(md,'fieldname','inversion.control_scaling_factors','size',[1 num_controls],'>',0,'NaN',1);
                         md = checkfield(md,'fieldname','inversion.cost_functions','size',[1 num_costfunc],'values',[101:105 201 501:506 601:604]);
                         md = checkfield(md,'fieldname','inversion.cost_functions_coefficients','size',[md.mesh.numberofvertices num_costfunc],'>=',0);
@@ -77 +82 @@
                         fielddisplay(obj,'incomplete_adjoint','1: linear viscosity, 0: non-linear viscosity');
                         fielddisplay(obj,'control_parameters','ex: {''FrictionCoefficient''}, or {''MaterialsRheologyBbar''}');
+                        fielddisplay(obj,'control_scaling_factors','order of magnitude of each control (useful for multi-parameter optimization)');
                         fielddisplay(obj,'cost_functions','indicate the type of response for each optimization step');
                         fielddisplay(obj,'cost_functions_coefficients','cost_functions_coefficients applied to the misfit of each vertex and for each control_parameter');
@@ -105 +111 @@
                         if ~obj.iscontrol, return; end
                         WriteData(fid,'object',obj,'class','inversion','fieldname','incomplete_adjoint','format','Boolean');
+                        WriteData(fid,'object',obj,'class','inversion','fieldname','control_scaling_factors','format','DoubleMat','mattype',3);
                         WriteData(fid,'object',obj,'class','inversion','fieldname','cost_functions_coefficients','format','DoubleMat','mattype',1);
                         WriteData(fid,'object',obj,'class','inversion','fieldname','min_parameters','format','DoubleMat','mattype',3);