source: issm/trunk/src/m/classes/taoinversion.m@ 17989

Last change on this file since 17989 was 17989, checked in by Mathieu Morlighem, 11 years ago

merged trunk-jpl and trunk for revision 17986
File size: 9.7 KB
Line 
1%TAOINVERSION class definition
2%
3% Usage:
4% taoinversion=taoinversion();
5
6classdef taoinversion
7 properties (SetAccess=public)
8 iscontrol = 0
9 incomplete_adjoint = 0
10 control_parameters = NaN
11 maxsteps = 0
12 maxiter = 0
13 fatol = 0
14 frtol = 0
15 gatol = 0
16 grtol = 0
17 gttol = 0
18 algorithm = ''
19 cost_functions = NaN
20 cost_functions_coefficients = NaN
21 min_parameters = NaN
22 max_parameters = NaN
23 vx_obs = NaN
24 vy_obs = NaN
25 vz_obs = NaN
26 vel_obs = NaN
27 thickness_obs = NaN
28 end
29 methods
30 function obj = taoinversion(varargin) % {{{
31 switch nargin
32 case 0
33 obj=setdefaultparameters(obj);
34 case 1
35 obj=structtoobj(taoinversion(),varargin{1});
36 otherwise
37 error('constructor not supported');
38 end
39 end % }}}
40 function obj = setdefaultparameters(obj) % {{{
41
42 %default is incomplete adjoint for now
43 obj.incomplete_adjoint=1;
44
45 %parameter to be inferred by control methods (only
46 %drag and B are supported yet)
47 obj.control_parameters={'FrictionCoefficient'};
48
49 %number of iterations and steps
50 obj.maxsteps=20;
51 obj.maxiter =30;
52
53 %default tolerances
54 obj.fatol = 0;
55 obj.frtol = 0;
56 obj.gatol = 0;
57 obj.grtol = 0;
58 obj.gttol = 1e-4;
59
60 %minimization algorithm
61 obj.algorithm = 'tao_blmvm';
62
63 %several responses can be used:
64 obj.cost_functions=101;
65
66 end % }}}
67 function md = checkconsistency(obj,md,solution,analyses) % {{{
68
69 %Early return
70 if ~obj.iscontrol, return; end
71
72 if ~IssmConfig('_HAVE_TAO_'),
73 md = checkmessage(md,['TAO has not been installed, ISSM needs to be reconfigured and recompiled with TAO']);
74 end
75
76 num_controls=numel(md.inversion.control_parameters);
77 num_costfunc=size(md.inversion.cost_functions,2);
78
79 md = checkfield(md,'fieldname','inversion.iscontrol','values',[0 1]);
80 md = checkfield(md,'fieldname','inversion.incomplete_adjoint','values',[0 1]);
81 md = checkfield(md,'fieldname','inversion.control_parameters','cell',1,'values',...
82 {'BalancethicknessThickeningRate' 'FrictionCoefficient' 'MaterialsRheologyBbar' 'MaterialsRheologyZbar' 'Vx' 'Vy' 'Thickness'});
83 md = checkfield(md,'fieldname','inversion.maxsteps','numel',1,'>=',0);
84 md = checkfield(md,'fieldname','inversion.maxiter','numel',1,'>=',0);
85 md = checkfield(md,'fieldname','inversion.fatol','numel',1,'>=',0);
86 md = checkfield(md,'fieldname','inversion.frtol','numel',1,'>=',0);
87 md = checkfield(md,'fieldname','inversion.gatol','numel',1,'>=',0);
88 md = checkfield(md,'fieldname','inversion.grtol','numel',1,'>=',0);
89 md = checkfield(md,'fieldname','inversion.gttol','numel',1,'>=',0);
90 md = checkfield(md,'fieldname','inversion.algorithm','values',{'tao_blmvm','tao_cg','tao_lmvm'});
91 md = checkfield(md,'fieldname','inversion.cost_functions','size',[1 num_costfunc],'values',[101:105 201 501:506]);
92 md = checkfield(md,'fieldname','inversion.cost_functions_coefficients','size',[md.mesh.numberofvertices num_costfunc],'>=',0);
93 md = checkfield(md,'fieldname','inversion.min_parameters','size',[md.mesh.numberofvertices num_controls]);
94 md = checkfield(md,'fieldname','inversion.max_parameters','size',[md.mesh.numberofvertices num_controls]);
95
96 if solution==BalancethicknessSolutionEnum()
97 md = checkfield(md,'fieldname','inversion.thickness_obs','size',[md.mesh.numberofvertices 1],'NaN',1);
98 elseif solution==BalancethicknessSoftSolutionEnum()
99 md = checkfield(md,'fieldname','inversion.thickness_obs','size',[md.mesh.numberofvertices 1],'NaN',1);
100 else
101 md = checkfield(md,'fieldname','inversion.vx_obs','size',[md.mesh.numberofvertices 1],'NaN',1);
102 md = checkfield(md,'fieldname','inversion.vy_obs','size',[md.mesh.numberofvertices 1],'NaN',1);
103 end
104 end % }}}
105 function disp(obj) % {{{
106 disp(sprintf(' taoinversion parameters:'));
107 fielddisplay(obj,'iscontrol','is inversion activated?');
108 fielddisplay(obj,'incomplete_adjoint','1: linear viscosity, 0: non-linear viscosity');
109 fielddisplay(obj,'control_parameters','ex: {''FrictionCoefficient''}, or {''MaterialsRheologyBbar''}');
110 fielddisplay(obj,'maxsteps','maximum number of iterations (gradient computation)');
111 fielddisplay(obj,'maxiter','maximum number of Function evaluation (forward run)');
112 fielddisplay(obj,'fatol','convergence criterion: f(X)-f(X*) (X: current iteration, X*: "true" solution, f: cost function)');
113 fielddisplay(obj,'frtol','convergence criterion: |f(X)-f(X*)|/|f(X*)|');
114 fielddisplay(obj,'gatol','convergence criterion: ||g(X)|| (g: gradient of the cost function)');
115 fielddisplay(obj,'grtol','convergence criterion: ||g(X)||/|f(X)|');
116 fielddisplay(obj,'gttol','convergence criterion: ||g(X)||/||g(X0)|| (g(X0): gradient at initial guess X0)');
117 fielddisplay(obj,'algorithm','minimization algorithm: ''tao_blmvm'', ''tao_cg'', ''tao_lmvm''');
118 fielddisplay(obj,'cost_functions','indicate the type of response for each optimization step');
119 fielddisplay(obj,'cost_functions_coefficients','cost_functions_coefficients applied to the misfit of each vertex and for each control_parameter');
120 fielddisplay(obj,'min_parameters','absolute minimum acceptable value of the inversed parameter on each vertex');
121 fielddisplay(obj,'max_parameters','absolute maximum acceptable value of the inversed parameter on each vertex');
122 fielddisplay(obj,'vx_obs','observed velocity x component [m/yr]');
123 fielddisplay(obj,'vy_obs','observed velocity y component [m/yr]');
124 fielddisplay(obj,'vel_obs','observed velocity magnitude [m/yr]');
125 fielddisplay(obj,'thickness_obs','observed thickness [m]');
126 disp('Available cost functions:');
127 disp(' 101: SurfaceAbsVelMisfit');
128 disp(' 102: SurfaceRelVelMisfit');
129 disp(' 103: SurfaceLogVelMisfit');
130 disp(' 104: SurfaceLogVxVyMisfit');
131 disp(' 105: SurfaceAverageVelMisfit');
132 disp(' 201: ThicknessAbsMisfit');
133 disp(' 501: DragCoefficientAbsGradient');
134 disp(' 502: RheologyBbarAbsGradient');
135 disp(' 503: ThicknessAbsGradient');
136 end % }}}
137 function marshall(obj,md,fid) % {{{
138
139 yts=365.0*24.0*3600.0;
140
141 WriteData(fid,'object',obj,'class','inversion','fieldname','iscontrol','format','Boolean');
142 WriteData(fid,'enum',InversionTypeEnum(),'data',1,'format','Integer');
143 if ~obj.iscontrol, return; end
144 WriteData(fid,'object',obj,'class','inversion','fieldname','incomplete_adjoint','format','Boolean');
145 WriteData(fid,'object',obj,'class','inversion','fieldname','maxsteps','format','Integer');
146 WriteData(fid,'object',obj,'class','inversion','fieldname','maxiter','format','Integer');
147 WriteData(fid,'object',obj,'class','inversion','fieldname','fatol','format','Double');
148 WriteData(fid,'object',obj,'class','inversion','fieldname','frtol','format','Double');
149 WriteData(fid,'object',obj,'class','inversion','fieldname','gatol','format','Double');
150 WriteData(fid,'object',obj,'class','inversion','fieldname','grtol','format','Double');
151 WriteData(fid,'object',obj,'class','inversion','fieldname','gttol','format','Double');
152 WriteData(fid,'object',obj,'class','inversion','fieldname','algorithm','format','String');
153 WriteData(fid,'object',obj,'class','inversion','fieldname','cost_functions_coefficients','format','DoubleMat','mattype',1);
154 WriteData(fid,'object',obj,'class','inversion','fieldname','min_parameters','format','DoubleMat','mattype',3);
155 WriteData(fid,'object',obj,'class','inversion','fieldname','max_parameters','format','DoubleMat','mattype',3);
156 WriteData(fid,'object',obj,'class','inversion','fieldname','vx_obs','format','DoubleMat','mattype',1,'scale',1./yts);
157 WriteData(fid,'object',obj,'class','inversion','fieldname','vy_obs','format','DoubleMat','mattype',1,'scale',1./yts);
158 WriteData(fid,'object',obj,'class','inversion','fieldname','vz_obs','format','DoubleMat','mattype',1,'scale',1./yts);
159 WriteData(fid,'object',obj,'class','inversion','fieldname','thickness_obs','format','DoubleMat','mattype',1);
160
161 %process control parameters
162 num_control_parameters=numel(obj.control_parameters);
163 data=zeros(1,num_control_parameters);
164 for i=1:num_control_parameters,
165 data(i)=StringToEnum(obj.control_parameters{i});
166 end
167 WriteData(fid,'data',data,'enum',InversionControlParametersEnum(),'format','DoubleMat','mattype',3);
168 WriteData(fid,'data',num_control_parameters,'enum',InversionNumControlParametersEnum(),'format','Integer');
169
170 %process cost functions
171 num_cost_functions=size(obj.cost_functions,2);
172 data=obj.cost_functions;
173 pos=find(obj.cost_functions==101); data(pos)=SurfaceAbsVelMisfitEnum();
174 pos=find(obj.cost_functions==102); data(pos)=SurfaceRelVelMisfitEnum();
175 pos=find(obj.cost_functions==103); data(pos)=SurfaceLogVelMisfitEnum();
176 pos=find(obj.cost_functions==104); data(pos)=SurfaceLogVxVyMisfitEnum();
177 pos=find(obj.cost_functions==105); data(pos)=SurfaceAverageVelMisfitEnum();
178 pos=find(obj.cost_functions==201); data(pos)=ThicknessAbsMisfitEnum();
179 pos=find(obj.cost_functions==501); data(pos)=DragCoefficientAbsGradientEnum();
180 pos=find(obj.cost_functions==502); data(pos)=RheologyBbarAbsGradientEnum();
181 pos=find(obj.cost_functions==503); data(pos)=ThicknessAbsGradientEnum();
182 pos=find(obj.cost_functions==504); data(pos)=ThicknessAlongGradientEnum();
183 pos=find(obj.cost_functions==505); data(pos)=ThicknessAcrossGradientEnum();
184 pos=find(obj.cost_functions==506); data(pos)=BalancethicknessMisfitEnum();
185 WriteData(fid,'data',data,'enum',InversionCostFunctionsEnum(),'format','DoubleMat','mattype',3);
186 WriteData(fid,'data',num_cost_functions,'enum',InversionNumCostFunctionsEnum(),'format','Integer');
187 end % }}}
188 end
189end
Note: See TracBrowser for help on using the repository browser.