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

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

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