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

Last change on this file since 19105 was 19105, checked in by Mathieu Morlighem, 10 years ago

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