Index: /issm/trunk-jpl/src/m/classes/inversionvalidation.m =================================================================== --- /issm/trunk-jpl/src/m/classes/inversionvalidation.m (revision 18008) +++ /issm/trunk-jpl/src/m/classes/inversionvalidation.m (revision 18008) @@ -0,0 +1,140 @@ +%INVERSIONVALIDATION class definition +% +% Usage: +% inversionvalidation=inversionvalidation(); + +classdef inversionvalidation + properties (SetAccess=public) + iscontrol = 0 + incomplete_adjoint = 0 + control_parameters = NaN + 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 + end + methods + function obj = inversionvalidation(varargin) % {{{ + switch nargin + case 0 + obj=setdefaultparameters(obj); + case 1 + obj=structtoobj(inversionvalidation(),varargin{1}); + otherwise + error('constructor not supported'); + end + end % }}} + function self = setdefaultparameters(self) % {{{ + + %default is incomplete adjoint for now + self.incomplete_adjoint=1; + + %parameter to be inferred by control methods (only + %drag and B are supported yet) + self.control_parameters={'FrictionCoefficient'}; + + %several responses can be used: + self.cost_functions=101; + end % }}} + function md = checkconsistency(obj,md,solution,analyses) % {{{ + + %Early return + if ~obj.iscontrol, return; 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.incomplete_adjoint','values',[0 1]); + md = checkfield(md,'fieldname','inversion.control_parameters','cell',1,'values',... + {'BalancethicknessThickeningRate' 'FrictionCoefficient' 'MaterialsRheologyBbar' 'MaterialsRheologyZbar' 'Vx' 'Vy' 'Thickness'}); + md = checkfield(md,'fieldname','inversion.cost_functions','size',[1 num_costfunc],'values',[101:105 201 501:506]); + 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); + 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(' inversionvalidation parameters:')); + fielddisplay(obj,'iscontrol','is inversion activated?'); + fielddisplay(obj,'incomplete_adjoint','1: linear viscosity, 0: non-linear viscosity'); + fielddisplay(obj,'control_parameters','ex: {''FrictionCoefficient''}, or {''MaterialsRheologyBbar''}'); + 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]'); + 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',3,'format','Integer'); + if ~obj.iscontrol, return; end + WriteData(fid,'object',obj,'class','inversion','fieldname','incomplete_adjoint','format','Boolean'); + 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); + WriteData(fid,'object',obj,'class','inversion','fieldname','thickness_obs','format','DoubleMat','mattype',1); + + %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(); + WriteData(fid,'data',data,'enum',InversionCostFunctionsEnum(),'format','DoubleMat','mattype',3); + WriteData(fid,'data',num_cost_functions,'enum',InversionNumCostFunctionsEnum(),'format','Integer'); + end % }}} + end +end