%TAOINVERSION class definition
%
%   Usage:
%      taoinversion=taoinversion();

classdef taoinversion
	properties (SetAccess=public) 
		iscontrol                   = 0
		incomplete_adjoint          = 0
		control_parameters          = NaN
		nsteps                      = 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
	end
	methods
		function obj = taoinversion(varargin) % {{{
			switch nargin
				case 0
					obj=setdefaultparameters(obj);
				case 1
					if isa(varargin{1},'inversion'),
						disp('converting inversion to taoinversion');
						in=varargin{1};
						obj.iscontrol                   = in.iscontrol;
						obj.incomplete_adjoint          = in.incomplete_adjoint;
						obj.control_parameters          = in.control_parameters;
						obj.nsteps                      = in.nsteps;
						obj.cost_functions              = in.cost_functions(1,:); %Keep first line only
						obj.cost_functions_coefficients = in.cost_functions_coefficients;
						obj.min_parameters              = in.min_parameters;
						obj.max_parameters              = in.max_parameters;
						obj.vx_obs                      = in.vx_obs;
						obj.vy_obs                      = in.vy_obs;
						obj.vz_obs                      = in.vz_obs;
						obj.vel_obs                     = in.vel_obs;
						obj.thickness_obs               = in.thickness_obs;
					end
				otherwise
					error('constructor not supported');
			end
		end % }}}
		function obj = setdefaultparameters(obj) % {{{

			%default is incomplete adjoint for now
			obj.incomplete_adjoint=1;

			%parameter to be inferred by control methods (only
			%drag and B are supported yet)
			obj.control_parameters={'FrictionCoefficient'};

			%number of steps in the control methods
			obj.nsteps=20;

			%several responses can be used:
			obj.cost_functions=101*ones(obj.nsteps,1);

		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,'inversion.iscontrol','values',[0 1]);
			md = checkfield(md,'inversion.incomplete_adjoint','values',[0 1]);
			md = checkfield(md,'inversion.control_parameters','cell',1,'values',...
				{'BalancethicknessThickeningRate' 'FrictionCoefficient' 'MaterialsRheologyBbar' 'MaterialsRheologyZbar' 'Vx' 'Vy' 'Thickness'});
			md = checkfield(md,'inversion.nsteps','numel',1,'>=',1);
			md = checkfield(md,'inversion.cost_functions','size',[1 num_costfunc],'values',[101:105 201 501:506]);
			md = checkfield(md,'inversion.cost_functions_coefficients','size',[md.mesh.numberofvertices num_costfunc],'>=',0);
			md = checkfield(md,'inversion.min_parameters','size',[md.mesh.numberofvertices num_controls]);
			md = checkfield(md,'inversion.max_parameters','size',[md.mesh.numberofvertices num_controls]);

			if solution==BalancethicknessSolutionEnum()
				md = checkfield(md,'inversion.thickness_obs','size',[md.mesh.numberofvertices 1],'NaN',1);
			elseif solution==BalancethicknessSoftSolutionEnum()
				md = checkfield(md,'inversion.thickness_obs','size',[md.mesh.numberofvertices 1],'NaN',1);
			else
				md = checkfield(md,'inversion.vx_obs','size',[md.mesh.numberofvertices 1],'NaN',1);
				md = checkfield(md,'inversion.vy_obs','size',[md.mesh.numberofvertices 1],'NaN',1);
			end
		end % }}}
		function disp(obj) % {{{
			disp(sprintf('   taoinversion 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,'nsteps','number of optimization searches');
			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',InversionTaoEnum(),'data',true,'format','Boolean');
			if ~obj.iscontrol, return; end
			WriteData(fid,'object',obj,'class','inversion','fieldname','incomplete_adjoint','format','Boolean');
			WriteData(fid,'object',obj,'class','inversion','fieldname','nsteps','format','Integer');
			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(data==101); data(pos)=SurfaceAbsVelMisfitEnum();
			pos=find(data==102); data(pos)=SurfaceRelVelMisfitEnum();
			pos=find(data==103); data(pos)=SurfaceLogVelMisfitEnum();
			pos=find(data==104); data(pos)=SurfaceLogVxVyMisfitEnum();
			pos=find(data==105); data(pos)=SurfaceAverageVelMisfitEnum();
			pos=find(data==201); data(pos)=ThicknessAbsMisfitEnum();
			pos=find(data==501); data(pos)=DragCoefficientAbsGradientEnum();
			pos=find(data==502); data(pos)=RheologyBbarAbsGradientEnum();
			pos=find(data==503); data(pos)=ThicknessAbsGradientEnum();
			pos=find(data==504); data(pos)=ThicknessAlongGradientEnum();
			pos=find(data==505); data(pos)=ThicknessAcrossGradientEnum();
			pos=find(data==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
