## Copyright (C) 2017 - Juan Pablo Carbajal
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program. If not, see <http://www.gnu.org/licenses/>.

## Author: Juan Pablo Carbajal <ajuanpi+dev@gmail.com>

## -*- texinfo -*-
## @deftypefn {} {[@var{P}, @var{L}] =} pod (@var{S}, @var{n})
## @deftypefnx {} {[@var{P}, @var{L}, @var{B}] =} pod (@var{S}, @var{n})
## @deftypefnx {} {[@dots{}] =} pod (@dots{}, @var{prop}, @var{value})
## Proper orthogonal decomposition basis and eigenvalues.
##
## Computes the POD basis @var{P} and the corresponding eigenvalues @var{L} of
## the M-by-N array @var{S}.
## The input argument @var{n} defines the number of basis elements
## (@var{n} <= @code{rank (@var{S})}) to be returned, either by specifying this
## number directly (in which case @var{n} should be an integer >=1)
## or by providing a relative reconstruction error in the Frobenius norm
## (in which case @var{n} should be in the interval (0,1)).
## The latter means that
##
## @code{norm (@var{P} * @var{B} - @var{S}) / norm (@var{S}) <= @var{n}}
##
## for some mixing coefficients @var{B}.
## The result is such that @code{@var{P.'} @var{P} == eye(size(W))}.
##
## The function accepts an inner product matrix W and a quadrture weights vector d
## such that
##
## @code{norm (@var{P} * @var{B} - @var{S}sqrtm (D)) / <= @var{n} norm (@var{S} sqrtm (D))}
##
## with @code{@var{P.'} W @var{P} == eye(size(W))} and @code{D = diag(d)}.
##
## @strong{Options}:
##
## @table @asis
## @item 'method'
## Defines the method to calculate the decompositon, valid methods are @asis{'svd'},
## @asis{'eig'}, and @asis{'eigs'}. If left empty (the default) a choice is made
## based on the @asis{'innerprod'} option described next.
##
## @item 'quadweights'
## Provide diagonal of D.
## It should be followed by a real positive vecotr. When given
## (and non-empty) it provides the weights of the quadrature formula used to
## link the continous POD the the discrete one.
##
## @item 'innerprod'
## Provide W matrix.
## It should be followed by a real symmetric positive definite matrix. When given
## (and non-empty) it provides the weighting inner product for the decomposition.
## In this case @asis{'method'} defaults to @asis{'eig'} as it is faster.
## @end table
##
## References:
## @indentedblock
## [1] S. Volkwein (2013). Proper Orthogonal Decomposition: Theory and Reduced-Order Modelling. Lecture notes.
##
## [2] Quarteroni, A., Manzoni, A., Negri, F. (2015). Reduced Basis Methods for Partial Differential Equations: An Introduction. Springer International Publishing.
## @end indentedblock
## @end deftypefn

function [P L varargout] = pod (S, l, varargin)

  if (nargin < 2)
    print_usage()
  endif

  ##### Parse options
  parser = inputParser ();
  parser.addParamValue ('innerprod', [], @isinner);
  parser.addParamValue ('method',  [], @isvmethod);
  parser.addParamValue ('quadweights',  [], @isquadw);
  
  parser.parse(varargin{:});
  res    = parser.Results;
  clear parser

  # Whether we need to calculate the weighted inner product version
  if !isempty (res.innerprod)
    if size(res.innerprod, 2) != size(S, 1)
      error ('pod:invalid-input-arg', "Wrong dimension of inner product\n")
    endif
    W = res.innerprod;
    is_weighted = true;
  else
    is_weighted = false;
  endif

  # Whether we need to calculate the weighting by quadrature weigts
  if !isempty (res.quadweights)
    if length(res.quadweights) != size(S, 2)
      error ('pod:invalid-input-arg', "Wrong number of quadrature weights\n")
    endif
    S = S * diag (sqrt (res.quadweights));
  endif
 
  # Whether the method was given
  if isempty (res.method)
    if is_weighted
      method = 'eig'; # The method of snapshots, or using eig is faster and more
                      # precise for the weighted case. See Ref[1] remark 1.3.3
    else
      method = 'svd'; # SVD is more robust for the unweighted case.
    endif
  else
    method = res.method;
  endif
  ##### End parse options

  switch method
    case 'svd'
      if is_weighted
        Wh      = real (sqrtm (W)); # if W is inner prod, the real here is
                                    # unnecessary, but we use it to avoid
                                    # complex eigs with very small imaginary part
        S       = Wh * S;
      endif
      [P s U] = svd (S, 1);
      L       = diag (s).^2 / (size (S, 1) - 1);

      if is_weighted
        P = Wh \ P;
      endif

    case {'eig','eigs'}
      if is_weighted
        K = S.' * W * S;
      else
        K = S.' * S;
      endif

      if isdimension (l) && l <= ceil ( 0.1 * min (size (S)) )
        [P L]  = eigs (K, l, 'lm');
      else
        [P L]  = eig (K);
      endif

      [L, o] = sort (sqrt (diag (L)), 'descend');
      P      = S * P(:,o) ./ L.';
  endswitch

  ## Reduce the basis according to l
  if !isdimension (l)
    if (l > 0 && l < 1) # l is reconstruction error
      ener = cumsum (L) ./ sum (L);
      l    = find (ener >= (1 - l^2), 1, 'first');
    else
      error ('Octave:invalid-input-arg',"invalid error parameter: %g\n", l);
    endif
  endif
  l = min (size (P,2), l);
  P = P(:, 1:l);
  L = L(1:l).';

  if nargout == 3
    # Compute mixing matrix
    switch method
      case 'svd'
        varargout{1} = s(1:l,1:l) * U(:,1:l).';
      case {'eig','eigs'}
        varargout{1} = P \ S;
    endswitch
  endif

endfunction

## Input validation functions
function tf = isquadw (x)
  tf = isreal (x) && (x > 0) && isvector (x);
  tf = tf || isempty (x);
endfunction

function tf = isinner (x)
  tf = isreal (x) && isdefinite (x);
  tf = tf || isempty (x);
endfunction

function tf = isvmethod (x)
  tf = ismember (x, {'svd', 'eig', 'eigs'});
  tf = tf || isempty (x);
endfunction

function tf = isdimension (x)
  tf = (x >= 1) && ( fix (x) - x == 0 );
endfunction

### TESTS
%!shared S
%! S = randn (100,10);

# invalid function calls
%!error(pod (S,-1));
%!error(pod (S,1,'innerprod', randn(100)));
%!error(pod (S,1,'method', 'other'));
%!error(pod (S,1,'innerprod', randn(100), 'quadweights', rand(1,5)));
%!error(pod (S,1,'quadweights', randn(1, 5)));

# weighted with identity same as unweighted
%!test
%! Pw = pod (S, 10, 'innerprod', eye(100));
%! P  = pod (S, 10);
%! s  = sign (diag (P.'*Pw)).';
%! assert (P,Pw.*s,sqrt(eps));
%!test
%! Pw = pod (S, 0.5, 'innerprod', eye(100));
%! P  = pod (S, 0.5);
%! s  = sign (diag (P.'*Pw)).';
%! assert (P,Pw.*s,sqrt(eps));

# eig, eigs, svd give same result
%!test
%! Pe = pod (S, 1, 'method', 'eig'); # uses eigs
%! P  = pod (S, 1, 'method', 'svd');
%! s  = sign (diag (P.'*Pe)).';
%! assert (P,Pe.*s,sqrt(eps));
%!test
%! Pe = pod (S, 1, 'method', 'eig');    # uses eigs
%! P  = pod (S, 0.99, 'method', 'eig'); # uses eig
%! s  = sign (diag (P.'*Pe)).';
%! assert (P,Pe.*s,sqrt(eps));

# tolerance works
%!test
%! tol = 1e-3;
%! [P,l,B] = pod (S, tol);
%! assert (norm (P*B - S) <= tol * norm (S));

%!demo
%! N     = 100;
%! t     = linspace (-1, 1, N).';
%! P_    = t.^[0:5];
%! M     = 50;
%! B_    = randn (size (P_, 2), M);
%! S     = P_ * B_;
%! [P,L,B] = pod (S, 6);
%! subplot(2,1,1)
%! plot (t, P);
%! axis tight
%! printf ("Error: %f\n", norm (P * B - S) / norm (S));
%! tol = 0.4;
%! [P,L,B] = pod (S, tol);
%! subplot(2,1,2)
%! plot (t, P);
%! axis tight
%! printf ("Error (given): %f (%f)\n", norm (P * B - S) / norm (S), tol);
