%SMBautoregression Class definition
%
%   Usage:
%      SMBautoregression=SMBautoregression();

classdef SMBautoregression
	properties (SetAccess=public)
		num_basins        = 0;
		beta0             = NaN;
		beta1             = NaN;
		ar_order          = 0;
		ar_initialtime    = 0;
		ar_timestep       = 0;
		phi               = NaN;
		basin_id          = NaN;
		lapserates        = NaN;
		elevationbins     = NaN;
		refelevation      = NaN;
		steps_per_step    = 1;
		averaging         = 0;
		requested_outputs = {};
	end
	methods
		function self = SMBautoregression(varargin) % {{{
			switch nargin
				case 0
					self=setdefaultparameters(self);
				otherwise
					error('constructor not supported');
			end
		end % }}}
		function self = extrude(self,md) % {{{
			%Nothing for now
		end % }}}
		function list = defaultoutputs(self,md) % {{{
			list = {''};
		end % }}}
		function self = initialize(self,md) % {{{
			if isnan(self.beta1)
				self.beta1 = zeros(1,self.num_basins); %no trend in SMB
				disp('      smb.beta1 (trend) not specified: value set to 0');
			end
			if (self.ar_order==0)
				self.ar_order = 1; %dummy 1 value for autoregression
				self.phi      = zeros(self.num_basins,self.ar_order); %autoregression coefficients all set to 0 
				disp('      smb.ar_order (order of autoregressive model) not specified: order of autoregressive model set to 0');
			end
			if (self.ar_initialtime==0)
				self.ar_initialtime = md.timestepping.start_time; %autoregression model has no prescribed initial time
				disp('      smb.ar_initialtime (initial time in the autoregressive model parameterization) not specified: set to md.timestepping.start_time');
			end
			if (self.ar_timestep==0)
				self.ar_timestep = md.timestepping.time_step; %autoregression model has no prescribed time step
				disp('      smb.ar_timestep (timestep of autoregressive model) not specified: set to md.timestepping.time_step');
			end
			if isnan(self.phi)
				self.phi = zeros(self.num_basins,self.ar_order); %autoregression model of order 0 
				disp('      smb.phi (lag coefficients) not specified: order of autoregressive model set to 0');
			end
		end % }}}
		function self = setdefaultparameters(self) % {{{
			self.ar_order    = 0.0; %autoregression model of order 0
		end % }}}
		function md = checkconsistency(self,md,solution,analyses) % {{{

			if ismember('MasstransportAnalysis',analyses),
				md = checkfield(md,'fieldname','smb.num_basins','numel',1,'NaN',1,'Inf',1,'>',0);
				md = checkfield(md,'fieldname','smb.basin_id','Inf',1,'>=',0,'<=',md.smb.num_basins,'size',[md.mesh.numberofelements,1]);
				md = checkfield(md,'fieldname','smb.beta0','NaN',1,'Inf',1,'size',[1,md.smb.num_basins],'numel',md.smb.num_basins); %scheme fails if passed as column vector
				md = checkfield(md,'fieldname','smb.beta1','NaN',1,'Inf',1,'size',[1,md.smb.num_basins],'numel',md.smb.num_basins); %scheme fails if passed as column vector
				md = checkfield(md,'fieldname','smb.ar_order','numel',1,'NaN',1,'Inf',1,'>=',0);
				md = checkfield(md,'fieldname','smb.ar_initialtime','numel',1,'NaN',1,'Inf',1); 
				md = checkfield(md,'fieldname','smb.ar_timestep','numel',1,'NaN',1,'Inf',1,'>=',md.timestepping.time_step); %autoregression time step cannot be finer than ISSM timestep
				md = checkfield(md,'fieldname','smb.phi','NaN',1,'Inf',1,'size',[md.smb.num_basins,md.smb.ar_order]);

				if (any(isnan(md.smb.refelevation)==0) || numel(md.smb.refelevation)>1)
               md = checkfield(md,'fieldname','smb.refelevation','NaN',1,'Inf',1,'>=',0,'size',[1,md.smb.num_basins],'numel',md.smb.num_basins);
            end
				[nbas,nbins] = size(md.smb.lapserates);
				if (any(isnan(reshape(md.smb.lapserates,[1,nbas*nbins]))==0) || numel(md.smb.lapserates)>1)
					md = checkfield(md,'fieldname','smb.lapserates','NaN',1,'Inf',1,'size',[md.smb.num_basins,nbins],'numel',md.smb.num_basins*nbins);
					md = checkfield(md,'fieldname','smb.elevationbins','NaN',1,'Inf',1,'size',[md.smb.num_basins,nbins-1],'numel',md.smb.num_basins*(nbins-1));
					if(issorted(md.smb.elevationbins,2)==0)
						error('md.smb.elevationbins should have rows in order of increasing elevation');
					end
				elseif (isnan(md.smb.elevationbins(1,1))==0 || numel(md.smb.elevationbins)>1)
					%elevationbins specified but not lapserates: this will inevitably lead to inconsistencies
					[nbas,nbins] = size(md.smb.elevationbins);
					nbins        = nbins+1;
					md = checkfield(md,'fieldname','smb.lapserates','NaN',1,'Inf',1,'size',[md.smb.num_basins,nbins],'numel',md.smb.num_basins*nbins);
					md = checkfield(md,'fieldname','smb.elevationbins','NaN',1,'Inf',1,'size',[md.smb.num_basins,nbins-1],'numel',md.smb.num_basins*(nbins-1));
				end
			end
			md = checkfield(md,'fieldname','smb.steps_per_step','>=',1,'numel',[1]);
			md = checkfield(md,'fieldname','smb.averaging','numel',[1],'values',[0 1 2]);
			md = checkfield(md,'fieldname','smb.requested_outputs','stringrow',1);
		end % }}}
		function disp(self) % {{{
			disp(sprintf('   surface forcings parameters:'));
			fielddisplay(self,'num_basins','number of different basins [unitless]');
			fielddisplay(self,'basin_id','basin number assigned to each element [unitless]');
			fielddisplay(self,'beta0','basin-specific intercept values [m ice eq./yr] (if beta_1==0 mean=beta_0/(1-sum(phi)))');
			fielddisplay(self,'beta1','basin-specific trend values [m ice eq. yr^(-2)]');
			fielddisplay(self,'ar_order','order of the autoregressive model [unitless]');
			fielddisplay(self,'ar_initialtime','initial time assumed in the autoregressive model parameterization [yr]');
			fielddisplay(self,'ar_timestep','time resolution of the autoregressive model [yr]');
			fielddisplay(self,'phi','basin-specific vectors of lag coefficients [unitless]');
			fielddisplay(self,'lapserates','basin-specific SMB lapse rates applied in each elevation bin, 1 row per basin, 1 column per bin [m ice eq yr^-1 m^-1] (default: no lapse rate)');
			fielddisplay(self,'elevationbins','basin-specific separations between elevation bins, 1 row per basin, 1 column per limit between bins [m] (default: no basin separation)');
			fielddisplay(self,'refelevation','basin-specific reference elevations at which SMB is calculated, and from which SMB is downscaled using lapserates (default: basin mean elevation) [m]');
			fielddisplay(self, 'steps_per_step', 'number of smb steps per time step');
			fielddisplay(self, 'averaging', 'averaging methods from short to long steps');
			disp(sprintf('%51s  0: Arithmetic (default)',' '));
			disp(sprintf('%51s  1: Geometric',' '));
			disp(sprintf('%51s  2: Harmonic',' '));
			fielddisplay(self,'requested_outputs','additional outputs requested');

		end % }}}
		function marshall(self,prefix,md,fid) % {{{

			yts=md.constants.yts;

			templapserates    = md.smb.lapserates;
			tempelevationbins = md.smb.elevationbins;
			temprefelevation  = md.smb.refelevation;
			[nbas,nbins]      = size(md.smb.lapserates);
			if(any(isnan(reshape(md.smb.lapserates,[1,nbas*nbins]))))
				templapserates = zeros(md.smb.num_basins,2);
				disp('      smb.lapserates not specified: set to 0');
			   tempelevationbins = zeros(md.smb.num_basins,1); %dummy elevation bins
			end
			if(any(isnan(md.smb.refelevation)))
				temprefelevation = zeros(1,md.smb.num_basins);
				areas = GetAreas(md.mesh.elements,md.mesh.x,md.mesh.y);
				for ii=1:md.smb.num_basins
					indices = find(md.smb.basin_id==ii);
					elemsh  = zeros(numel(indices),1);
					for jj=1:numel(indices)
						elemsh(jj) = mean(md.geometry.surface(md.mesh.elements(indices(jj),:)));
					end
					temprefelevation(ii) = sum(areas(indices).*elemsh)/sum(areas(indices));
				end
				if(any(reshape(md.smb.lapserates,[1,nbas*nbins])~=0))
					disp('      smb.refelevation not specified: Reference elevations set to mean surface elevation of basins');
				end
			end
			[nbas,nbins] = size(templapserates);

			WriteData(fid,prefix,'name','md.smb.model','data',13,'format','Integer');
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','num_basins','format','Integer');
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','ar_order','format','Integer');
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','ar_initialtime','format','Double','scale',yts);
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','ar_timestep','format','Double','scale',yts);
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','basin_id','data',self.basin_id-1,'name','md.smb.basin_id','format','IntMat','mattype',2); %0-indexed
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','beta0','format','DoubleMat','name','md.smb.beta0','scale',1./yts,'yts',yts);
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','beta1','format','DoubleMat','name','md.smb.beta1','scale',1./(yts^2),'yts',yts);
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','phi','format','DoubleMat','name','md.smb.phi','yts',yts); 
			WriteData(fid,prefix,'data',templapserates,'format','DoubleMat','name','md.smb.lapserates','scale',1./yts,'yts',yts);
			WriteData(fid,prefix,'data',tempelevationbins,'format','DoubleMat','name','md.smb.elevationbins');
			WriteData(fid,prefix,'data',temprefelevation,'format','DoubleMat','name','md.smb.refelevation');
			WriteData(fid,prefix,'data',nbins,'format','Integer','name','md.smb.num_bins');
			WriteData(fid,prefix,'object',self,'fieldname','steps_per_step','format','Integer');
			WriteData(fid,prefix,'object',self,'fieldname','averaging','format','Integer');

			%process requested outputs
			outputs = self.requested_outputs;
			pos  = find(ismember(outputs,'default'));
			if ~isempty(pos),
				outputs(pos) = [];                         %remove 'default' from outputs
				outputs      = [outputs defaultoutputs(self,md)]; %add defaults
			end
			WriteData(fid,prefix,'data',outputs,'name','md.smb.requested_outputs','format','StringArray');

		end % }}}
	end
end
