source: issm/trunk/src/m/classes/m1qn3inversion.py@ 21341

Last change on this file since 21341 was 21341, checked in by Mathieu Morlighem, 8 years ago

merged trunk-jpl and trunk for revision 21337
File size: 10.2 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 m1qn3inversion(object):
11 '''
12 M1QN3 class definition
13
14 Usage:
15 m1qn3inversion=m1qn3inversion()
16 '''
17
18 def __init__(self,*args): # {{{
19
20 if not len(args):
21 print 'empty init'
22 self.iscontrol = 0
23 self.incomplete_adjoint = 0
24 self.control_parameters = float('NaN')
25 self.control_scaling_factors = float('NaN')
26 self.maxsteps = 0
27 self.maxiter = 0
28 self.dxmin = 0.
29 self.gttol = 0.
30 self.cost_functions = float('NaN')
31 self.cost_functions_coefficients = float('NaN')
32 self.min_parameters = float('NaN')
33 self.max_parameters = float('NaN')
34 self.vx_obs = float('NaN')
35 self.vy_obs = float('NaN')
36 self.vz_obs = float('NaN')
37 self.vel_obs = float('NaN')
38 self.thickness_obs = float('NaN')
39
40 #set defaults
41 self.setdefaultparameters()
42 elif len(args)==1 and args[0].__module__=='inversion':
43 print 'converting inversion to m1qn3inversion'
44 inv=args[0]
45 #first call setdefaultparameters:
46 self.setdefaultparameters()
47
48 #then go fish whatever is available in the inversion object provided to the constructor
49 self.iscontrol = inv.iscontrol
50 self.incomplete_adjoint = inv.incomplete_adjoint
51 self.control_parameters = inv.control_parameters
52 self.maxsteps = inv.nsteps
53 self.cost_functions = inv.cost_functions
54 self.cost_functions_coefficients = inv.cost_functions_coefficients
55 self.min_parameters = inv.min_parameters
56 self.max_parameters = inv.max_parameters
57 self.vx_obs = inv.vx_obs
58 self.vy_obs = inv.vy_obs
59 self.vz_obs = inv.vz_obs
60 self.vel_obs = inv.vel_obs
61 self.thickness_obs = inv.thickness_obs
62 else:
63 raise Exception('constructor not supported')
64 #}}}
65 def __repr__(self): # {{{
66 string=' m1qn3inversion parameters:'
67 string="%s\n%s"%(string,fielddisplay(self,'iscontrol','is inversion activated?'))
68 string="%s\n%s"%(string,fielddisplay(self,'incomplete_adjoint','1: linear viscosity, 0: non-linear viscosity'))
69 string="%s\n%s"%(string,fielddisplay(self,'control_parameters','ex: [''FrictionCoefficient''], or [''MaterialsRheologyBbar'']'))
70 string="%s\n%s"%(string,fielddisplay(self,'control_scaling_factors','order of magnitude of each control (useful for multi-parameter optimization)'))
71 string="%s\n%s"%(string,fielddisplay(self,'maxsteps','maximum number of iterations (gradient computation)'))
72 string="%s\n%s"%(string,fielddisplay(self,'maxiter','maximum number of Function evaluation (forward run)'))
73 string="%s\n%s"%(string,fielddisplay(self,'dxmin','convergence criterion: two points less than dxmin from eachother (sup-norm) are considered identical'))
74 string="%s\n%s"%(string,fielddisplay(self,'gttol','||g(X)||/||g(X0)|| (g(X0): gradient at initial guess X0)'))
75 string="%s\n%s"%(string,fielddisplay(self,'cost_functions','indicate the type of response for each optimization step'))
76 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'))
77 string="%s\n%s"%(string,fielddisplay(self,'min_parameters','absolute minimum acceptable value of the inversed parameter on each vertex'))
78 string="%s\n%s"%(string,fielddisplay(self,'max_parameters','absolute maximum acceptable value of the inversed parameter on each vertex'))
79 string="%s\n%s"%(string,fielddisplay(self,'vx_obs','observed velocity x component [m/yr]'))
80 string="%s\n%s"%(string,fielddisplay(self,'vy_obs','observed velocity y component [m/yr]'))
81 string="%s\n%s"%(string,fielddisplay(self,'vel_obs','observed velocity magnitude [m/yr]'))
82 string="%s\n%s"%(string,fielddisplay(self,'thickness_obs','observed thickness [m]'))
83 string="%s\n%s"%(string,'Available cost functions:')
84 string="%s\n%s"%(string,' 101: SurfaceAbsVelMisfit')
85 string="%s\n%s"%(string,' 102: SurfaceRelVelMisfit')
86 string="%s\n%s"%(string,' 103: SurfaceLogVelMisfit')
87 string="%s\n%s"%(string,' 104: SurfaceLogVxVyMisfit')
88 string="%s\n%s"%(string,' 105: SurfaceAverageVelMisfit')
89 string="%s\n%s"%(string,' 201: ThicknessAbsMisfit')
90 string="%s\n%s"%(string,' 501: DragCoefficientAbsGradient')
91 string="%s\n%s"%(string,' 502: RheologyBbarAbsGradient')
92 string="%s\n%s"%(string,' 503: ThicknessAbsGradient')
93 return string
94 #}}}
95 def extrude(self,md): # {{{
96 self.vx_obs=project3d(md,'vector',self.vx_obs,'type','node')
97 self.vy_obs=project3d(md,'vector',self.vy_obs,'type','node')
98 self.vel_obs=project3d(md,'vector',self.vel_obs,'type','node')
99 self.thickness_obs=project3d(md,'vector',self.thickness_obs,'type','node')
100 if not np.any(np.isnan(self.cost_functions_coefficients)):
101 self.cost_functions_coefficients=project3d(md,'vector',self.cost_functions_coefficients,'type','node')
102 if not np.any(np.isnan(self.min_parameters)):
103 self.min_parameters=project3d(md,'vector',self.min_parameters,'type','node')
104 if not np.any(np.isnan(self.max_parameters)):
105 self.max_parameters=project3d(md,'vector',self.max_parameters,'type','node')
106 return self
107 #}}}
108 def setdefaultparameters(self): # {{{
109
110 #default is incomplete adjoint for now
111 self.incomplete_adjoint=1
112
113 #parameter to be inferred by control methods (only
114 #drag and B are supported yet)
115 self.control_parameters='FrictionCoefficient'
116
117 #Scaling factor for each control
118 self.control_scaling_factors=1
119
120 #number of iterations
121 self.maxsteps=20
122 self.maxiter=40
123
124 #several responses can be used:
125 self.cost_functions=101
126
127 #m1qn3 parameters
128 self.dxmin = 0.1
129 self.gttol = 1e-4
130
131 return self
132 #}}}
133 def checkconsistency(self,md,solution,analyses): # {{{
134
135 #Early return
136 if not self.iscontrol:
137 return md
138
139 num_controls=np.size(md.inversion.control_parameters)
140 num_costfunc=np.size(md.inversion.cost_functions)
141
142 md = checkfield(md,'fieldname','inversion.iscontrol','values',[0,1])
143 md = checkfield(md,'fieldname','inversion.incomplete_adjoint','values',[0,1])
144 md = checkfield(md,'fieldname','inversion.control_parameters','cell',1,'values',supportedcontrols())
145 md = checkfield(md,'fieldname','inversion.control_scaling_factors','size',[num_controls],'>',0,'NaN',1,'Inf',1)
146 md = checkfield(md,'fieldname','inversion.maxsteps','numel',[1],'>=',0)
147 md = checkfield(md,'fieldname','inversion.maxiter','numel',[1],'>=',0)
148 md = checkfield(md,'fieldname','inversion.dxmin','numel',[1],'>',0.)
149 md = checkfield(md,'fieldname','inversion.gttol','numel',[1],'>',0.)
150 md = checkfield(md,'fieldname','inversion.cost_functions','size',[num_costfunc],'values',supportedcostfunctions())
151 md = checkfield(md,'fieldname','inversion.cost_functions_coefficients','size',[md.mesh.numberofvertices,num_costfunc],'>=',0)
152 md = checkfield(md,'fieldname','inversion.min_parameters','size',[md.mesh.numberofvertices,num_controls])
153 md = checkfield(md,'fieldname','inversion.max_parameters','size',[md.mesh.numberofvertices,num_controls])
154
155 if solution=='BalancethicknessSolution':
156 md = checkfield(md,'fieldname','inversion.thickness_obs','size',[md.mesh.numberofvertices],'NaN',1,'Inf',1)
157 else:
158 md = checkfield(md,'fieldname','inversion.vx_obs','size',[md.mesh.numberofvertices],'NaN',1,'Inf',1)
159 md = checkfield(md,'fieldname','inversion.vy_obs','size',[md.mesh.numberofvertices],'NaN',1,'Inf',1)
160
161 return md
162 # }}}
163 def marshall(self,prefix,md,fid): # {{{
164
165 yts=md.constants.yts
166
167 WriteData(fid,prefix,'object',self,'class','inversion','fieldname','iscontrol','format','Boolean')
168 WriteData(fid,prefix,'name','md.inversion.type','data',2,'format','Integer')
169 if not self.iscontrol:
170 return
171 WriteData(fid,prefix,'object',self,'class','inversion','fieldname','incomplete_adjoint','format','Boolean')
172 WriteData(fid,prefix,'object',self,'class','inversion','fieldname','control_scaling_factors','format','DoubleMat','mattype',3)
173 WriteData(fid,prefix,'object',self,'class','inversion','fieldname','maxsteps','format','Integer')
174 WriteData(fid,prefix,'object',self,'class','inversion','fieldname','maxiter','format','Integer')
175 WriteData(fid,prefix,'object',self,'class','inversion','fieldname','dxmin','format','Double')
176 WriteData(fid,prefix,'object',self,'class','inversion','fieldname','gttol','format','Double')
177 WriteData(fid,prefix,'object',self,'class','inversion','fieldname','cost_functions_coefficients','format','DoubleMat','mattype',1)
178 WriteData(fid,prefix,'object',self,'class','inversion','fieldname','min_parameters','format','DoubleMat','mattype',3)
179 WriteData(fid,prefix,'object',self,'class','inversion','fieldname','max_parameters','format','DoubleMat','mattype',3)
180 WriteData(fid,prefix,'object',self,'class','inversion','fieldname','vx_obs','format','DoubleMat','mattype',1,'scale',1./yts)
181 WriteData(fid,prefix,'object',self,'class','inversion','fieldname','vy_obs','format','DoubleMat','mattype',1,'scale',1./yts)
182 WriteData(fid,prefix,'object',self,'class','inversion','fieldname','vz_obs','format','DoubleMat','mattype',1,'scale',1./yts)
183 WriteData(fid,prefix,'object',self,'class','inversion','fieldname','thickness_obs','format','DoubleMat','mattype',1)
184
185 #process control parameters
186 num_control_parameters=len(self.control_parameters)
187 WriteData(fid,prefix,'object',self,'fieldname','control_parameters','format','StringArray')
188 WriteData(fid,prefix,'data',num_control_parameters,'name','md.inversion.num_control_parameters','format','Integer')
189
190 #process cost functions
191 num_cost_functions=np.size(self.cost_functions)
192 data=marshallcostfunctions(self.cost_functions)
193 WriteData(fid,prefix,'data',data,'name','md.inversion.cost_functions','format','StringArray')
194 WriteData(fid,prefix,'data',num_cost_functions,'name','md.inversion.num_cost_functions','format','Integer')
195 # }}}
Note: See TracBrowser for help on using the repository browser.