%INVERSION class definition
%
%   Usage:
%      inversion=inversion();

classdef inversion
	properties (SetAccess=public) 
		iscontrol                   = 0
		incomplete_adjoint          = 0
		control_parameters          = NaN
		nsteps                      = 0
		maxiter_per_step            = NaN
		cost_functions              = NaN
		cost_functions_coefficients = NaN
		gradient_scaling            = NaN
		cost_function_threshold     = 0
		min_parameters              = NaN
		max_parameters              = NaN
		step_threshold              = NaN
		vx_obs                      = NaN
		vy_obs                      = NaN
		vz_obs                      = NaN
		vel_obs                     = NaN
		thickness_obs               = NaN
	end
	methods
		function createxml(obj,fid) % {{{
			fprintf(fid, '<!-- inversion -->\n');            

			% inversion parameters
			fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="inversion parameters">','<section name="inversion" />');                    
			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="iscontrol" type="',class(obj.iscontrol),'" default="',convert2str(obj.iscontrol),'">','     <section name="inversion" />','     <help> is inversion activated? </help>','  </parameter>');

			% incompleteadjoing drop-down (0 or 1)
			fprintf(fid,'%s\n%s\n%s\n%s\n','  <parameter key ="incomplete_adjoint" type="alternative" optional="false">','     <section name="inversion" />','     <help> 1: linear viscosity, 0: non-linear viscosity </help>');
			fprintf(fid,'%s\n','       <option value="0" type="string" default="true"> </option>');
			fprintf(fid,'%s\n%s\n','       <option value="1" type="string" default="false"> </option>','</parameter>');

			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="control_parameters" type="',class(obj.control_parameters),'" default="',convert2str(obj.control_parameters),'">','     <section name="inversion" />','     <help> ex: {''FrictionCoefficient''}, or {''MaterialsRheologyBbar''} </help>','  </parameter>');

			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="nsteps" type="',class(obj.nsteps),'" default="',convert2str(obj.nsteps),'">','     <section name="inversion" />','     <help> number of optimization searches </help>','  </parameter>');
			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="cost_functions" type="',class(obj.cost_functions),'" default="',convert2str(obj.cost_functions),'">','     <section name="inversion" />','     <help> indicate the type of response for each optimization step  </help>','  </parameter>');
			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="cost_functions_coefficients" type="',class(obj.cost_functions_coefficients),'" default="',convert2str(obj.cost_functions_coefficients),'">','     <section name="inversion" />','     <help> cost_functions_coefficients applied to the misfit of each vertex and for each control_parameter </help>','  </parameter>');

			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="cost_function_threshold" type="',class(obj.cost_function_threshold),'" default="',convert2str(obj.cost_function_threshold),'">','     <section name="inversion" />','     <help> misfit convergence criterion. Default is 1%, NaN if not applied </help>','  </parameter>');
			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="maxiter_per_step" type="',class(obj.maxiter_per_step),'" default="',convert2str(obj.maxiter_per_step),'">','     <section name="inversion" />','     <help> maximum iterations during each optimization step  </help>','  </parameter>');
			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="gradient_scaling" type="',class(obj.gradient_scaling),'" default="',convert2str(obj.gradient_scaling),'">','     <section name="inversion" />','     <help> scaling factor on gradient direction during optimization, for each optimization step </help>','  </parameter>');

			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="step_threshold" type="',class(obj.step_threshold),'" default="',convert2str(obj.step_threshold),'">','     <section name="inversion" />','     <help> decrease threshold for misfit, default is 30% </help>','  </parameter>');
			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="min_parameters" type="',class(obj.min_parameters),'" default="',convert2str(obj.min_parameters),'">','     <section name="inversion" />','     <help> absolute minimum acceptable value of the inversed parameter on each vertex </help>','  </parameter>');
			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="max_parameters" type="',class(obj.max_parameters),'" default="',convert2str(obj.max_parameters),'">','     <section name="inversion" />','     <help> absolute maximum acceptable value of the inversed parameter on each vertex </help>','  </parameter>');

			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="vx_obs" type="',class(obj.vx_obs),'" default="',convert2str(obj.vx_obs),'">','     <section name="inversion" />','     <help> observed velocity x component [m/yr] </help>','  </parameter>');
			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="vy_obs" type="',class(obj.vy_obs),'" default="',convert2str(obj.vy_obs),'">','     <section name="inversion" />','     <help> observed velocity y component [m/yr]  </help>','  </parameter>');
			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="vel_obs" type="',class(obj.vel_obs),'" default="',convert2str(obj.vel_obs),'">','     <section name="inversion" />','     <help> observed velocity magnitude [m/yr] </help>','  </parameter>');
			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="thickness_obs" type="',class(obj.thickness_obs),'" default="',convert2str(obj.thickness_obs),'">','     <section name="inversion" />','     <help> observed thickness [m]) </help>','  </parameter>');

			fprintf(fid,'%s\n%s\n','</frame>');    

			fprintf(fid,'%s\n%s\n%s\n','<frame key="2" label="Available cost functions">','<section name="inversion" />');                    
			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="SurfaceAbsVelMisfit" type="','string','" default="','101','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="SurfaceRelVelMisfit" type="','string','" default="','102','">','     <section name="inversion" />','     <help>   </help>','  </parameter>');
			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="SurfaceLogVelMisfit" type="','string','" default="','103','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');

			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="SurfaceLogVxVyMisfit" type="','string','" default="','104','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="SurfaceAverageVelMisfit" type="','string','" default="','105','">','     <section name="inversion" />','     <help>   </help>','  </parameter>');
			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="ThicknessAbsMisfit" type="','string','" default="','106','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');

			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="DragCoefficientAbsGradient" type="','string','" default="','107','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="RheologyBbarAbsGradient" type="','string','" default="','108','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="ThicknessAbsGradient" type="','string','" default="','109','">','     <section name="inversion" />','     <help> </help>','  </parameter>');

			fprintf(fid,'%s\n%s\n','</frame>');    

		end % }}}       
		function obj = inversion(varargin) % {{{
			switch nargin
				case 0
					obj=setdefaultparameters(obj);
				case 1
					obj=structtoobj(inversion(),varargin{1});
				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;

			%maximum number of iteration in the optimization algorithm for
			%each step
			obj.maxiter_per_step=20*ones(obj.nsteps,1);

			%the inversed parameter is updated as follows:
			%new_par=old_par + gradient_scaling(n)*C*gradient with C in [0 1];
			%usually the gradient_scaling must be of the order of magnitude of the 
			%inversed parameter (10^8 for B, 50 for drag) and can be decreased
			%after the first iterations
			obj.gradient_scaling=50*ones(obj.nsteps,1);

			%several responses can be used:
			obj.cost_functions=101;

			%step_threshold is used to speed up control method. When
			%misfit(1)/misfit(0) < obj.step_threshold, we go directly to
			%the next step
			obj.step_threshold=.7*ones(obj.nsteps,1); %30 per cent decrement.

			%cost_function_threshold is a criteria to stop the control methods.
			%if J[n]-J[n-1]/J[n] < criteria, the control run stops
			%NaN if not applied
			obj.cost_function_threshold=NaN; %not activated

		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' 'DamageDbar' 'Vx' 'Vy' 'Thickness' 'BalancethicknessNux' 'BalancethicknessNuy' 'BalancethicknessApparentMassbalance'});
			md = checkfield(md,'fieldname','inversion.nsteps','numel',1,'>=',0);
			md = checkfield(md,'fieldname','inversion.maxiter_per_step','size',[md.inversion.nsteps 1],'>=',0);
			md = checkfield(md,'fieldname','inversion.step_threshold','size',[md.inversion.nsteps 1]);
			md = checkfield(md,'fieldname','inversion.cost_functions','size',[1 num_costfunc],'values',[101:105 201 501:507]);
			md = checkfield(md,'fieldname','inversion.cost_functions_coefficients','size',[md.mesh.numberofvertices num_costfunc],'>=',0);
			md = checkfield(md,'fieldname','inversion.gradient_scaling','size',[md.inversion.nsteps num_controls]);
			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]);

			%Only SSA, HO and FS are supported right now
			if solution==StressbalanceSolutionEnum()
				if ~(md.flowequation.isSSA || md.flowequation.isHO || md.flowequation.isFS),
					md = checkmessage(md,['inversion can only be performed for SSA, HO or FS ice flow models']);
				end
			end

			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);
				md = checkfield(md,'fieldname','inversion.vy_obs','size',[md.mesh.numberofvertices 1],'NaN',1);
			end
		end % }}}
		function disp(obj) % {{{
			disp(sprintf('   inversion 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,'cost_function_threshold','misfit convergence criterion. Default is 1%, NaN if not applied');
			fielddisplay(obj,'maxiter_per_step','maximum iterations during each optimization step');
			fielddisplay(obj,'gradient_scaling','scaling factor on gradient direction during optimization, for each optimization step');
			fielddisplay(obj,'step_threshold','decrease threshold for misfit, default is 30%');
			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,'enum',InversionTypeEnum(),'data',0,'format','Integer');
			WriteData(fid,'object',obj,'fieldname','iscontrol','format','Boolean');
			WriteData(fid,'object',obj,'fieldname','incomplete_adjoint','format','Boolean');
			if ~obj.iscontrol, return; end
			WriteData(fid,'object',obj,'fieldname','nsteps','format','Integer');
			WriteData(fid,'object',obj,'fieldname','maxiter_per_step','format','IntMat','mattype',3);
			WriteData(fid,'object',obj,'fieldname','cost_functions_coefficients','format','DoubleMat','mattype',1);
			WriteData(fid,'object',obj,'fieldname','gradient_scaling','format','DoubleMat','mattype',3);
			WriteData(fid,'object',obj,'fieldname','cost_function_threshold','format','Double');
			WriteData(fid,'object',obj,'fieldname','min_parameters','format','DoubleMat','mattype',3);
			WriteData(fid,'object',obj,'fieldname','max_parameters','format','DoubleMat','mattype',3);
			WriteData(fid,'object',obj,'fieldname','step_threshold','format','DoubleMat','mattype',3);
			WriteData(fid,'object',obj,'fieldname','vx_obs','format','DoubleMat','mattype',1,'scale',1./yts);
			WriteData(fid,'object',obj,'fieldname','vy_obs','format','DoubleMat','mattype',1,'scale',1./yts);
			WriteData(fid,'object',obj,'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);


			%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==507); data(pos)=Balancethickness2MisfitEnum();
			WriteData(fid,'data',data,'enum',InversionCostFunctionsEnum(),'format','DoubleMat','mattype',3);
			WriteData(fid,'data',num_cost_functions,'enum',InversionNumCostFunctionsEnum(),'format','Integer');
		end % }}}
	end
end
