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

Last change on this file since 21303 was 21303, checked in by bdef, 8 years ago

CHG: major cosmetic clean up we now ask shape (N,) rather than (N,1) also uniformised numpy imports
File size: 10.1 KB
Line 
1import numpy as np
2from project3d import project3d
3from fielddisplay import fielddisplay
4from checkfield import checkfield
5from WriteData import WriteData
6from supportedcontrols import supportedcontrols
7from supportedcostfunctions import supportedcostfunctions
8from marshallcostfunctions import marshallcostfunctions
9
10class inversion(object):
11 """
12 INVERSION class definition
13
14 Usage:
15 inversion=inversion()
16 """
17
18 def __init__(self): # {{{
19 self.iscontrol = 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 = ''
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 self.surface_obs = float('NaN')
37
38 #set defaults
39 self.setdefaultparameters()
40
41 #}}}
42 def __repr__(self): # {{{
43 string=' inversion parameters:'
44 string="%s\n%s"%(string,fielddisplay(self,'iscontrol','is inversion activated?'))
45 string="%s\n%s"%(string,fielddisplay(self,'incomplete_adjoint','1: linear viscosity, 0: non-linear viscosity'))
46 string="%s\n%s"%(string,fielddisplay(self,'control_parameters','ex: {''FrictionCoefficient''}, or {''MaterialsRheologyBbar''}'))
47 string="%s\n%s"%(string,fielddisplay(self,'nsteps','number of optimization searches'))
48 string="%s\n%s"%(string,fielddisplay(self,'cost_functions','indicate the type of response for each optimization step'))
49 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'))
50 string="%s\n%s"%(string,fielddisplay(self,'cost_function_threshold','misfit convergence criterion. Default is 1%, NaN if not applied'))
51 string="%s\n%s"%(string,fielddisplay(self,'maxiter_per_step','maximum iterations during each optimization step'))
52 string="%s\n%s"%(string,fielddisplay(self,'gradient_scaling','scaling factor on gradient direction during optimization, for each optimization step'))
53 string="%s\n%s"%(string,fielddisplay(self,'step_threshold','decrease threshold for misfit, default is 30%'))
54 string="%s\n%s"%(string,fielddisplay(self,'min_parameters','absolute minimum acceptable value of the inversed parameter on each vertex'))
55 string="%s\n%s"%(string,fielddisplay(self,'max_parameters','absolute maximum acceptable value of the inversed parameter on each vertex'))
56 string="%s\n%s"%(string,fielddisplay(self,'vx_obs','observed velocity x component [m/yr]'))
57 string="%s\n%s"%(string,fielddisplay(self,'vy_obs','observed velocity y component [m/yr]'))
58 string="%s\n%s"%(string,fielddisplay(self,'vel_obs','observed velocity magnitude [m/yr]'))
59 string="%s\n%s"%(string,fielddisplay(self,'thickness_obs','observed thickness [m]'))
60 string="%s\n%s"%(string,fielddisplay(self,'surface_obs','observed surface elevation [m]'))
61 string="%s\n%s"%(string,'Available cost functions:')
62 string="%s\n%s"%(string,' 101: SurfaceAbsVelMisfit')
63 string="%s\n%s"%(string,' 102: SurfaceRelVelMisfit')
64 string="%s\n%s"%(string,' 103: SurfaceLogVelMisfit')
65 string="%s\n%s"%(string,' 104: SurfaceLogVxVyMisfit')
66 string="%s\n%s"%(string,' 105: SurfaceAverageVelMisfit')
67 string="%s\n%s"%(string,' 201: ThicknessAbsMisfit')
68 string="%s\n%s"%(string,' 501: DragCoefficientAbsGradient')
69 string="%s\n%s"%(string,' 502: RheologyBbarAbsGradient')
70 string="%s\n%s"%(string,' 503: ThicknessAbsGradient')
71 return string
72 #}}}
73 def extrude(self,md): # {{{
74 self.vx_obs=project3d(md,'vector',self.vx_obs,'type','node')
75 self.vy_obs=project3d(md,'vector',self.vy_obs,'type','node')
76 self.vel_obs=project3d(md,'vector',self.vel_obs,'type','node')
77 self.thickness_obs=project3d(md,'vector',self.thickness_obs,'type','node')
78 if not np.any(np.isnan(self.cost_functions_coefficients)):
79 self.cost_functions_coefficients=project3d(md,'vector',self.cost_functions_coefficients,'type','node')
80 if not np.any(np.isnan(self.min_parameters)):
81 self.min_parameters=project3d(md,'vector',self.min_parameters,'type','node')
82 if not np.any(np.isnan(self.max_parameters)):
83 self.max_parameters=project3d(md,'vector',self.max_parameters,'type','node')
84 return self
85 #}}}
86 def setdefaultparameters(self): # {{{
87
88 #default is incomplete adjoint for now
89 self.incomplete_adjoint=1
90
91 #parameter to be inferred by control methods (only
92 #drag and B are supported yet)
93 self.control_parameters='FrictionCoefficient'
94
95 #number of steps in the control methods
96 self.nsteps=20
97
98 #maximum number of iteration in the optimization algorithm for
99 #each step
100 self.maxiter_per_step=20*np.ones(self.nsteps)
101
102 #the inversed parameter is updated as follows:
103 #new_par=old_par + gradient_scaling(n)*C*gradient with C in [0 1];
104 #usually the gradient_scaling must be of the order of magnitude of the
105 #inversed parameter (10^8 for B, 50 for drag) and can be decreased
106 #after the first iterations
107 self.gradient_scaling=50*np.ones((self.nsteps,1))
108
109 #several responses can be used:
110 self.cost_functions=[101,]
111
112 #step_threshold is used to speed up control method. When
113 #misfit(1)/misfit(0) < self.step_threshold, we go directly to
114 #the next step
115 self.step_threshold=.7*np.ones(self.nsteps) #30 per cent decrement
116
117 #cost_function_threshold is a criteria to stop the control methods.
118 #if J[n]-J[n-1]/J[n] < criteria, the control run stops
119 #NaN if not applied
120 self.cost_function_threshold=float('NaN') #not activated
121
122 return self
123 #}}}
124 def checkconsistency(self,md,solution,analyses): # {{{
125
126 #Early return
127 if not self.iscontrol:
128 return md
129
130 num_controls=np.size(md.inversion.control_parameters)
131 num_costfunc=np.size(md.inversion.cost_functions)
132
133 md = checkfield(md,'fieldname','inversion.iscontrol','values',[0,1])
134 md = checkfield(md,'fieldname','inversion.incomplete_adjoint','values',[0,1])
135 md = checkfield(md,'fieldname','inversion.control_parameters','cell',1,'values',supportedcontrols())
136 md = checkfield(md,'fieldname','inversion.nsteps','numel',[1],'>=',0)
137 md = checkfield(md,'fieldname','inversion.maxiter_per_step','size',[md.inversion.nsteps],'>=',0)
138 md = checkfield(md,'fieldname','inversion.step_threshold','size',[md.inversion.nsteps])
139 md = checkfield(md,'fieldname','inversion.cost_functions','size',[num_costfunc],'values',supportedcostfunctions())
140 md = checkfield(md,'fieldname','inversion.cost_functions_coefficients','size',[md.mesh.numberofvertices,num_costfunc],'>=',0)
141 md = checkfield(md,'fieldname','inversion.gradient_scaling','size',[md.inversion.nsteps,num_controls])
142 md = checkfield(md,'fieldname','inversion.min_parameters','size',[md.mesh.numberofvertices,num_controls])
143 md = checkfield(md,'fieldname','inversion.max_parameters','size',[md.mesh.numberofvertices,num_controls])
144
145 #Only SSA, HO and FS are supported right now
146 if solution=='StressbalanceSolution':
147 if not (md.flowequation.isSSA or md.flowequation.isHO or md.flowequation.isFS or md.flowequation.isL1L2):
148 md.checkmessage("'inversion can only be performed for SSA, HO or FS ice flow models");
149
150 if solution=='BalancethicknessSolution':
151 md = checkfield(md,'fieldname','inversion.thickness_obs','size',[md.mesh.numberofvertices],'NaN',1,'Inf',1)
152 else:
153 md = checkfield(md,'fieldname','inversion.vx_obs','size',[md.mesh.numberofvertices],'NaN',1,'Inf',1)
154 md = checkfield(md,'fieldname','inversion.vy_obs','size',[md.mesh.numberofvertices],'NaN',1,'Inf',1)
155
156 return md
157 # }}}
158 def marshall(self,prefix,md,fid): # {{{
159
160 yts=md.constants.yts
161
162 WriteData(fid,prefix,'name','md.inversion.type','data',0,'format','Integer')
163 WriteData(fid,prefix,'object',self,'fieldname','iscontrol','format','Boolean')
164 WriteData(fid,prefix,'object',self,'fieldname','incomplete_adjoint','format','Boolean')
165 if not self.iscontrol:
166 return
167 WriteData(fid,prefix,'object',self,'fieldname','nsteps','format','Integer')
168 WriteData(fid,prefix,'object',self,'fieldname','maxiter_per_step','format','DoubleMat','mattype',3)
169 WriteData(fid,prefix,'object',self,'fieldname','cost_functions_coefficients','format','DoubleMat','mattype',1)
170 WriteData(fid,prefix,'object',self,'fieldname','gradient_scaling','format','DoubleMat','mattype',3)
171 WriteData(fid,prefix,'object',self,'fieldname','cost_function_threshold','format','Double')
172 WriteData(fid,prefix,'object',self,'fieldname','min_parameters','format','DoubleMat','mattype',3)
173 WriteData(fid,prefix,'object',self,'fieldname','max_parameters','format','DoubleMat','mattype',3)
174 WriteData(fid,prefix,'object',self,'fieldname','step_threshold','format','DoubleMat','mattype',3)
175 WriteData(fid,prefix,'object',self,'fieldname','vx_obs','format','DoubleMat','mattype',1,'scale',1./yts)
176 WriteData(fid,prefix,'object',self,'fieldname','vy_obs','format','DoubleMat','mattype',1,'scale',1./yts)
177 WriteData(fid,prefix,'object',self,'fieldname','vz_obs','format','DoubleMat','mattype',1,'scale',1./yts)
178 WriteData(fid,prefix,'object',self,'fieldname','thickness_obs','format','DoubleMat','mattype',1)
179 WriteData(fid,prefix,'object',self,'fieldname','surface_obs','format','DoubleMat','mattype',1)
180
181 #process control parameters
182 num_control_parameters=len(self.control_parameters)
183 WriteData(fid,prefix,'object',self,'fieldname','control_parameters','format','StringArray')
184 WriteData(fid,prefix,'data',num_control_parameters,'name','md.inversion.num_control_parameters','format','Integer')
185
186 #process cost functions
187 num_cost_functions=np.size(self.cost_functions)
188 data=marshallcostfunctions(self.cost_functions)
189 WriteData(fid,prefix,'data',data,'name','md.inversion.cost_functions','format','StringArray')
190 WriteData(fid,prefix,'data',num_cost_functions,'name','md.inversion.num_cost_functions','format','Integer')
191 # }}}
Note: See TracBrowser for help on using the repository browser.