## 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>

## Singular value decomposition of the outputs of an ODE
# We construct a snapshot matrix from a 2nd order ODE
# (damped pendulum) and compute its SVd decomposition

## Dynamical system and snapshot generation
# We define the equations of motion of the angle repect to the vertical of a
# pendulum of unit length. Units are chosen such that the acceleration
# due to gravity is ten.

sys = @(x,t) [x(2); -10*sin(x(1))-0.8*x(2)];

# The parameters we will use are the initial angle and speed
nx      = 15;
x0      = linspace (-pi, pi, nx) * 2/3;
nv      = 10;
v0      = linspace (-pi, pi, nv);
[X0 V0] = meshgrid (x0, v0);
Q0      = [X0(:) V0(:)];

# Calculate sizes and create the time vector
ns  = size (Q0, 1);              # number of snapshots
nT  = 100;                       # number of time samples
T   = 10;                        # duration of simulation
t   = linspace (0, T, nT).';

if !exist ('Y','var')
  Y   = zeros (nT, ns);
  Q   = zeros (nT, ns, 2);
  printf ("Generating snapshots ...\n");
  fflush (stdout);
  tic;
  for i=1:ns
    tmp      = lsode (sys, Q0(i,:), t);
    Y(:,i)   = tmp(:,1);
    Q(:,i,:) = tmp;
  endfor
  toc;
endif

## POD on snapshot matrix
tol         = 0.2;
[phi lam B] = pod (Y, tol);
N           = size (phi, 2);
err         = norm (phi * B - Y) / norm (Y);
printf ("Size of basis (error: %.2f%%): %d\n", err*100, N);
fflush (stdout);
figure (1)
plot (t, phi);
axis tight
xlabel ("Time")
ylabel ("Basis")

## Surfaces to interpolate
BB = zeros (nv, nx, N);

if N > 3
  m = 2;
  k = ceil (N / m);
else
  m = 1;
  k = N;
end

figure (2)
clf
colormap (cubehelix);
for i=1:N
  subplot (m,k,i);
  BB(:,:,i) = reshape (B(i,:), nv, nx);
  surf (X0, V0, BB(:,:,i));
  axis tight
  xlabel ("angle");
  ylabel ("speed")
  zlabel (sprintf ("Coeff %d", i));
  view (40, 30);
  shading interp
endfor

## Try a single prediction
xp = min(x0) + ( max (x0) - min (x0) ) * rand ();
vp = min(v0) + ( max (v0) - min (v0) ) * rand ();
bp = zeros (N,1);
printf ("Emulation via interpolation ...\n");
fflush (stdout);
tic
for i=1:N
  bp(i) = interp2 (X0, V0, BB(:,:,i), xp, vp);
endfor
toc
yp  = phi * bp;
printf ("Simulation ...\n");
fflush (stdout);
tic
yp_ = lsode (sys, [xp vp], t);
toc
figure (3)
clf
plot (t, yp,'bo;emulated;', t, yp_(:,1), 'r-;simulated;')
axis tight
xlabel ("Time")
ylabel ("Angle")
