%
%                    ROM_code.m
%
%

  clear all

%
%  Set limits used to approximate forced heat equation
%

  alpha = 2.0;
  L = 2;
  J = 40;
  K = 5;
  h = L/J;
  xj = 0:h:L;
  Nq = 41;
  hq = h/(Nq-1); 
  tf = 3;
  dt = 0.1;
  tvals = 0:dt:tf;
  xvals = xj;
  
%
% Create the right hand side vector
%

  xm1 = 0;
  xp1 = L;

  w = hq*ones(Nq,1); 
  w(1)  = 0.5*hq; 
  w(Nq) = 0.5*hq; 

  for j=1:J-1
    xm1 = xj(j+1)-h;
    xp1 = xj(j+1)+h;
    xm = xm1:hq:xj(j+1);
    xp = xj(j+1):hq:xp1;
    Gm = xm.*((xm-2).^4).*(xm-xm1)/h;
    Gp = xp.*((xp-2).^4).*(xp1-xp)/h;
    Im = Gm*w;
    Ip = Gp*w;
    g_vec(j,1) = Im+Ip;

    initm = sin(3*pi*xm/2).*(xm-xm1)/h;
    initp = sin(3*pi*xp/2).*(xp1-xp)/h;
    Initm = initm*w;
    Initp = initp*w;
    init_vec(j,1) = Initm + Initp;
  end

%
% Create the mass matrix M_mat, stiffness matrix K_mat, system matrix A_mat and vector G_vec
%

  M_mat = h*(diag((2/3)*ones(J-1,1)) + diag((1/6)*ones(J-2,1),-1) + diag((1/6)*ones(J-2,1),1));
  Minv = inv(M_mat);
  K_mat = (1/h)*(diag(2*ones(J-1,1)) + diag((-1)*ones(J-2,1),-1) + diag((-1)*ones(J-2,1),1));

  init_proj = Minv*init_vec;

%
% Create the vector of initial conditions
%

  r_init = sin(3*pi*xj(2:J)/2);
  r_init_full = sin(3*pi*xj/2);

  figure(8) 
  plot(xj(2:J),init_proj,xj(2:J),r_init)

%
% Numerically integrate the ODE system 
%

  A_mat = -alpha*Minv*K_mat;
  G_vec = Minv*g_vec;
  ode_options = odeset('RelTol',1e-6);
  [t,uvals] = ode45(@ODE_rhs,tvals,r_init,ode_options,A_mat,G_vec);
  M = length(tvals);
  Uvals = zeros(M,J+1);
  Uvals(:,2:J) = uvals;

%
% Plot the solution
%  

  figure(1)
  plot(xvals,Uvals(M,:))

  figure(2)
  mesh(Uvals)

%
% Consruct the covariance matrix and solve the eigenvalue problem.
%

  A = Uvals';
  K = (1/M)*A'*A;
  [V,D] = eig(K);

  r = rank(D);
  lambda = diag(D(1:r,1:r));
  V_red = V(1:5,:);
  lambda_tot = sum(lambda);
  for j=1:r
    ratio(j) = sum(lambda(1:j))/lambda_tot;
  end 
  Jr = 3;

  for j=1:Jr
    Phi(:,j) = (1/sqrt(M*lambda(j)))*A*V(:,j);
  end

%
% Plot the basis functions
%

  figure(3)
  plot(xvals,Phi(:,1),'-k',xvals,Phi(:,2),'--b',xvals,Phi(:,3),'-.r',xvals,0*xvals,'k-','linewidth',3)
  set(gca,'Fontsize',[20]);
  legend(' \phi_1^R(x)', ' \phi_2^R(x)',' \phi_3^R(x)','Location','NorthEast')
  xlabel('x-axis')
  ylabel('Reduced-Basis Functions')

%
% Create the derivative of the basis functions
%

  for i=1:J
    for j=1:Jr
      Phi_der(i,j) = (1/h)*(Phi(i+1,j) - Phi(i,j));
    end
  end

%
% Create the mass matrix M_red, stiffness matrix K_red rhs, and initial condition for the reduced-basis functions
%

  fvec = xvals.*(2-xvals).^4; 
  for i=1:Jr
    for j=1:Jr
      M_red(i,j) = h*sum(Phi(:,i).*Phi(:,j));
      K_red(i,j) = h*sum(Phi_der(:,i).*Phi_der(:,j));
    end
    gvec_red(i,1) = h*sum(Phi(:,i).*fvec');
  end
  
  M_red_inv = inv(M_red);
  for i=1:Jr
    init_vals(i,1) = h*sum(Phi(:,i).*r_init_full'); 
  end
  init_red = M_red_inv*init_vals;

  init_proj = init_red(1)*Phi(:,1) + init_red(2)*Phi(:,2) + init_red(3)*Phi(:,3);
  figure(4)
  plot(xvals,r_init_full,'-b',xvals,init_proj,'--r')

  A_red = -alpha*M_red_inv*K_red;
  G_red = M_red_inv*gvec_red;
  ode_options = odeset('RelTol',1e-6);
  [t,uvals_red] = ode45(@ODE_rhs,tvals,init_red,ode_options,A_red,G_red);

  Uvals_red = uvals_red(M,1)*Phi(:,1) + uvals_red(M,2)*Phi(:,2) + uvals_red(M,3)*Phi(:,3);

  figure(5)
  plot(xvals,Uvals(M,:),'-k',xvals,Uvals_red,'--b','linewidth',3)
  set(gca,'Fontsize',[20]);
  %legend(' u^J(T,x)', ' u^{{J_R}}(T,x)','Location','NorthEast')
  legend(' Full-Order','Reduced-Order','Location','NorthEast')
  xlabel('x-axis')
  ylabel('Solutions at Time T=3')


%
% Compute MC samples for a direct evaluation
%

  pt1 = 1.0;
  pt2 = 2.0;
  R_MC = 500;
  qr_MC = pt1 + (pt2-pt1)*rand(R_MC,1);
  wr_MC = (1/R_MC)*ones(size(qr_MC));

%
% Numerically integrate the full-order ODE system using random points for MC
%

  tic
  rvals_MC = zeros(R_MC,J+1);
  for r = 1:R_MC
    A_mat = -qr_MC(r)*Minv*K_mat;
    G_vec = Minv*g_vec;
    ode_options = odeset('RelTol',1e-6);
    [t,Rvals_MC] = ode45(@ODE_rhs,tvals,r_init,ode_options,A_mat,G_vec);
    rvals_MC(r,2:J) = Rvals_MC(length(t),:);
  end
  toc

%
% Numerically integrate the reduced-order ODE system using random points for MC
%

  tic
  for r = 1:R_MC
    A_red = -qr_MC(r)*M_red_inv*K_red;
    G_red = M_red_inv*gvec_red;
    ode_options = odeset('RelTol',1e-6);
    [t,uvals_red] = ode45(@ODE_rhs,tvals,init_red,ode_options,A_red,G_red);
    Uvals_red = uvals_red(M,1)*Phi(:,1) + uvals_red(M,2)*Phi(:,2) + uvals_red(M,3)*Phi(:,3);
    rvals_MC_red(r,:) = Uvals_red';
  end
  toc

%
% Compute the mean and variance using MC for full-order solution
%

  rmean_MC = wr_MC'*rvals_MC;
  var_MC = zeros(1,J+1);
  for r=1:R_MC
    var_MC = var_MC + wr_MC(r)*((rvals_MC(r,:) - rmean_MC).^2);
  end

%
% Uvals_red = uvals_red(M,1)*Phi(:,1) + uvals_red(M,2)*Phi(:,2) + uvals_red(M,3)*Phi(:,3);
%

  rmean_MC_red = wr_MC'*rvals_MC_red;
  var_MC_red = zeros(1,J+1);
  for r=1:R_MC
    var_MC_red = var_MC_red + wr_MC(r)*((rvals_MC_red(r,:) - rmean_MC_red).^2);
  end


%
% Plot the solution
%

  figure(6)
  plot(xj,var_MC,'k-',xj,var_MC_red,'b--','linewidth',3)
  set(gca,'Fontsize',[22]);
  legend(' Full-Order MC',' Reduced-Order MC','Location','Northeast')
  xlabel('x-axis')
  ylabel('Variance')