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

Last change on this file since 13093 was 13093, checked in by jschierm, 13 years ago

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