function [Y,varargout] = x12(X,STARTDATE,DUMMY,OPT)
% x12  [Not a public function] Matlab interface for X12-ARIMA.
%
% Backend IRIS function.
% No help provided.

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

switch lower(OPT.mode)
    case {0,'mult','m'}
        OPT.mode = 'mult';
    case {1,'add','a'}
        OPT.mode = 'add';
    case {2,'pseudo','pseudoadd','p'}
        OPT.mode = 'pseudoadd';
    case {3,'log','logadd','l'}
        OPT.mode = 'logadd';
    otherwise
        OPT.mode = 'auto';
end

if ischar(OPT.output)
    OPT.output = {OPT.output};
elseif ~iscellstr(OPT.output)
    OPT.output = {'d11'};
end
noutput = length(OPT.output);

% Get the entire path to this file.
thisdir = fileparts(mfilename('fullpath'));

if strcmpi(OPT.specfile,'default')
    OPT.specfile = fullfile(thisdir,'default.spc');
end

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

kb = OPT.backcast;
kf = OPT.forecast;
nper = size(X,1);
nx = size(X,2);

% Output data.
[varargout{1:noutput}] = deal(nan(nper,nx));
% Output file(s).
varargout{noutput+1}(1:nx) = {''};
% Error file(s).
varargout{noutput+2}(1:nx) = {''};

freq = datfreq(STARTDATE);
if freq ~= 4 && freq ~= 12
    utils.warning('x12', ...
        'X12 runs only on quarterly or monthly data.');
    return
end

threeyearswarn = false;
seventyyearswarn = false;
bcast15yearswarn = false;
nanwarn = false;
mdl = struct('mode',NaN,'spec',[],'ar',[],'ma',[]);
mdl = mdl(ones(1,nx));
outp = cell(1,nx);
err = cell(1,nx);
Y = nan(nper+kb+kf,nx);
for i = 1 : nx
    [~,tmptitle] = fileparts(tempname(pwd()));
    first = find(~isnan(X(:,i)),1);
    last = find(~isnan(X(:,i)),1,'last');
    data = X(first:last,i);
    if length(data) < 3*freq
        threeyearswarn = true;
    elseif length(data) > 70*freq
        seventyyearswarn = true;
    elseif any(isnan(data))
        nanwarn = true;
    else
        if length(data) > 15*freq && kb > 0
            bcast15yearswarn = true;
        end
        offset = first - 1;
        [aux,fcast,bcast] = ...
            xxrunx12(thisdir,tmptitle,data,STARTDATE+offset,DUMMY,OPT);
        for j = 1 : noutput
            varargout{j}(first:last,i) = aux(:,j);
        end
        % Append forecasts and backcasts to original data.
        Y(first:last+kb+kf,i) = [bcast;data;fcast];
        % Catch output file.
        if exist([tmptitle,'.out'],'file')
            outp{i} = xxreadoutputfile([tmptitle,'.out']);
        end
        % Catch error file.
        if exist([tmptitle,'.err'],'file')
            err{i} = xxreadoutputfile([tmptitle,'.err']);
        end
        % Catch ARIMA model specification.
        if exist([tmptitle,'.mdl'],'file')
            mdl(i) = xxreadmodel([tmptitle,'.mdl'],outp{i});
        end
        % Delete all X12 files.
        if ~isempty(OPT.saveas)
            dosaveas();
        end
        if OPT.cleanup
            stat = warning();
            warning('off'); %#ok<WNOFF>
            delete([tmptitle,'.*']);
            warning(stat);
        end
    end
end

dowarnings();

varargout{noutput+1} = outp;
varargout{noutput+2} = err;
varargout{noutput+3} = mdl;

% Nested functions.

%**************************************************************************
    function dowarnings()
        if threeyearswarn
            utils.warning('x12', ...
                'X12 requires three or more years of observations.');
        end
        if seventyyearswarn
            utils.warning('x12', ...
                'X12 cannot handle more than 70 years of observations.');
        end
        if bcast15yearswarn
            utils.warning('x12', ...
                'X12 does not produce backcasts for time seris longer than 15 years.');
        end
        if nanwarn
            utils.warning('x12', ...
                'Input data contain within-sample NaNs.');
        end
    end
% dowarnings().

%**************************************************************************
    function dosaveas()
        [fpath,ftitle] = fileparts(OPT.saveas);
        list = dir([tmptitle,'.*']);
        for ii = 1 : length(list)
            [~,~,fext] = fileparts(list(ii).name);
            copyfile(list(ii).name,fullfile(fpath,[ftitle,fext]));
        end
    end
% dosaveas().

end

% Subfunctions.

%**************************************************************************
function [DATA,FCAST,BCAST] = ...
    xxrunx12(THISDIR,TMPTITLE,DATA,STARTDATE,DUMMY,OPT)
% Flip sign if all values are negative
% so that multiplicative mode is possible.
flipsign = false;
if all(DATA < 0)
    DATA = -DATA;
    flipsign = true;
end

nonpositive = any(DATA <= 0);
if strcmp(OPT.mode,'auto')
    if nonpositive
        OPT.mode = 'add';
    else
        OPT.mode = 'mult';
    end
elseif strcmp(OPT.mode,'mult') && nonpositive
    utils.warning('x12', ...
        ['Unable to use multiplicative mode because of ', ...
        'input data combine positive and non-positive numbers. ', ...
        'Switching to additive mode.']);
    OPT.mode = 'add';
end

xxspecfile(TMPTITLE,DATA,STARTDATE,DUMMY,OPT);

% Set up a system command to run X12a.exe, enclosing the command in
% double quotes.
if ispc()
    command = ['"',fullfile(THISDIR,'x12a.exe'),'" ',TMPTITLE];
elseif isunix()
    command = ['"',fullfile(THISDIR,'x12a'),'" ',TMPTITLE];
else
    command = ['wine "',fullfile(THISDIR,'x12a.exe'),'" ',TMPTITLE];
end

[status,result] = system(command);
if OPT.display
    disp(result);
end

% Return NaNs if X12 fails.
if status ~= 0
    DATA(:) = NaN;
    utils.warning('x12','Unable to execute X12.');
    return
end

% Read in-sample results.
nper = length(DATA);
[DATA,flagdata] = xxgetoutputdata(TMPTITLE,nper,OPT.output,2);

% Read forecasts.
FCAST = zeros(0,1);
flagfcast = true;
kf = OPT.forecast;
if kf > 0
    [FCAST,flagfcast] = xxgetoutputdata(TMPTITLE,kf,{'fct'},4);
end

% Read backcasts.
BCAST = zeros(0,1);
kb = OPT.backcast;
if kb > 0
    BCAST = xxgetoutputdata(TMPTITLE,kb,{'bct'},4);
end

if ~flagdata || ~flagfcast
    utils.warning('x12', ...
        ['Unable to read X12 output file(s). This is most likely because ', ...
        'X12 failed to estimate an appropriate seasonal model or failed to ', ...
        'converge. Run X12 with three output arguments ', ...
        'to capture the X12 output and error messages.']);
    return
end

if flipsign
    DATA = -DATA;
    FCAST = -FCAST;
    BCAST = -BCAST;
end

end
% xxrunx12().

%**************************************************************************
function xxspecfile(TMPTITLE,DATA,STARTDATE,DUMMY,OPT)
% xxspecfile  Create and save SPC file based on a template.

[datayear,dataper] = dat2ypf(STARTDATE);
[dummyyear,dummyper] = dat2ypf(STARTDATE-OPT.backcast);

spec = file2char(OPT.specfile);

% Time series specs.

% Check for the required placeholders $series_data$ and $x11_save$:
if length(strfind(spec,'$series_data$')) ~= 1 ...
        || length(strfind(spec,'$x11_save$')) ~= 1
    utils.error('x12', ...
        ['Invalid X12 spec file. Some of the required placeholders, ', ...
        '$series_data$ and $x11_save$, are missing or used more than once.']);
end

% Data.
spec = strrep(spec,'$series_data$',sprintf('    %.8f\r\n',DATA));
% Seasonal period.
spec = strrep(spec,'$series_freq$',sprintf('%g',datfreq(STARTDATE)));
% Start date.
spec = strrep(spec,'$series_startyear$',sprintf('%g',round(datayear)));
spec = strrep(spec,'$series_startper$',sprintf('%g',round(dataper)));

% Transform specs.
if any(strcmp(OPT.mode,{'mult','pseudoadd','logadd'}))
    replace = 'log';
else
    replace = 'none';
end
spec = strrep(spec,'$transform_function$',replace);

% AUTOMDL specs.
spec = strrep(spec,'$maxorder$',sprintf('%g %g',round(OPT.maxorder)));

% FORECAST specs.
spec = strrep(spec,'$forecast_maxback$',sprintf('%g',OPT.backcast));
spec = strrep(spec,'$forecast_maxlead$',sprintf('%g',OPT.forecast));

% REGRESSION specs. If there's no REGRESSSION variable, we cannot include
% the spec in the spec file because X12 would complain. In that case, we
% keep the entire spec commented out. If tdays is requested but no user
% dummies are specified, we need to keep the dummy section commented out,
% and vice versa.
if OPT.tdays || ~isempty(DUMMY)
    spec = strrep(spec,'#regression ','');
    if OPT.tdays
        spec = strrep(spec,'#tdays ','');
        spec = strrep(spec,'$tdays$','td');
    end
    if ~isempty(DUMMY)
        spec = strrep(spec,'#dummy ','');
        ndummy = size(DUMMY,2);
        name = '';
        for i = 1 : ndummy
            name = [name,sprintf('dummy%g ',i)]; %#ok<AGROW>
        end
        spec = strrep(spec,'$dummy_type$',lower(OPT.dummytype));
        spec = strrep(spec,'$dummy_name$',name);
        spec = strrep(spec,'$dummy_data$',sprintf('    %.8f\r\n',DUMMY));
        spec = strrep(spec,'$dummy_startyear$', ...
            sprintf('%g',round(dummyyear)));
        spec = strrep(spec,'$dummy_startper$', ...
            sprintf('%g',round(dummyper)));
    end
end

% ESTIMATION specs.
spec = strrep(spec,'$maxiter$',sprintf('%g',round(OPT.maxiter)));
spec = strrep(spec,'$tolerance$',sprintf('%e',OPT.tolerance));

% X11 specs.
spec = strrep(spec,'$x11_mode$',OPT.mode);
spec = strrep(spec,'$x11_save$',sprintf('%s ',OPT.output{:}));

char2file(spec,[TMPTITLE,'.spc']);

end
% xxspecfile().

%**************************************************************************
function [DATA,FLAG] = xxgetoutputdata(TMPTITLE,NPER,OUTPUT,NCOL)

if ischar(OUTPUT)
    OUTPUT = {OUTPUT};
end
FLAG = true;
DATA = zeros(NPER,0);
format = repmat(' %f',[1,NCOL]);
for ioutput = 1 : length(OUTPUT)
    output = OUTPUT{ioutput};
    fid = fopen(sprintf('%s.%s',TMPTITLE,output),'r');
    if fid > -1
        fgetl(fid); % skip first 2 lines
        fgetl(fid);
        read = fscanf(fid,format);
        fclose(fid);
    else
        read = [];
    end
    if length(read) == NCOL*NPER
        read = reshape(read,[NCOL,NPER]).';
        DATA(:,end+1) = read(:,2); %#ok<AGROW>
    else
        DATA(:,end+1) = NaN; %#ok<AGROW>
        FLAG = false;
    end
end
end
% xxgetoutput().

%**************************************************************************
function C = xxreadoutputfile(FNAME)
C = file2char(FNAME);
C = strfun.converteols(C);
C = strfun.removeltel(C);
C = regexprep(C,'\n\n+','\n\n');
end
% xxreadoutputfile().

%**************************************************************************
function MDL = xxreadmodel(FNAME,OUTPFILE)

C = file2char(FNAME);
MDL = struct('mode',NaN,'spec',[],'ar',[],'ma',[]);

% ARIMA spec block.
arima = regexp(C,'arima\s*\{\s*model\s*=([^\}]+)\}','once','tokens');
if isempty(arima) || isempty(arima{1})
    return
end
arima = arima{1};

% Non-seasonal and seasonal ARIMA model specification
spec = regexp(arima,'\((.*?)\)\s*\((.*?)\)','once','tokens');
if isempty(spec) || length(spec) ~= 2 ...
        || isempty(spec{1}) || isempty(spec{2})
    return
end
specar = sscanf(spec{1},'%g').';
specma = sscanf(spec{2},'%g').';
if isempty(specar) || isempty(specma)
    return
end

% Estimated AR and MA coefficients.
estar = regexp(arima,'ar\s*=\s*\((.*?)\)','once','tokens');
estma = regexp(arima,'ma\s*=\s*\((.*?)\)','once','tokens');
if isempty(estar) && isempty(estma)
    return
end
try
    estar = sscanf(estar{1},'%g').';
catch %#ok<CTCH>
    estar = [];
end
try
    estma = sscanf(estma{1},'%g').';
catch %#ok<CTCH>
    estma = [];
end
if isempty(estar) && isempty(estma)
    return
end

mode = NaN;
if ~isempty(OUTPFILE) && ischar(OUTPFILE)
    tok = regexp(OUTPFILE,'Type of run\s*-\s*([\w\-]+)','tokens','once');
    if ~isempty(tok) && ~isempty(tok{1})
        mode = tok{1};
    end
end

% Create output struct only after we make sure all pieces have been read in
% all right.
MDL.mode = mode;
MDL.spec = {specar,specma};
MDL.ar = estar;
MDL.ma = estma;

end
% xxreadmodel().