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

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

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