Index: /issm/trunk-jpl/src/m/classes/adinversion.m =================================================================== --- /issm/trunk-jpl/src/m/classes/adinversion.m (revision 18875) +++ /issm/trunk-jpl/src/m/classes/adinversion.m (revision 18875) @@ -0,0 +1,192 @@ +%ADINVERSION class definition +% +% Usage: +% adinversion=adinversion(); + +classdef adinversion + properties (SetAccess=public) + iscontrol = 0 + control_parameters = NaN + control_scaling_factors = NaN + maxsteps = 0 + maxiter = 0 + dxmin = 0 + gttol = 0 + cost_functions = NaN + cost_functions_coefficients = NaN + min_parameters = NaN + max_parameters = NaN + vx_obs = NaN + vy_obs = NaN + vz_obs = NaN + vel_obs = NaN + thickness_obs = NaN + surface_obs = NaN + + end + methods + function obj = adinversion(varargin) % {{{ + switch nargin + case 0 + obj=setdefaultparameters(obj); + case 1 + obj=structtoobj(adinversion(),varargin{1}); + otherwise + error('constructor not supported'); + end + end % }}} + function self = setdefaultparameters(self) % {{{ + + %parameter to be inferred by control methods (only + %drag and B are supported yet) + self.control_parameters={'FrictionCoefficient'}; + + %Scaling factor for each control + self.control_scaling_factors=1; + + %number of iterations + self.maxsteps=20; + self.maxiter=40; + + %several responses can be used: + self.cost_functions={'FrictionCoefficient'}; + + %m1qn3 parameters + self.dxmin = 0.1; + self.gttol = 1e-4; + + end % }}} + function md = checkconsistency(obj,md,solution,analyses) % {{{ + + %Early return + if ~obj.iscontrol, return; end + + if ~IssmConfig('_HAVE_M1QN3_'), + md = checkmessage(md,['M1QN3 has not been installed, ISSM needs to be reconfigured and recompiled with AD']); + end + + num_controls=numel(md.inversion.control_parameters); + num_costfunc=size(md.inversion.cost_functions,2); + + md = checkfield(md,'fieldname','inversion.iscontrol','values',[0 1]); + md = checkfield(md,'fieldname','inversion.control_parameters','cell',1,'values',... + {'BalancethicknessThickeningRate' 'FrictionCoefficient' 'MaterialsRheologyBbar' 'DamageDbar',... + 'Vx' 'Vy' 'Thickness' 'BalancethicknessOmega' 'BalancethicknessApparentMassbalance'}); + md = checkfield(md,'fieldname','inversion.control_scaling_factors','size',[1 num_controls],'>',0,'NaN',1); + md = checkfield(md,'fieldname','inversion.maxsteps','numel',1,'>=',0); + md = checkfield(md,'fieldname','inversion.maxiter','numel',1,'>=',0); + md = checkfield(md,'fieldname','inversion.dxmin','numel',1,'>',0); + md = checkfield(md,'fieldname','inversion.gttol','numel',1,'>',0); + md = checkfield(md,'fieldname','inversion.cost_functions','size',[1 num_costfunc],'values',[101:105 201 501:506 601:604 1001:1010]); + md = checkfield(md,'fieldname','inversion.cost_functions_coefficients','size',[md.mesh.numberofvertices num_costfunc],'>=',0); + md = checkfield(md,'fieldname','inversion.min_parameters','size',[md.mesh.numberofvertices num_controls]); + md = checkfield(md,'fieldname','inversion.max_parameters','size',[md.mesh.numberofvertices num_controls]); + + if solution==BalancethicknessSolutionEnum() + md = checkfield(md,'fieldname','inversion.thickness_obs','size',[md.mesh.numberofvertices 1],'NaN',1); + md = checkfield(md,'fieldname','inversion.surface_obs','size',[md.mesh.numberofvertices 1],'NaN',1); + elseif solution==BalancethicknessSoftSolutionEnum() + md = checkfield(md,'fieldname','inversion.thickness_obs','size',[md.mesh.numberofvertices 1],'NaN',1); + else + md = checkfield(md,'fieldname','inversion.vx_obs','size',[md.mesh.numberofvertices 1],'NaN',1); + if ~strcmp(domaintype(md.mesh),'2Dvertical'), + md = checkfield(md,'fieldname','inversion.vy_obs','size',[md.mesh.numberofvertices 1],'NaN',1); + end + end + end % }}} + function disp(obj) % {{{ + disp(sprintf(' adinversion parameters:')); + fielddisplay(obj,'iscontrol','is inversion activated?'); + fielddisplay(obj,'control_parameters','ex: {''FrictionCoefficient''}, or {''MaterialsRheologyBbar''}'); + fielddisplay(obj,'control_scaling_factors','order of magnitude of each control (useful for multi-parameter optimization)'); + fielddisplay(obj,'maxsteps','maximum number of iterations (gradient computation)'); + fielddisplay(obj,'maxiter','maximum number of Function evaluation (forward run)'); + fielddisplay(obj,'dxmin','convergence criterion: two points less than dxmin from eachother (sup-norm) are considered identical'); + fielddisplay(obj,'gttol','convergence criterion: ||g(X)||/||g(X0)|| (g(X0): gradient at initial guess X0)'); + fielddisplay(obj,'cost_functions','indicate the type of response for each optimization step'); + fielddisplay(obj,'cost_functions_coefficients','cost_functions_coefficients applied to the misfit of each vertex and for each control_parameter'); + fielddisplay(obj,'min_parameters','absolute minimum acceptable value of the inversed parameter on each vertex'); + fielddisplay(obj,'max_parameters','absolute maximum acceptable value of the inversed parameter on each vertex'); + fielddisplay(obj,'vx_obs','observed velocity x component [m/yr]'); + fielddisplay(obj,'vy_obs','observed velocity y component [m/yr]'); + fielddisplay(obj,'vel_obs','observed velocity magnitude [m/yr]'); + fielddisplay(obj,'thickness_obs','observed thickness [m]'); + fielddisplay(obj,'surface_obs','observed surface elevation [m]'); + + disp('Available cost functions:'); + disp(' 101: SurfaceAbsVelMisfit'); + disp(' 102: SurfaceRelVelMisfit'); + disp(' 103: SurfaceLogVelMisfit'); + disp(' 104: SurfaceLogVxVyMisfit'); + disp(' 105: SurfaceAverageVelMisfit'); + disp(' 201: ThicknessAbsMisfit'); + disp(' 501: DragCoefficientAbsGradient'); + disp(' 502: RheologyBbarAbsGradient'); + disp(' 503: ThicknessAbsGradient'); + end % }}} + function marshall(obj,md,fid) % {{{ + + yts=365.0*24.0*3600.0; + + WriteData(fid,'object',obj,'class','inversion','fieldname','iscontrol','format','Boolean'); + WriteData(fid,'enum',InversionTypeEnum(),'data',4,'format','Integer'); + if ~obj.iscontrol, return; end + WriteData(fid,'object',obj,'class','inversion','fieldname','control_scaling_factors','format','DoubleMat','mattype',3); + WriteData(fid,'object',obj,'class','inversion','fieldname','maxsteps','format','Integer'); + WriteData(fid,'object',obj,'class','inversion','fieldname','maxiter','format','Integer'); + WriteData(fid,'object',obj,'class','inversion','fieldname','dxmin','format','Double'); + WriteData(fid,'object',obj,'class','inversion','fieldname','gttol','format','Double'); + WriteData(fid,'object',obj,'class','inversion','fieldname','cost_functions_coefficients','format','DoubleMat','mattype',1); + WriteData(fid,'object',obj,'class','inversion','fieldname','min_parameters','format','DoubleMat','mattype',3); + WriteData(fid,'object',obj,'class','inversion','fieldname','max_parameters','format','DoubleMat','mattype',3); + WriteData(fid,'object',obj,'class','inversion','fieldname','vx_obs','format','DoubleMat','mattype',1,'scale',1./yts); + WriteData(fid,'object',obj,'class','inversion','fieldname','vy_obs','format','DoubleMat','mattype',1,'scale',1./yts); + WriteData(fid,'object',obj,'class','inversion','fieldname','vz_obs','format','DoubleMat','mattype',1,'scale',1./yts); + if(numel(obj.thickness_obs)==md.mesh.numberofelements), + mattype=2; + else + mattype=1; + end + WriteData(fid,'object',obj,'class','inversion','fieldname','thickness_obs','format','DoubleMat','mattype',mattype); + WriteData(fid,'object',obj,'class','inversion','fieldname','surface_obs','format','DoubleMat','mattype',mattype); + + %process control parameters + num_control_parameters=numel(obj.control_parameters); + data=zeros(1,num_control_parameters); + for i=1:num_control_parameters, + data(i)=StringToEnum(obj.control_parameters{i}); + end + WriteData(fid,'data',data,'enum',InversionControlParametersEnum(),'format','DoubleMat','mattype',3); + WriteData(fid,'data',num_control_parameters,'enum',InversionNumControlParametersEnum(),'format','Integer'); + + %process cost functions + num_cost_functions=size(obj.cost_functions,2); + data=obj.cost_functions; + pos=find(obj.cost_functions==101); data(pos)=SurfaceAbsVelMisfitEnum(); + pos=find(obj.cost_functions==102); data(pos)=SurfaceRelVelMisfitEnum(); + pos=find(obj.cost_functions==103); data(pos)=SurfaceLogVelMisfitEnum(); + pos=find(obj.cost_functions==104); data(pos)=SurfaceLogVxVyMisfitEnum(); + pos=find(obj.cost_functions==105); data(pos)=SurfaceAverageVelMisfitEnum(); + pos=find(obj.cost_functions==201); data(pos)=ThicknessAbsMisfitEnum(); + pos=find(obj.cost_functions==501); data(pos)=DragCoefficientAbsGradientEnum(); + pos=find(obj.cost_functions==502); data(pos)=RheologyBbarAbsGradientEnum(); + pos=find(obj.cost_functions==503); data(pos)=ThicknessAbsGradientEnum(); + pos=find(obj.cost_functions==504); data(pos)=ThicknessAlongGradientEnum(); + pos=find(obj.cost_functions==505); data(pos)=ThicknessAcrossGradientEnum(); + pos=find(obj.cost_functions==506); data(pos)=BalancethicknessMisfitEnum(); + pos=find(obj.cost_functions==601); data(pos)=SurfaceAbsMisfitEnum(); + pos=find(obj.cost_functions==1001); data(pos)=Outputdefinition1Enum(); + pos=find(obj.cost_functions==1002); data(pos)=Outputdefinition2Enum(); + pos=find(obj.cost_functions==1003); data(pos)=Outputdefinition3Enum(); + pos=find(obj.cost_functions==1004); data(pos)=Outputdefinition4Enum(); + pos=find(obj.cost_functions==1005); data(pos)=Outputdefinition5Enum(); + pos=find(obj.cost_functions==1006); data(pos)=Outputdefinition6Enum(); + pos=find(obj.cost_functions==1007); data(pos)=Outputdefinition7Enum(); + pos=find(obj.cost_functions==1008); data(pos)=Outputdefinition8Enum(); + pos=find(obj.cost_functions==1009); data(pos)=Outputdefinition8Enum(); + pos=find(obj.cost_functions==1010); data(pos)=Outputdefinition10Enum(); + WriteData(fid,'data',data,'enum',InversionCostFunctionsEnum(),'format','DoubleMat','mattype',3); + WriteData(fid,'data',num_cost_functions,'enum',InversionNumCostFunctionsEnum(),'format','Integer'); + end % }}} + end +end