source: issm/trunk-jpl/src/m/classes/inversion.py@ 16764

Last change on this file since 16764 was 16764, checked in by Eric.Larour, 11 years ago

CHG: checkfield could not run for massatfluggate or misfit classes, because these classes do not
have access to their own fields through the model! So we modified checkfield.m and checkfield.py to accept
either directly a field, or indirectly a fieldname which then is used to retrieve from the model.
File size: 10.0 KB
Line 
1import numpy
2import copy
3from fielddisplay import fielddisplay
4from EnumDefinitions import *
5from StringToEnum import StringToEnum
6from checkfield import *
7from WriteData import *
8
9class inversion(object):
10 """
11 INVERSION class definition
12
13 Usage:
14 inversion=inversion();
15 """
16
17 def __init__(self): # {{{
18 self.iscontrol = 0
19 self.tao = 0
20 self.incomplete_adjoint = 0
21 self.control_parameters = float('NaN')
22 self.nsteps = 0
23 self.maxiter_per_step = float('NaN')
24 self.cost_functions = float('NaN')
25 self.cost_functions_coefficients = float('NaN')
26 self.gradient_scaling = float('NaN')
27 self.cost_function_threshold = 0
28 self.min_parameters = float('NaN')
29 self.max_parameters = float('NaN')
30 self.step_threshold = float('NaN')
31 self.vx_obs = float('NaN')
32 self.vy_obs = float('NaN')
33 self.vz_obs = float('NaN')
34 self.vel_obs = float('NaN')
35 self.thickness_obs = float('NaN')
36
37 #set defaults
38 self.setdefaultparameters()
39
40 #}}}
41 def __repr__(self): # {{{
42 string=' inversion parameters:'
43 string="%s\n%s"%(string,fielddisplay(self,'iscontrol','is inversion activated?'))
44 string="%s\n%s"%(string,fielddisplay(self,'incomplete_adjoint','1: linear viscosity, 0: non-linear viscosity'))
45 string="%s\n%s"%(string,fielddisplay(self,'control_parameters','ex: {''FrictionCoefficient''}, or {''MaterialsRheologyBbar''}'))
46 string="%s\n%s"%(string,fielddisplay(self,'nsteps','number of optimization searches'))
47 string="%s\n%s"%(string,fielddisplay(self,'cost_functions','indicate the type of response for each optimization step'))
48 string="%s\n%s"%(string,fielddisplay(self,'cost_functions_coefficients','cost_functions_coefficients applied to the misfit of each vertex and for each control_parameter'))
49 string="%s\n%s"%(string,fielddisplay(self,'cost_function_threshold','misfit convergence criterion. Default is 1%, NaN if not applied'))
50 string="%s\n%s"%(string,fielddisplay(self,'maxiter_per_step','maximum iterations during each optimization step'))
51 string="%s\n%s"%(string,fielddisplay(self,'gradient_scaling','scaling factor on gradient direction during optimization, for each optimization step'))
52 string="%s\n%s"%(string,fielddisplay(self,'step_threshold','decrease threshold for misfit, default is 30%'))
53 string="%s\n%s"%(string,fielddisplay(self,'min_parameters','absolute minimum acceptable value of the inversed parameter on each vertex'))
54 string="%s\n%s"%(string,fielddisplay(self,'max_parameters','absolute maximum acceptable value of the inversed parameter on each vertex'))
55 string="%s\n%s"%(string,fielddisplay(self,'vx_obs','observed velocity x component [m/yr]'))
56 string="%s\n%s"%(string,fielddisplay(self,'vy_obs','observed velocity y component [m/yr]'))
57 string="%s\n%s"%(string,fielddisplay(self,'vel_obs','observed velocity magnitude [m/yr]'))
58 string="%s\n%s"%(string,fielddisplay(self,'thickness_obs','observed thickness [m]'))
59 string="%s\n%s"%(string,'Available cost functions:')
60 string="%s\n%s"%(string,' 101: SurfaceAbsVelMisfit')
61 string="%s\n%s"%(string,' 102: SurfaceRelVelMisfit')
62 string="%s\n%s"%(string,' 103: SurfaceLogVelMisfit')
63 string="%s\n%s"%(string,' 104: SurfaceLogVxVyMisfit')
64 string="%s\n%s"%(string,' 105: SurfaceAverageVelMisfit')
65 string="%s\n%s"%(string,' 201: ThicknessAbsMisfit')
66 string="%s\n%s"%(string,' 501: DragCoefficientAbsGradient')
67 string="%s\n%s"%(string,' 502: RheologyBbarAbsGradient')
68 string="%s\n%s"%(string,' 503: ThicknessAbsGradient')
69 return string
70 #}}}
71 def setdefaultparameters(self): # {{{
72
73 #default is incomplete adjoint for now
74 self.incomplete_adjoint=1
75
76 #parameter to be inferred by control methods (only
77 #drag and B are supported yet)
78 self.control_parameters='FrictionCoefficient'
79
80 #number of steps in the control methods
81 self.nsteps=20
82
83 #maximum number of iteration in the optimization algorithm for
84 #each step
85 self.maxiter_per_step=20*numpy.ones(self.nsteps)
86
87 #the inversed parameter is updated as follows:
88 #new_par=old_par + gradient_scaling(n)*C*gradient with C in [0 1];
89 #usually the gradient_scaling must be of the order of magnitude of the
90 #inversed parameter (10^8 for B, 50 for drag) and can be decreased
91 #after the first iterations
92 self.gradient_scaling=50*numpy.ones((self.nsteps,1))
93
94 #several responses can be used:
95 self.cost_functions=101
96
97 #step_threshold is used to speed up control method. When
98 #misfit(1)/misfit(0) < self.step_threshold, we go directly to
99 #the next step
100 self.step_threshold=.7*numpy.ones(self.nsteps) #30 per cent decrement
101
102 #cost_function_threshold is a criteria to stop the control methods.
103 #if J[n]-J[n-1]/J[n] < criteria, the control run stops
104 #NaN if not applied
105 self.cost_function_threshold=float('NaN') #not activated
106
107 return self
108 #}}}
109 def checkconsistency(self,md,solution,analyses): # {{{
110
111 #Early return
112 if not self.iscontrol:
113 return md
114
115 num_controls=numpy.size(md.inversion.control_parameters)
116 num_costfunc=numpy.size(md.inversion.cost_functions)
117
118 md = checkfield(md,'fieldname','inversion.iscontrol','values',[0,1])
119 md = checkfield(md,'fieldname','inversion.tao','values',[0,1])
120 md = checkfield(md,'fieldname','inversion.incomplete_adjoint','values',[0,1])
121 md = checkfield(md,'fieldname','inversion.control_parameters','cell',1,'values',['BalancethicknessThickeningRate','FrictionCoefficient','MaterialsRheologyBbar','DamageDbar','Vx','Vy'])
122 md = checkfield(md,'fieldname','inversion.nsteps','numel',[1],'>=',0)
123 md = checkfield(md,'fieldname','inversion.maxiter_per_step','size',[md.inversion.nsteps],'>=',0)
124 md = checkfield(md,'fieldname','inversion.step_threshold','size',[md.inversion.nsteps])
125 md = checkfield(md,'fieldname','inversion.cost_functions','size',[num_costfunc],'values',[101,102,103,104,105,201,501,502,503,504,505])
126 md = checkfield(md,'fieldname','inversion.cost_functions_coefficients','size',[md.mesh.numberofvertices,num_costfunc],'>=',0)
127 md = checkfield(md,'fieldname','inversion.gradient_scaling','size',[md.inversion.nsteps,num_controls])
128 md = checkfield(md,'fieldname','inversion.min_parameters','size',[md.mesh.numberofvertices,num_controls])
129 md = checkfield(md,'fieldname','inversion.max_parameters','size',[md.mesh.numberofvertices,num_controls])
130
131 #Only SSA, HO and FS are supported right now
132 if solution==StressbalanceSolutionEnum():
133 if not (md.flowequation.isSSA or md.flowequation.isHO or md.flowequation.isFS):
134 md.checkmessage("'inversion can only be performed for SSA, HO or FS ice flow models");
135
136 if solution==BalancethicknessSolutionEnum():
137 md = checkfield(md,'fieldname','inversion.thickness_obs','size',[md.mesh.numberofvertices],'NaN',1)
138 else:
139 md = checkfield(md,'fieldname','inversion.vx_obs','size',[md.mesh.numberofvertices],'NaN',1)
140 md = checkfield(md,'fieldname','inversion.vy_obs','size',[md.mesh.numberofvertices],'NaN',1)
141
142 return md
143 # }}}
144 def marshall(self,md,fid): # {{{
145
146 yts=365.0*24.0*3600.0
147
148 WriteData(fid,'object',self,'fieldname','iscontrol','format','Boolean')
149 WriteData(fid,'object',self,'fieldname','tao','format','Boolean')
150 WriteData(fid,'object',self,'fieldname','incomplete_adjoint','format','Boolean')
151 if not self.iscontrol:
152 return
153 WriteData(fid,'object',self,'fieldname','nsteps','format','Integer')
154 WriteData(fid,'object',self,'fieldname','maxiter_per_step','format','DoubleMat','mattype',3)
155 WriteData(fid,'object',self,'fieldname','cost_functions_coefficients','format','DoubleMat','mattype',1)
156 WriteData(fid,'object',self,'fieldname','gradient_scaling','format','DoubleMat','mattype',3)
157 WriteData(fid,'object',self,'fieldname','cost_function_threshold','format','Double')
158 WriteData(fid,'object',self,'fieldname','min_parameters','format','DoubleMat','mattype',3)
159 WriteData(fid,'object',self,'fieldname','max_parameters','format','DoubleMat','mattype',3)
160 WriteData(fid,'object',self,'fieldname','step_threshold','format','DoubleMat','mattype',3)
161 WriteData(fid,'object',self,'fieldname','vx_obs','format','DoubleMat','mattype',1,'scale',1./yts)
162 WriteData(fid,'object',self,'fieldname','vy_obs','format','DoubleMat','mattype',1,'scale',1./yts)
163 WriteData(fid,'object',self,'fieldname','vz_obs','format','DoubleMat','mattype',1,'scale',1./yts)
164 WriteData(fid,'object',self,'fieldname','thickness_obs','format','DoubleMat','mattype',1)
165
166 #process control parameters
167 num_control_parameters=len(self.control_parameters)
168 data=numpy.array([StringToEnum(control_parameter)[0] for control_parameter in self.control_parameters]).reshape(1,-1)
169 WriteData(fid,'data',data,'enum',InversionControlParametersEnum(),'format','DoubleMat','mattype',3)
170 WriteData(fid,'data',num_control_parameters,'enum',InversionNumControlParametersEnum(),'format','Integer')
171
172 #process cost functions
173 num_cost_functions=numpy.size(self.cost_functions)
174 data=copy.deepcopy(self.cost_functions)
175 data=[SurfaceAbsVelMisfitEnum() if x==101 else x for x in data]
176 data=[SurfaceRelVelMisfitEnum() if x==102 else x for x in data]
177 data=[SurfaceLogVelMisfitEnum() if x==103 else x for x in data]
178 data=[SurfaceLogVxVyMisfitEnum() if x==104 else x for x in data]
179 data=[SurfaceAverageVelMisfitEnum() if x==105 else x for x in data]
180 data=[ThicknessAbsMisfitEnum() if x==201 else x for x in data]
181 data=[DragCoefficientAbsGradientEnum() if x==501 else x for x in data]
182 data=[RheologyBbarAbsGradientEnum() if x==502 else x for x in data]
183 data=[ThicknessAbsGradientEnum() if x==503 else x for x in data]
184 data=[ThicknessAlongGradientEnum() if x==504 else x for x in data]
185 data=[ThicknessAcrossGradientEnum() if x==505 else x for x in data]
186 WriteData(fid,'data',numpy.array(data).reshape(1,-1),'enum',InversionCostFunctionsEnum(),'format','DoubleMat','mattype',3)
187 WriteData(fid,'data',num_cost_functions,'enum',InversionNumCostFunctionsEnum(),'format','Integer')
188 # }}}
Note: See TracBrowser for help on using the repository browser.