Index: ../trunk-jpl/src/m/classes/adinversion.m
===================================================================
--- ../trunk-jpl/src/m/classes/adinversion.m	(revision 0)
+++ ../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
