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