source: issm/trunk/src/m/classes/adinversion.py@ 20500

Last change on this file since 20500 was 20500, checked in by Mathieu Morlighem, 9 years ago

merged trunk-jpl and trunk for revision 20497
File size: 10.8 KB
RevLine 
[19379]1"""
2== == == == == == == == == == == == == == == == == == ==
3Auto generated python script for ISSM: /home/andrei/issm/trunk-jpl/src/m/classes/adinversion.m
4Created on 2015-05-15 via translateToPy.py Ver 1.0 by andrei
5== == == == == == == == == == == == == == == == == == ==
6
7Matlab script conversion into python
8translateToPy.py Author: Michael Pellegrin
9translateToPy.py Date: 09/24/12
10== == == == == == == == == == == == == == == == == == ==
11"""
12
13from MatlabFuncs import *
14
15from EnumDefinitions import *
16from numpy import *
17
18# ADINVERSION class definition
19
20#
21
22# Usage:
23
24# adinversion=adinversion();
25
26
27
28class adinversion:
29 def __init__(self):
30 iscontrol = 0
31 control_parameters = float('Nan')
32 control_scaling_factors = float('Nan')
33 maxsteps = 0
34 maxiter = 0
35 dxmin = 0
36 gttol = 0
37 cost_functions = float('Nan')
38 cost_functions_coefficients = float('Nan')
39 min_parameters = float('Nan')
40 max_parameters = float('Nan')
41 vx_obs = float('Nan')
42 vy_obs = float('Nan')
43 vz_obs = float('Nan')
44 vel_obs = float('Nan')
45 thickness_obs = float('Nan')
46 surface_obs = float('Nan')
47
48 def setdefaultparameters(self):
49
50 self.control_parameters=['FrictionCoefficient']
51
52
53# Scaling factor for each control
54 self.control_scaling_factors=1
55
56# number of iterations
57 self.maxsteps=20
58 self.maxiter=40
59
60# several responses can be used:
61 self.cost_functions=['FrictionCoefficient']
62
63# m1qn3 parameters
64 self.dxmin = 0.1
65 self.gttol = 1e-4
66
67 return self
68
69 def checkconsistency(self, md, solution, analyses):
70
71# Early return
72 if not self.iscontrol:
73 return
74
75 if not IssmConfig('_HAVE_M1QN3_'):
76 md = checkmessage(md,['M1QN3 has not been installed, ISSM needs to be reconfigured and recompiled with AD'])
77
78
79 num_controls=numpy.numel(md.inversion.control_parameters)
80 num_costfunc=numpy.size(md.inversion.cost_functions,2)
81
82
83 md = checkfield(md,'fieldname','inversion.iscontrol','values',[0, 1])
84 md = checkfield(md,'fieldname','inversion.control_parameters','cell',1,'values',\
85 ['BalancethicknessThickeningRate' 'FrictionCoefficient' 'MaterialsRheologyBbar' 'DamageDbar',\
86 'Vx' 'Vy' 'Thickness' 'BalancethicknessOmega' 'BalancethicknessApparentMassbalance'])
87 md = checkfield(md,'fieldname','inversion.control_scaling_factors','size',[1, num_controls],'>',0,float('Nan'),1)
88 md = checkfield(md,'fieldname','inversion.maxsteps','numel',1,'>=',0)
89 md = checkfield(md,'fieldname','inversion.maxiter','numel',1,'>=',0)
90 md = checkfield(md,'fieldname','inversion.dxmin','numel',1,'>',0)
91 md = checkfield(md,'fieldname','inversion.gttol','numel',1,'>',0)
92 md = checkfield(md,'fieldname','inversion.cost_functions','size',[1, num_costfunc],'values', [i for i in range(101,106)]+[201]+[i for i in range(501,507)]+[i for i in range(601,605)]+[i for i in range(1001, 1011)])
93 md = checkfield(md,'fieldname','inversion.cost_functions_coefficients','size',[md.mesh.numberofvertices, num_costfunc],'>=',0)
94 md = checkfield(md,'fieldname','inversion.min_parameters','size',[md.mesh.numberofvertices, num_controls])
95 md = checkfield(md,'fieldname','inversion.max_parameters','size',[md.mesh.numberofvertices, num_controls])
96
97
98 if solution==BalancethicknessSolutionEnum():
99 md = checkfield(md,'fieldname','inversion.thickness_obs','size',[md.mesh.numberofvertices, 1],float('Nan'),1)
100 md = checkfield(md,'fieldname','inversion.surface_obs','size',[md.mesh.numberofvertices, 1], float('Nan'),1)
101 elif solution==BalancethicknessSoftSolutionEnum():
102 md = checkfield(md,'fieldname','inversion.thickness_obs','size',[md.mesh.numberofvertices, 1],float('Nan'),1)
103 else:
104 md = checkfield(md,'fieldname','inversion.vx_obs','size',[md.mesh.numberofvertices, 1],float('Nan'),1)
105 if not numpy.strcmp(domaintype(md.mesh),'2Dvertical'):
106 md = checkfield(md,'fieldname','inversion.vy_obs','size',[md.mesh.numberofvertices, 1],float('Nan'),1)
107 return md
108
109 def __repr__(self):
110 string = ' adinversion parameters:'
111 string ="%s\n\%s"%(string, fielddisplay(self,'iscontrol','is inversion activated?'))
112 string ="%s\n\%s"%(string, fielddisplay(self,'control_parameters','ex: [''FrictionCoefficient''], or [''MaterialsRheologyBbar'']'))
113 string ="%s\n\%s"%(string, fielddisplay(self,'control_scaling_factors','order of magnitude of each control (useful for multi-parameter optimization)'))
114 string ="%s\n\%s"%(string, fielddisplay(self,'maxsteps','maximum number of iterations (gradient computation)'))
115 string ="%s\n\%s"%(string, fielddisplay(self,'maxiter','maximum number of Function evaluation (forward run)'))
116 string ="%s\n\%s"%(string, fielddisplay(self,'dxmin','convergence criterion: two points less than dxmin from eachother (sup-norm) are considered identical'))
117 string ="%s\n\%s"%(string, fielddisplay(self,'gttol','convergence criterion: ||g(X)||/||g(X0)|| (g(X0): gradient at initial guess X0)'))
118 string ="%s\n\%s"%(string, fielddisplay(self,'cost_functions','indicate the type of response for each optimization step'))
119 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'))
120 string ="%s\n\%s"%(string, fielddisplay(self,'min_parameters','absolute minimum acceptable value of the inversed parameter on each vertex'))
121 string ="%s\n\%s"%(string, fielddisplay(self,'max_parameters','absolute maximum acceptable value of the inversed parameter on each vertex'))
122 string ="%s\n\%s"%(string, fielddisplay(self,'vx_obs','observed velocity x component [m/yr]'))
123 string ="%s\n\%s"%(string, fielddisplay(self,'vy_obs','observed velocity y component [m/yr]'))
124 string ="%s\n\%s"%(string, fielddisplay(self,'vel_obs','observed velocity magnitude [m/yr]'))
125 string ="%s\n\%s"%(string, fielddisplay(self,'thickness_obs','observed thickness [m]'))
126 string ="%s\n\%s"%(string, fielddisplay(self,'surface_obs','observed surface elevation [m]'))
127 string ="%s\n%s"%(string,'Available cost functions:');
128 string ="%s\n%s"%(string,' 101: SurfaceAbsVelMisfit');
129 string ="%s\n%s"%(string,' 102: SurfaceRelVelMisfit');
130 string ="%s\n%s"%(string,' 103: SurfaceLogVelMisfit');
131 string ="%s\n%s"%(string,' 104: SurfaceLogVxVyMisfit');
132 string ="%s\n%s"%(string,' 105: SurfaceAverageVelMisfit');
133 string ="%s\n%s"%(string,' 201: ThicknessAbsMisfit');
134 string ="%s\n%s"%(string,' 501: DragCoefficientAbsGradient');
135 string ="%s\n%s"%(string,' 502: RheologyBbarAbsGradient');
136 string ="%s\n%s"%(string,' 503: ThicknessAbsGradient');
137
138 return string
139
140 def marshall(self):
141
142 yts=365.0*24.0*3600.0;
143
144 WriteData(fid,'object',self,'class','inversion','fieldname','iscontrol','format','Boolean');
145 WriteData(fid,'enum',InversionTypeEnum(),'data',4,'format','Integer');
146 if not self.iscontrol:
147 return
148 WriteData(fid,'object',self,'class','inversion','fieldname','control_scaling_factors','format','DoubleMat','mattype',3);
149 WriteData(fid,'object',self,'class','inversion','fieldname','maxsteps','format','Integer');
150 WriteData(fid,'object',self,'class','inversion','fieldname','maxiter','format','Integer');
151 WriteData(fid,'object',self,'class','inversion','fieldname','dxmin','format','Double');
152 WriteData(fid,'object',self,'class','inversion','fieldname','gttol','format','Double');
153 WriteData(fid,'object',self,'class','inversion','fieldname','cost_functions_coefficients','format','DoubleMat','mattype',1);
154 WriteData(fid,'object',self,'class','inversion','fieldname','min_parameters','format','DoubleMat','mattype',3);
155 WriteData(fid,'object',self,'class','inversion','fieldname','max_parameters','format','DoubleMat','mattype',3);
156 WriteData(fid,'object',self,'class','inversion','fieldname','vx_obs','format','DoubleMat','mattype',1,'scale',1./yts);
157 WriteData(fid,'object',self,'class','inversion','fieldname','vy_obs','format','DoubleMat','mattype',1,'scale',1./yts);
158 WriteData(fid,'object',self,'class','inversion','fieldname','vz_obs','format','DoubleMat','mattype',1,'scale',1./yts);
159 if(numel(self.thickness_obs)==md.mesh.numberofelements):
160 mattype=2;
161 else:
162 mattype=1;
163
164 WriteData(fid,'object',self,'class','inversion','fieldname','thickness_obs','format','DoubleMat','mattype',mattype);
165 WriteData(fid,'object',self,'class','inversion','fieldname','surface_obs','format','DoubleMat','mattype',mattype);
166
167 #process control parameters
168 num_control_parameters = numpy.numel(self.control_parameters);
169 data = numpy.array([StringToEnum(self.control_parameter[0]) for control_parameter in self.control_parameters]).reshape(1,-1)
170
171 WriteData(fid,'data',data,'enum',InversionControlParametersEnum(),'format','DoubleMat','mattype',3);
172 WriteData(fid,'data',num_control_parameters,'enum',InversionNumControlParametersEnum(),'format','Integer');
173
174 #process cost functions
175 num_cost_functions=numpy.size(self.cost_functions,2);
176 data=copy.deepcopy(self.cost_functions)
177 data[numpy.nonzero(self.cost_functions==101)] =SurfaceAbsVelMisfitEnum();
178 data[numpy.nonzero(self.cost_functions==102)]=SurfaceRelVelMisfitEnum();
179 data[numpy.nonzero(self.cost_functions==103)]=SurfaceLogVelMisfitEnum();
180 data[numpy.nonzero(self.cost_functions==104)]=SurfaceLogVxVyMisfitEnum();
181 data[numpy.nonzero(self.cost_functions==105)]=SurfaceAverageVelMisfitEnum();
182 data[numpy.nonzero(self.cost_functions==201)]=ThicknessAbsMisfitEnum();
183 data[numpy.nonzero(self.cost_functions==501)]=DragCoefficientAbsGradientEnum();
184 data[numpy.nonzero(self.cost_functions==502)]=RheologyBbarAbsGradientEnum();
185 data[numpy.nonzero(self.cost_functions==503)]=ThicknessAbsGradientEnum();
186 data[numpy.nonzero(self.cost_functions==504)]=ThicknessAlongGradientEnum();
187 data[numpy.nonzero(self.cost_functions==505)]=ThicknessAcrossGradientEnum();
188 data[numpy.nonzero(self.cost_functions==506)]=BalancethicknessMisfitEnum();
189 data[numpy.nonzero(self.cost_functions==601)]=SurfaceAbsMisfitEnum();
190 data[numpy.nonzero(self.cost_functions==1001)]=Outputdefinition1Enum();
191 data[numpy.nonzero(self.cost_functions==1002)]=Outputdefinition2Enum();
192 data[numpy.nonzero(self.cost_functions==1003)]=Outputdefinition3Enum();
193 data[numpy.nonzero(self.cost_functions==1004)]=Outputdefinition4Enum();
194 data[numpy.nonzero(self.cost_functions==1005)]=Outputdefinition5Enum();
195 data[numpy.nonzero(self.cost_functions==1006)]=Outputdefinition6Enum();
196 data[numpy.nonzero(self.cost_functions==1007)]=Outputdefinition7Enum();
197 data[numpy.nonzero(self.cost_functions==1008)]=Outputdefinition8Enum();
198 data[numpy.nonzero(self.cost_functions==1009)]=Outputdefinition8Enum();
199 data[numpy.nonzero(self.cost_functions==1010)]=Outputdefinition10Enum();
200 WriteData(fid,'data',data,'enum',InversionCostFunctionsEnum(),'format','DoubleMat','mattype',3);
201 WriteData(fid,'data',num_cost_functions,'enum',InversionNumCostFunctionsEnum(),'format','Integer');
202
Note: See TracBrowser for help on using the repository browser.