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

Last change on this file since 25836 was 25836, checked in by Mathieu Morlighem, 4 years ago

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