% LISOVFIT uses LISO to fit a pole/zero model to the input frequency-series.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% DESCRIPTION: LISOVFIT uses LISO to fit a pole/zero model to the input
%              frequency-series.
%
% CALL:        >> pzm = lisovfit(a,pl)
%
% INPUTS:      pl   - a parameter list
%              a    - input analysis object
%
% OUTPUTS:
%              pzm  - the fitted pzmodel.
%
% <a href="matlab:utils.helper.displayMethodInfo('ao', 'lisovfit')">Parameters Description</a>
%
% VERSION:     $Id: lisovfit.m,v 1.14 2011/04/08 08:56:14 hewitson Exp $
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function varargout = lisovfit(varargin)

  % Check if this is a call for parameters
  if utils.helper.isinfocall(varargin{:})
    varargout{1} = getInfo(varargin{3});
    return
  end

  import utils.const.*
  utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename);

  % Collect input variable names
  in_names = cell(size(varargin));
  for ii = 1:nargin,in_names{ii} = inputname(ii);end

  % Collect all AOs and plists
  [bs, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
  [pl, pl_invars] = utils.helper.collect_objects(varargin(:), 'plist', in_names);

  % combine plists
  pl = parse(pl, getDefaultPlist());

  % Extract parameters
  pzm0  = find(pl, 'PZM0');
  pzml  = find(pl, 'PZML');
  pzmu  = find(pl, 'PZMU');
  delay = find(pl, 'delay');
  f1    = find(pl, 'f1');
  f2    = find(pl, 'f2');
  nf    = find(pl, 'nf');
  method    = find(pl, 'method');
  np1    = find(pl, 'np1');
  np2    = find(pl, 'np2');

  
  % Check inputs
  if isempty(f1) || isempty(f2) || isempty(nf)
    error('### Specify the full frequency range with f1, f2, and nf');
  end
  if f1 > f2
    error('### The starting frequency should be less than the end frequency');
  end
  if numel(delay) ~= 3
    error('### The delay must be specified as a 3 element numerical vector: [limit start upper]');
  end
  if ~(delay(1) < delay(2))
    error('### The lower limit for the delay must be less than the starting guess.');
  end
  if ~(delay(2) < delay(3))
    error('### The upper limit for the delay must be greater than the starting guess.');
  end
  if numel(pzm0.poles) ~= numel(pzml.poles)
    error('### The starting, lower, and upper models must have the same number of poles');
  end
  if numel(pzm0.poles) ~= numel(pzmu.poles)
    error('### The starting, lower, and upper models must have the same number of poles');
  end
  if numel(pzm0.zeros) ~= numel(pzml.zeros)
    error('### The starting, lower, and upper models must have the same number of zeros');
  end
  if numel(pzm0.zeros) ~= numel(pzmu.zeros)
    error('### The starting, lower, and upper models must have the same number of zeros');
  end

  % Loop over input AOs
  pzms = [];
  pzmls = [];
  pzmus = [];
  for j=1:numel(bs)
    if isa(bs(j).data, 'fsdata')
      % Generate temp file names
      outfile  = [tempname '.fil'];
      datafile = [tempname '.dat'];
      % Write LISO fit file
      switch lower(method)
        case 'fit'
          if isreal(bs(j).data.getY)
            writeLISOfitFile(pzm0, pzml, pzmu, delay, f1, f2, nf, outfile, datafile, false);
          else
            writeLISOfitFile(pzm0, pzml, pzmu, delay, f1, f2, nf, outfile, datafile, true);
          end
        case 'vfit'
          writeLISOvfitFile(np1,np2,outfile, datafile);
        otherwise
          error('### method should be ''vfit'' or ''fit''.');
      end
      % Export AO data file
      export(bs(j), datafile);
      % call fil
      utils.bin.fil(outfile);
      % get fitted model
      pzmfit = pzmodel(outfile);
      % pzmodel returns 3 models, set name for the first one
      pzmfit(1).name = sprintf('fit(%s)', pzm0.name);
      % Set units
      [ounits, iunits] = factor(bs(j).data.yunits);
      pzmfit(1).setIunits(iunits);
      pzmfit(1).setOunits(ounits);
      pzmfit(2).setIunits(iunits);
      pzmfit(2).setOunits(ounits);
      pzmfit(3).setIunits(iunits);
      pzmfit(3).setOunits(ounits);
      % add history
      pzmfit(1).addHistory(getInfo('None'), pl, ao_invars(j), bs(j).hist);
      pzmfit(2).addHistory(getInfo('None'), pl, ao_invars(j), bs(j).hist);
      pzmfit(3).addHistory(getInfo('None'), pl, ao_invars(j), bs(j).hist);
      % add to output
      pzms = [pzms pzmfit(1)];
      pzmls = [pzmls pzmfit(2)];
      pzmus = [pzmus pzmfit(3)];
    else
      error('### unknown data type.');
    end
  end

  % Set outputs
  if nargout == 1
    varargout{1} = pzms;
  elseif nargout == 3
    varargout{1} = pzms;
    varargout{2} = pzmls;
    varargout{3} = pzmus;
  else
    error('### Incorrect number of output arguments.');
  end
end

%----------------------------------------------------------
% Write a liso vector fitting  file
%
function writeLISOvfitFile(np1, np2, outfile, datafile)
  fd = fopen(outfile, 'w+');

  fprintf(fd, '# Temporary LISO fitting file \n');

  % Write fit command
   fprintf(fd, 'vfit %s reim rel %d %d \n', datafile,np1,np2);
   fprintf(fd, '\n');

  % Close file
  fclose(fd);
end

%----------------------------------------------------------
% Write a liso fit file
%
function writeLISOfitFile(pzm0, pzml, pzmu, delay, f1, f2, nf, outfile, datafile, complexData)
  fd = fopen(outfile, 'w+');

  fprintf(fd, '# Temporary LISO fitting file from pzmodel: %s\n\n', pzm0.name);

  % first output poles
  Np = numel(pzm0.poles);
  for k=1:Np
    pole = pzm0.poles(k);
    lpole = pzml.poles(k);
    upole = pzmu.poles(k);
    fprintf(fd, '# POLE %d\n', k);
    % write pole start
    if isnan(pole.q)
      if ~isnan(lpole.q) || ~isnan(upole.q)
        error('### Poles in start, lower, and upper models must be of the same for (real, complex)');
      end
      fprintf(fd, 'pole %g\n',  pole.f);
    else
      fprintf(fd, 'pole %g %g\n',  pole.f, pole.q);
    end
    % write param range
    %   param pole2:f 1.0746e-06 0.00010746
    fprintf(fd, 'param pole%d:f %g %g\n', k-1, lpole.f, upole.f);
    if ~isnan(pole.q)
      fprintf(fd, 'param pole%d:q %g %g\n', k-1, lpole.q, upole.q);
    end
    fprintf(fd, '\n');
  end

  % then output zeros
  Nz = numel(pzm0.zeros);
  for k=1:Nz
    zero = pzm0.zeros(k);
    lzero = pzml.zeros(k);
    uzero = pzmu.zeros(k);
    fprintf(fd, '# ZERO %d\n', k);
    % write pole start
    if isnan(zero.q)
      if ~isnan(lzero.q) || ~isnan(uzero.q)
        error('### Zeros in start, lower, and upper models must be of the same for (real, complex)');
      end
      fprintf(fd, 'zero %g\n',  zero.f);
    else
      fprintf(fd, 'zero %g %g\n',  zero.f, zero.q);
    end
    % write param range
    %   param pole2:f 1.0746e-06 0.00010746
    fprintf(fd, 'param zero%d:f %g %g\n', k-1, lzero.f, uzero.f);
    if ~isnan(zero.q)
      fprintf(fd, 'param zero%d:q %g %g\n', k-1, lzero.q, uzero.q);
    end
    fprintf(fd, '\n');
  end

  % Write delay out
  if ~isempty(delay)
    fprintf(fd, 'delay %g\n', delay(2));
    fprintf(fd, 'param delay %g %g\n', delay(1), delay(3));
  end
  fprintf(fd, '\n');

  % Write factor
  fprintf(fd, 'factor %g\n', pzm0.gain);
  fprintf(fd, 'param factor %g %g\n', pzml.gain, pzmu.gain);
  fprintf(fd, '\n');

  % Write frequency command
  fprintf(fd, 'freq lin %g %g %g\n', f1, f2, nf);
  fprintf(fd, '\n');

  % Write fit command
  if complexData
    fprintf(fd, 'fit %s reim semi\n', datafile);
  else
    fprintf(fd, 'fit %s abs rel\n', datafile);
  end
  fprintf(fd, 'rewrite samebetter\n');
  fprintf(fd, '\n');

  % Close file
  fclose(fd);
end

%--------------------------------------------------------------------------
% Get Info Object
%--------------------------------------------------------------------------
function ii = getInfo(varargin)
  if nargin == 1 && strcmpi(varargin{1}, 'None')
    sets = {};
    pls  = [];
  else
    sets = {'Default'};
    pls  = getDefaultPlist;
  end
  % Build info object
  ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: lisovfit.m,v 1.14 2011/04/08 08:56:14 hewitson Exp $', sets, pls);
  ii.setModifier(false);
end

%--------------------------------------------------------------------------
% Get Default Plist
%--------------------------------------------------------------------------
function plout = getDefaultPlist()
  persistent pl;  
  if exist('pl', 'var')==0 || isempty(pl)
    pl = buildplist();
  end
  plout = pl;  
end

function pl = buildplist()
  
  pl = plist();
  
  % PZM0
  p = param({'PZM0', 'A pzmodel describing the starting guess.'}, {1, {pzmodel}, paramValue.OPTIONAL});
  pl.append(p);
  
  % PZML
  p = param({'PZMU', 'A pzmodel describing the upper-bound.'}, {1, {pzmodel}, paramValue.OPTIONAL});
  pl.append(p);
  
  % PZMU
  p = param({'PZML', 'A pzmodel describing the lower-bound.'}, {1, {pzmodel}, paramValue.OPTIONAL});
  pl.append(p);
  
  % Delay
  p = param({'Delay', 'A 3-element numeric vector describing<br>a time-delay to include in the fit.'}, {1, {[0 1 10]}, paramValue.OPTIONAL});
  pl.append(p);
  
  % F1
  p = param({'F1', 'A start freqeuency to fit over'}, {1, {0}, paramValue.OPTIONAL});
  pl.append(p);
  
  % F2
  p = param({'F2', 'The end freqeuency to fit over'}, {1, {1}, paramValue.OPTIONAL});
  pl.append(p);

  % Nf
  p = param({'NF','The number of frequencies to include in the fit.'}, {1, {100}, paramValue.OPTIONAL});
  pl.append(p);
  
  % method
  p = param({'method', 'The fitting method.'}, {1, {'fit', 'vfit'}, paramValue.SINGLE});
  pl.append(p);
  
  % NP1
  p = param({'NP1', 'The minimum number of poles for vfit.'}, {1, {2}, paramValue.OPTIONAL});
  pl.append(p);
  
  % NP2
  p = param({'NP2', 'The maximum number of poles for vfit.'}, {1, {10}, paramValue.OPTIONAL});
  pl.append(p);
  
  

end

% PARAMETERS:
%              'PZM0'  - a pzmodel describing the starting guess.
%              'PZML'  - a pzmodel describing the lower-bound.
%              'PZMU'  - a pzmodel describing the upper-bound.
%              'DELAY' - a 3-element numeric vector describing a time-delay
%                        to include in the fit: [lower start upper]
%              'f1'    - start freqeuency to fit over [default: taken from input AO]
%              'f2'    - end freqeuency to fit over [default: taken from input AO]
%              'nf'    - number of frequencies to fit [default: taken from input AO]
%              'method'- 'fit' for fitting or 'vfit' for vector fitting [default: fit]
%              'np1'    - lower num. poles for vfit [default: 2]
%              'np2'    - upper num. poles for vfit [default: 10]