function this = convert(this,freq2,range,varargin)
% convert  Convert tseries object to a different frequency.
%
% Syntax
% =======
%
%     Y = convert(X,NEWFREQ)
%     Y = convert(X,NEWFREQ,RANGE,...)
%
% Input arguments
% ================
%
% * `X` [ tseries ] - Input tseries object that will be converted to a new
% frequency, `freq`, aggregating or intrapolating the data.
%
% * `NEWFREQ` [ numeric | char ] - New frequency to which the input data
% will be converted: `1` or `'A'` for annual, `2` or `'H'` for half-yearly,
% `4` or `'Q'` for quarterly, `6` or `'B'` for bi-monthly, and `12` or
% `'M'` for monthly.
%
% * `RANGE` [ numeric ] - Date range on which the input data will be
% converted.
%
% Output arguments
% =================
%
% * `Y` [ tseries ] - Output tseries created by converting `X` to the new
% frequency.
%
% Options
% ========
%
% * `'ignoreNaN='` [ `true` | *`false`* ] - Exclude NaNs from agreggation.
%
% * `'missing='` [ numeric | *`NaN`* | `'last'` ] - Replace missing
% observations with this value.
%
% Options for high- to low-frequency conversion (aggregation)
% ============================================================
%
% * `'method='` [ function_handle | `'first'` | `'last'` | *`@mean`* ] -
% Method that will be used to aggregate the high frequency data.
%
% * `'select='` [ numeric | *`Inf`* ] - Select only these high-frequency
% observations within each low-frequency period; Inf means all observations
% will be used.
%
% Options for low- to high-frequency conversion (interpolation)
% ==============================================================
%
% * `'method='` [ char | *`'cubic'`* | `'quadsum'` | `'quadavg'` ] -
% Interpolation method; any option available in the built-in `interp1`
% function can be used.
%
% * `'position='` [ *`'centre'`* | `'start'` | `'end'` ] - Position of the
% low-frequency date grid.
%
% Description
% ============
%
% The function handle that you pass in through the 'method' option when you
% aggregate the data (convert higher frequency to lower frequency) should
% behave like the built-in functions `mean`, `sum` etc. In other words, it
% is expected to accept two input arguments:
%
% * the data to be aggregated,
% * the dimension along which the aggregation is calculated.
%
% The function will be called with the second input argument set to 1, as
% the data are processed en block columnwise. If this call fails, `convert`
% will attempt to call the function with just one input argument, the data,
% but this is not a safe option under some circumstances since dimension
% mismatch may occur.
%
% Example
% ========

% -IRIS Toolbox.
% -Copyright (c) 2007-2012 Jaromir Benes.

freq1 = datfreq(this.start);
freq2 = recognisefreq_(freq2);

if isempty(freq2)
    utils.error('tseries', ...
        'Cannot determine requested output frequency.');
end

if nargin < 3
    range = Inf;
end

if freq1 == freq2
    % No conversion.
    if isequal(range,Inf)
        return
    else
        this.data = rangedata(this,range);
        this.start = range(1);
        this = mytrim(this);
        return
    end
elseif freq1 == 0
    % Conversion of daily series.
    options = passvalopt('tseries.convertdaily',varargin{:});
    call = @dailyaggregate_;
else
    % Conversion of Y, Z, Q, B, or M series.
    if freq1 > freq2
        % Aggregate.
        options = passvalopt('tseries.convertaggregate',varargin{:});
        if ~isempty(options.function)
            options.method = options.function;
        end
        call = @aggregate_;
    else
        % Interpolate.
        options = passvalopt('tseries.convertinterp',varargin{:});
        if any(strcmpi(options.method,{'quadsum','quadavg'}))
            % Quadratic interpolation matching sum or average.
            call = @interpmatch_;
        else
            % Built-in interp1.
            call = @interp_;
        end
    end
end

%**************************************************************************

this = call(this,range,freq1,freq2,options);

end

% Subfunctions follow.

% $ ***********************************************************************
function freq = recognisefreq_(freq)
freqNum = [1,2,4,6,12];
if ischar(freq)
    if ~isempty(freq)
        freqLetter = 'yzqbm';
        freq = lower(freq(1));
        if freq == 'a'
            % Dual options for annual frequency: Y or A.
            freq = 'y';
        elseif freq == 's'
            % Dual options for semi-annual frequency: Z or S.
            freq = 'z';
        end
        freq = freqNum(freq == freqLetter);
    else
        freq = [];
    end
elseif ~any(freq == freqNum)
    freq = [];
end
end
% $ recognisefreq_().

% $ ***********************************************************************
function x = aggregate_(x,range,fromfreq,tofreq,options)
if isa(options.method,'function_handle') ...
        && any(strcmp(char(options.method),{'first','last'}))
    options.method = char(options.method);
end

if ischar(options.method)
    options.method = str2func(options.method);
end

if isnan(x.start) && isempty(x.data)
    return
end

if isempty(range)
    x = empty(x);
    return
end

if ~any(isinf(range))
    x = resize(x,range);
end

datfunc = {@yy,@zz,@qq,@bb,@mm};
fromdate = datfunc{fromfreq == [1,2,4,6,12]};
todate = datfunc{tofreq == [1,2,4,6,12]};

startyear = dat2ypf(get(x,'start'));
endyear = dat2ypf(get(x,'end'));

fromdata = mygetdata(x,fromdate(startyear,1):fromdate(endyear,fromfreq));
fromdatasize = size(fromdata);
nper = fromdatasize(1);
fromdata = fromdata(:,:);
nfromdata = size(fromdata,2);
factor = fromfreq/tofreq;
todata = nan(nper/factor,nfromdata);
for i = 1 : size(fromdata,2)
    tmpdata = reshape(fromdata(:,i),[factor,nper/factor]);
    if ~isequal(options.select,Inf)
        tmpdata = tmpdata(options.select,:);
    end
    if options.ignorenan && any(isnan(fromdata(:,i)))
        for j = 1 : size(tmpdata,2)
            index = ~isnan(tmpdata(:,j));
            if any(index)
                try
                    todata(j,i) = options.method(tmpdata(index,j),1);
                catch %#ok<CTCH>
                    todata(j,i) = NaN;
                end
            else
                todata(j,i) = NaN;
            end
        end
    else
        try
            todata(:,i) = options.method(tmpdata,1);
        catch %#ok<CTCH>
            try %#ok<TRYNC>
                todata(:,i) = options.method(tmpdata);
            end
        end
    end
end
todata = reshape(todata,[nper/factor,fromdatasize(2:end)]);

x.start = todate(startyear,1);
x.data = todata;
x = mytrim(x);
end
% $ aggregate_().

% $ ***********************************************************************
% Subfunction first().
function x = first(x,varargin)
x = x(1,:);
end
% $ first().

% $ ***********************************************************************
function x = last(x,varargin)
x = x(end,:);
end
% $ last().


% $ ***********************************************************************
function this = ...
    dailyaggregate_(this,userrange,freq1,freq2,options) %#ok<INUSL>
if ischar(options.method)
    options.method = str2func(options.method);
end

if any(isinf(userrange))
    userrange = this.start + (0 : size(this.data,1)-1);
else
    userrange = userrange(1) : userrange(end);
end

periodfcn = @(month) ceil(freq2*month/12);
datefcn = @(year,period) datcode(freq2,year,period);

start = this.start;
data = this.data;

range = start + (0 : size(data,1)-1);
if isempty(range)
    return
end

tmpsize = size(data);
data = data(:,:);

tmp = datevec(userrange);
useryear = tmp(:,1);
userperiod = periodfcn(tmp(:,2));

tmp = datevec(range);
year = tmp(:,1);
period = periodfcn(tmp(:,2));

% Treat missing observations.
for t = 2 : size(data,1)
    index = isnan(data(t,:));
    if any(index)
        switch options.missing
            case 'last'
                data(t,index) = data(t-1,index);
            otherwise
                data(t,index) = options.missing;
        end
    end
end

start2 = datefcn(useryear(1),userperiod(1));
data2 = [];
while ~isempty(useryear)
    index = year == useryear(1) & period == userperiod(1);
    x = data(index,:);
    nx = size(x,2);
    xagg = nan(1,nx);
    for i = 1 : nx
        tmp = x(:,i);
        if options.ignorenan
            tmp = tmp(~isnan(tmp));
        end
        if isempty(tmp)
            xagg(1,i) = NaN;
        else
            try
                xagg(1,i) = options.method(tmp,1);
            catch %#ok<CTCH>
                try %#ok<TRYNC>
                    xagg(1,i) = options.method(tmp);
                end
            end
        end
    end
    data2 = [data2;xagg]; %#ok<AGROW>
    year(index) = [];
    period(index) = [];
    data(index,:) = [];
    index = useryear == useryear(1) & userperiod == userperiod(1);
    useryear(index) = [];
    userperiod(index) = [];
end

data2 = reshape(data2,[size(data2,1),tmpsize(2:end)]);

this.start = start2;
this.data = data2;
this = mytrim(this);
end
% $ dailyaggregate_().

% $ ***********************************************************************
function x = interp_(x,range1,freq1,freq2,options)
if isnan(x.start) && isempty(x.data)
    return
end
if isempty(range1)
    x = empty(x);
    return
end
if ~any(isinf(range1))
    range1 = range1(1) : range1(end);
end

[xdata,range1] = mygetdata(x,range1);
xsize = size(xdata);
xdata = xdata(:,:);

[startyear1,startper1] = dat2ypf(range1(1));
[endyear1,endper1] = dat2ypf(range1(end));

startyear2 = startyear1;
endyear2 = endyear1;
% Find the earliest freq2 period contained (at least partially) in freq1
% start period.
startper2 = 1 + floor((startper1-1)*freq2/freq1);
% Find the latest freq2 period contained (at least partially) in freq1 end
% period.
endper2 = ceil((endper1)*freq2/freq1);
range2 = ...
    datcode(freq2,startyear2,startper2) : datcode(freq2,endyear2,endper2);

grid1 = dat2grid(range1,options.position);
grid2 = dat2grid(range2,options.position);
xdata2 = interp1(grid1,xdata,grid2,options.method,'extrap');
if size(xdata2,1) == 1 && size(xdata2,2) == length(range2)
    xdata2 = xdata2(:);
else
    xdata2 = reshape(xdata2,[size(xdata2,1),xsize(2:end)]);
end
x.start = range2(1);
x.data = xdata2;
x = mytrim(x);
end
% $ interp_().

% $ ***********************************************************************
function x = interpmatch_(x,range1,freq1,freq2,options)

if isnan(x.start) && isempty(x.data)
    return
end
if isempty(range1)
    x = empty(x);
    return
end
if ~any(isinf(range1))
    range1 = range1(1) : range1(end);
end

n = freq2/freq1;
if n ~= round(n)
    error('iris:tseris',...
        'Source and target frequencies incompatible for ''%s'' interpolation.',...
        options.method);
end

[xdata,range1] = mygetdata(x,range1);
xsize = size(xdata);
xdata = xdata(:,:);

[startyear1,startper1] = dat2ypf(range1(1));
[endyear1,endper1] = dat2ypf(range1(end));

startyear2 = startyear1;
endyear2 = endyear1;
% Find the earliest freq2 period contained (at least partially) in freq1
% start period.
startper2 = 1 + floor((startper1-1)*freq2/freq1);
% Find the latest freq2 period contained (at least partially) in freq1 end
% period.
endper2 = ceil((endper1)*freq2/freq1);
range2 = ...
    datcode(freq2,startyear2,startper2) : datcode(freq2,endyear2,endper2);

[xdata2,flag] = interpmatchevaluate_(xdata,n);
if ~flag
    warning('iris:tseries',...
        'Cannot compute ''%s'' interpolation for series with within-sample NaNs.',...
        options.method);
end
if strcmpi(options.method,'quadavg')
    xdata2 = xdata2*n;
end

xdata2 = reshape(xdata2,[size(xdata2,1),xsize(2:end)]);
x.start = range2(1);
x.data = xdata2;
x = mytrim(x);

end

% $ ***********************************************************************
function [y2,flag] = interpmatchevaluate_(y1,n)
[nobs,ny] = size(y1);
y2 = nan(nobs*n,ny);

t1 = (1 : n)';
t2 = (n+1 : 2*n)';
t3 = (2*n+1 : 3*n)';
M = [...
    n, sum(t1), sum(t1.^2);...
    n, sum(t2), sum(t2.^2);...
    n, sum(t3), sum(t3.^2);...
    ];

flag = true;
for i = 1 : ny
    y1i = y1(:,i);
    [samplei,flagi] = getsample(y1i');
    flag = flag && flagi;
    if ~any(samplei)
        continue
    end
    y1i = y1i(samplei);
    nobsi = numel(y1i);
    yy = [ y1i(1:end-2), y1i(2:end-1), y1i(3:end) ]';
    b = nan(3,nobsi);
    b(:,2:end-1) = M \ yy;
    y2i = nan(n,nobsi);
    for t = 2 : nobsi-1
        y2i(:,t) = b(1,t)*ones(n,1) + b(2,t)*t2 + b(3,t)*t2.^2;
    end
    y2i(:,1) = b(1,2) + b(2,2)*t1 + b(3,2)*t1.^2;
    y2i(:,end) = b(1,end-1) + b(2,end-1)*t3 + b(3,end-1)*t3.^2;
    samplei = samplei(ones(1,n),:);
    samplei = samplei(:);
    y2(samplei,i) = y2i(:);
end
end
% $ interpmatchevaluate_().