| 1 | import numpy as np
|
|---|
| 2 |
|
|---|
| 3 | from checkfield import checkfield
|
|---|
| 4 | from fielddisplay import fielddisplay
|
|---|
| 5 | from marshallcostfunctions import marshallcostfunctions
|
|---|
| 6 | from project3d import project3d
|
|---|
| 7 | from supportedcontrols import supportedcontrols
|
|---|
| 8 | from supportedcostfunctions import supportedcostfunctions
|
|---|
| 9 | from WriteData import WriteData
|
|---|
| 10 |
|
|---|
| 11 |
|
|---|
| 12 | class inversion(object):
|
|---|
| 13 | """INVERSION class definition
|
|---|
| 14 |
|
|---|
| 15 | Usage:
|
|---|
| 16 | inversion = inversion()
|
|---|
| 17 | """
|
|---|
| 18 |
|
|---|
| 19 | def __init__(self): # {{{
|
|---|
| 20 | self.iscontrol = 0
|
|---|
| 21 | self.incomplete_adjoint = 0
|
|---|
| 22 | self.control_parameters = np.nan
|
|---|
| 23 | self.nsteps = 0
|
|---|
| 24 | self.maxiter_per_step = np.nan
|
|---|
| 25 | self.cost_functions = ''
|
|---|
| 26 | self.cost_functions_coefficients = np.nan
|
|---|
| 27 | self.gradient_scaling = np.nan
|
|---|
| 28 | self.cost_function_threshold = 0
|
|---|
| 29 | self.min_parameters = np.nan
|
|---|
| 30 | self.max_parameters = np.nan
|
|---|
| 31 | self.step_threshold = np.nan
|
|---|
| 32 | self.vx_obs = np.nan
|
|---|
| 33 | self.vy_obs = np.nan
|
|---|
| 34 | self.vz_obs = np.nan
|
|---|
| 35 | self.vel_obs = np.nan
|
|---|
| 36 | self.thickness_obs = np.nan
|
|---|
| 37 | self.surface_obs = np.nan
|
|---|
| 38 |
|
|---|
| 39 | self.setdefaultparameters()
|
|---|
| 40 | # }}}
|
|---|
| 41 |
|
|---|
| 42 | def __repr__(self): # {{{
|
|---|
| 43 | s = ' inversion parameters:\n'
|
|---|
| 44 | s += '{}\n'.format(fielddisplay(self, 'iscontrol', 'is inversion activated?'))
|
|---|
| 45 | s += '{}\n'.format(fielddisplay(self, 'incomplete_adjoint', '1: linear viscosity, 0: non - linear viscosity'))
|
|---|
| 46 | s += '{}\n'.format(fielddisplay(self, 'control_parameters', 'ex: {''FrictionCoefficient''}, or {''MaterialsRheologyBbar''}'))
|
|---|
| 47 | s += '{}\n'.format(fielddisplay(self, 'nsteps', 'number of optimization searches'))
|
|---|
| 48 | s += '{}\n'.format(fielddisplay(self, 'cost_functions', 'indicate the type of response for each optimization step'))
|
|---|
| 49 | s += '{}\n'.format(fielddisplay(self, 'cost_functions_coefficients', 'cost_functions_coefficients applied to the misfit of each vertex and for each control_parameter'))
|
|---|
| 50 | s += '{}\n'.format(fielddisplay(self, 'cost_function_threshold', 'misfit convergence criterion. Default is 1%, NaN if not applied'))
|
|---|
| 51 | s += '{}\n'.format(fielddisplay(self, 'maxiter_per_step', 'maximum iterations during each optimization step'))
|
|---|
| 52 | s += '{}\n'.format(fielddisplay(self, 'gradient_scaling', 'scaling factor on gradient direction during optimization, for each optimization step'))
|
|---|
| 53 | s += '{}\n'.format(fielddisplay(self, 'step_threshold', 'decrease threshold for misfit, default is 30%'))
|
|---|
| 54 | s += '{}\n'.format(fielddisplay(self, 'min_parameters', 'absolute minimum acceptable value of the inversed parameter on each vertex'))
|
|---|
| 55 | s += '{}\n'.format(fielddisplay(self, 'max_parameters', 'absolute maximum acceptable value of the inversed parameter on each vertex'))
|
|---|
| 56 | s += '{}\n'.format(fielddisplay(self, 'vx_obs', 'observed velocity x component [m/yr]'))
|
|---|
| 57 | s += '{}\n'.format(fielddisplay(self, 'vy_obs', 'observed velocity y component [m/yr]'))
|
|---|
| 58 | s += '{}\n'.format(fielddisplay(self, 'vel_obs', 'observed velocity magnitude [m/yr]'))
|
|---|
| 59 | s += '{}\n'.format(fielddisplay(self, 'thickness_obs', 'observed thickness [m]'))
|
|---|
| 60 | s += '{}\n'.format(fielddisplay(self, 'surface_obs', 'observed surface elevation [m]'))
|
|---|
| 61 | s += '{}\n'.format('Available cost functions:')
|
|---|
| 62 | s += '{}\n'.format(' 101: SurfaceAbsVelMisfit')
|
|---|
| 63 | s += '{}\n'.format(' 102: SurfaceRelVelMisfit')
|
|---|
| 64 | s += '{}\n'.format(' 103: SurfaceLogVelMisfit')
|
|---|
| 65 | s += '{}\n'.format(' 104: SurfaceLogVxVyMisfit')
|
|---|
| 66 | s += '{}\n'.format(' 105: SurfaceAverageVelMisfit')
|
|---|
| 67 | s += '{}\n'.format(' 201: ThicknessAbsMisfit')
|
|---|
| 68 | s += '{}\n'.format(' 501: DragCoefficientAbsGradient')
|
|---|
| 69 | s += '{}\n'.format(' 502: RheologyBbarAbsGradient')
|
|---|
| 70 | s += '{}\n'.format(' 503: ThicknessAbsGradient')
|
|---|
| 71 | return s
|
|---|
| 72 | # }}}
|
|---|
| 73 |
|
|---|
| 74 | def setdefaultparameters(self): # {{{
|
|---|
| 75 | #default is incomplete adjoint for now
|
|---|
| 76 | self.incomplete_adjoint = 1
|
|---|
| 77 | #parameter to be inferred by control methods (only
|
|---|
| 78 | #drag and B are supported yet)
|
|---|
| 79 | self.control_parameters = 'FrictionCoefficient'
|
|---|
| 80 | #number of steps in the control methods
|
|---|
| 81 | self.nsteps = 20
|
|---|
| 82 | #maximum number of iteration in the optimization algorithm for
|
|---|
| 83 | #each step
|
|---|
| 84 | self.maxiter_per_step = 20 * np.ones(self.nsteps)
|
|---|
| 85 | #the inversed parameter is updated as follows:
|
|---|
| 86 | #new_par = old_par + gradient_scaling(n) * C * gradient with C in [0 1]
|
|---|
| 87 | #usually the gradient_scaling must be of the order of magnitude of the
|
|---|
| 88 | #inversed parameter (1.0e8 for B, 50 for drag) and can be decreased
|
|---|
| 89 | #after the first iterations
|
|---|
| 90 | self.gradient_scaling = 50 * np.ones((self.nsteps, 1))
|
|---|
| 91 | #several responses can be used:
|
|---|
| 92 | self.cost_functions = [101, ]
|
|---|
| 93 | #step_threshold is used to speed up control method. When
|
|---|
| 94 | #misfit(1) / misfit(0) < self.step_threshold, we go directly to
|
|---|
| 95 | #the next step
|
|---|
| 96 | self.step_threshold = 0.7 * np.ones(self.nsteps) #30 per cent decrement
|
|---|
| 97 | #cost_function_threshold is a criteria to stop the control methods.
|
|---|
| 98 | #if J[n] - J[n - 1] / J[n] < criteria, the control run stops
|
|---|
| 99 | #NaN if not applied
|
|---|
| 100 | self.cost_function_threshold = np.nan #not activated
|
|---|
| 101 | return self
|
|---|
| 102 | # }}}
|
|---|
| 103 |
|
|---|
| 104 | def extrude(self, md): # {{{
|
|---|
| 105 | self.vx_obs = project3d(md, 'vector', self.vx_obs, 'type', 'node')
|
|---|
| 106 | self.vy_obs = project3d(md, 'vector', self.vy_obs, 'type', 'node')
|
|---|
| 107 | self.vel_obs = project3d(md, 'vector', self.vel_obs, 'type', 'node')
|
|---|
| 108 | self.thickness_obs = project3d(md, 'vector', self.thickness_obs, 'type', 'node')
|
|---|
| 109 | if not np.any(np.isnan(self.cost_functions_coefficients)):
|
|---|
| 110 | self.cost_functions_coefficients = project3d(md, 'vector', self.cost_functions_coefficients, 'type', 'node')
|
|---|
| 111 | if not np.any(np.isnan(self.min_parameters)):
|
|---|
| 112 | self.min_parameters = project3d(md, 'vector', self.min_parameters, 'type', 'node')
|
|---|
| 113 | if not np.any(np.isnan(self.max_parameters)):
|
|---|
| 114 | self.max_parameters = project3d(md, 'vector', self.max_parameters, 'type', 'node')
|
|---|
| 115 | return self
|
|---|
| 116 | # }}}
|
|---|
| 117 |
|
|---|
| 118 | def checkconsistency(self, md, solution, analyses): # {{{
|
|---|
| 119 | # Early return
|
|---|
| 120 | if not self.iscontrol:
|
|---|
| 121 | return md
|
|---|
| 122 |
|
|---|
| 123 | num_controls = np.size(md.inversion.control_parameters)
|
|---|
| 124 | num_costfunc = np.size(md.inversion.cost_functions)
|
|---|
| 125 |
|
|---|
| 126 | md = checkfield(md, 'fieldname', 'inversion.iscontrol', 'values', [0, 1])
|
|---|
| 127 | md = checkfield(md, 'fieldname', 'inversion.incomplete_adjoint', 'values', [0, 1])
|
|---|
| 128 | md = checkfield(md, 'fieldname', 'inversion.control_parameters', 'cell', 1, 'values', supportedcontrols())
|
|---|
| 129 | md = checkfield(md, 'fieldname', 'inversion.nsteps', 'numel', [1], '>=', 0)
|
|---|
| 130 | md = checkfield(md, 'fieldname', 'inversion.maxiter_per_step', 'size', [md.inversion.nsteps], '>=', 0)
|
|---|
| 131 | md = checkfield(md, 'fieldname', 'inversion.step_threshold', 'size', [md.inversion.nsteps])
|
|---|
| 132 | md = checkfield(md, 'fieldname', 'inversion.cost_functions', 'size', [num_costfunc], 'values', supportedcostfunctions())
|
|---|
| 133 | if num_costfunc == 1:
|
|---|
| 134 | md.inversion.cost_functions_coefficients = np.squeeze(md.inversion.cost_functions_coefficients)
|
|---|
| 135 | md = checkfield(md, 'fieldname', 'inversion.cost_functions_coefficients', 'size', [md.mesh.numberofvertices], '>=', 0)
|
|---|
| 136 | else:
|
|---|
| 137 | md = checkfield(md, 'fieldname', 'inversion.cost_functions_coefficients', 'size', [md.mesh.numberofvertices, num_costfunc], '>=', 0)
|
|---|
| 138 |
|
|---|
| 139 | if num_controls == 1:
|
|---|
| 140 | md.inversion.gradient_scaling = np.squeeze(md.inversion.gradient_scaling)
|
|---|
| 141 | md.inversion.min_parameters = np.squeeze(md.inversion.min_parameters)
|
|---|
| 142 | md.inversion.max_parameters = np.squeeze(md.inversion.max_parameters)
|
|---|
| 143 | md = checkfield(md, 'fieldname', 'inversion.gradient_scaling', 'size', [md.inversion.nsteps])
|
|---|
| 144 | md = checkfield(md, 'fieldname', 'inversion.min_parameters', 'size', [md.mesh.numberofvertices])
|
|---|
| 145 | md = checkfield(md, 'fieldname', 'inversion.max_parameters', 'size', [md.mesh.numberofvertices])
|
|---|
| 146 | else:
|
|---|
| 147 | md = checkfield(md, 'fieldname', 'inversion.gradient_scaling', 'size', [md.inversion.nsteps, num_controls])
|
|---|
| 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 | # Only SSA, HO and FS are supported right now
|
|---|
| 152 | if solution == 'StressbalanceSolution':
|
|---|
| 153 | if not (md.flowequation.isSSA or md.flowequation.isMOLHO or md.flowequation.isHO or md.flowequation.isFS or md.flowequation.isL1L2):
|
|---|
| 154 | md.checkmessage("'inversion can only be performed for SSA, MOLHO, HO or FS ice flow models")
|
|---|
| 155 | if solution == 'BalancethicknessSolution':
|
|---|
| 156 | md = checkfield(md, 'fieldname', 'inversion.thickness_obs', 'size', [md.mesh.numberofvertices], 'NaN', 1, 'Inf', 1)
|
|---|
| 157 | elif solution == 'BalancethicknessSoftSolution':
|
|---|
| 158 | md = checkfield(md, 'fieldname', 'inversion.thickness_obs', 'size', [md.mesh.numberofvertices], 'NaN', 1, 'Inf', 1)
|
|---|
| 159 | else:
|
|---|
| 160 | md = checkfield(md, 'fieldname', 'inversion.vx_obs', 'size', [md.mesh.numberofvertices], 'NaN', 1, 'Inf', 1)
|
|---|
| 161 | md = checkfield(md, 'fieldname', 'inversion.vy_obs', 'size', [md.mesh.numberofvertices], 'NaN', 1, 'Inf', 1)
|
|---|
| 162 | return md
|
|---|
| 163 | # }}}
|
|---|
| 164 |
|
|---|
| 165 | def marshall(self, prefix, md, fid): # {{{
|
|---|
| 166 | yts = md.constants.yts
|
|---|
| 167 |
|
|---|
| 168 | WriteData(fid, prefix, 'name', 'md.inversion.type', 'data', 0, 'format', 'Integer')
|
|---|
| 169 | WriteData(fid, prefix, 'object', self, 'fieldname', 'iscontrol', 'format', 'Boolean')
|
|---|
| 170 | WriteData(fid, prefix, 'object', self, 'fieldname', 'incomplete_adjoint', 'format', 'Boolean')
|
|---|
| 171 | WriteData(fid, prefix, 'object', self, 'fieldname', 'vel_obs', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / yts)
|
|---|
| 172 | if not self.iscontrol:
|
|---|
| 173 | return
|
|---|
| 174 | WriteData(fid, prefix, 'object', self, 'fieldname', 'nsteps', 'format', 'Integer')
|
|---|
| 175 | WriteData(fid, prefix, 'object', self, 'fieldname', 'maxiter_per_step', 'format', 'IntMat', 'mattype', 3)
|
|---|
| 176 | WriteData(fid, prefix, 'object', self, 'fieldname', 'cost_functions_coefficients', 'format', 'DoubleMat', 'mattype', 1)
|
|---|
| 177 | WriteData(fid, prefix, 'object', self, 'fieldname', 'gradient_scaling', 'format', 'DoubleMat', 'mattype', 3)
|
|---|
| 178 | WriteData(fid, prefix, 'object', self, 'fieldname', 'cost_function_threshold', 'format', 'Double')
|
|---|
| 179 | WriteData(fid, prefix, 'object', self, 'fieldname', 'min_parameters', 'format', 'DoubleMat', 'mattype', 3)
|
|---|
| 180 | WriteData(fid, prefix, 'object', self, 'fieldname', 'max_parameters', 'format', 'DoubleMat', 'mattype', 3)
|
|---|
| 181 | WriteData(fid, prefix, 'object', self, 'fieldname', 'step_threshold', 'format', 'DoubleMat', 'mattype', 3)
|
|---|
| 182 | WriteData(fid, prefix, 'object', self, 'fieldname', 'vx_obs', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / yts)
|
|---|
| 183 | WriteData(fid, prefix, 'object', self, 'fieldname', 'vy_obs', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / yts)
|
|---|
| 184 | WriteData(fid, prefix, 'object', self, 'fieldname', 'vz_obs', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / yts)
|
|---|
| 185 | WriteData(fid, prefix, 'object', self, 'fieldname', 'thickness_obs', 'format', 'DoubleMat', 'mattype', 1)
|
|---|
| 186 | WriteData(fid, prefix, 'object', self, 'fieldname', 'surface_obs', 'format', 'DoubleMat', 'mattype', 1)
|
|---|
| 187 |
|
|---|
| 188 | # Process control parameters
|
|---|
| 189 | num_control_parameters = len(self.control_parameters)
|
|---|
| 190 | WriteData(fid, prefix, 'object', self, 'fieldname', 'control_parameters', 'format', 'StringArray')
|
|---|
| 191 | WriteData(fid, prefix, 'data', num_control_parameters, 'name', 'md.inversion.num_control_parameters', 'format', 'Integer')
|
|---|
| 192 |
|
|---|
| 193 | # Process cost functions
|
|---|
| 194 | num_cost_functions = np.size(self.cost_functions)
|
|---|
| 195 | data = marshallcostfunctions(self.cost_functions)
|
|---|
| 196 | WriteData(fid, prefix, 'data', data, 'name', 'md.inversion.cost_functions', 'format', 'StringArray')
|
|---|
| 197 | WriteData(fid, prefix, 'data', num_cost_functions, 'name', 'md.inversion.num_cost_functions', 'format', 'Integer')
|
|---|
| 198 | # }}}
|
|---|
![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
