%% Experiment with preconditioning for the NTW problem


%% Setup
addpath('~/Dropbox/research/pmpack');

%% Construct a small problem

% these parameters make the problem small and amenable to interactive
% computaiton
d = 2;
p_order = 2; % we'll assume that 5th order expansion is okay.

P = ntw_cgikl(d); % create a 2d NTW problem

tol=1e-6;
maxiter=5000;

%% Solve the small system and plot convergence
S = instrumented_solver('minres','maxiter',maxiter,'tol',tol);
X = spectral_galerkin(P.A,P.b,P.s,p_order,'verbose',1,'lowmem',1,...
    'solver',S.solver);
assert(S.stats()==0,'linear system solve did not converge'); % check the solve
h = S.histguess();
semilogy(h(:,2),h(:,1),'.-'); xlabel('time (s)'); ylabel('relres');

%% Add a preconditioner
s0 = midpoint(P.s);
A0 = P.A(s0);
R = chol(A0);

pcon1=@(x) galerkin_preconditioner(R',x);
pcon2=@(x) galerkin_preconditioner(R,x);

S = instrumented_solver('minres','maxiter',maxiter,'tol',tol,...
        'precond',{pcon1,pcon2});
X = spectral_galerkin(P.A,P.b,P.s,p_order,'verbose',1,'lowmem',1,...
    'solver',S.solver);
assert(S.stats()==0,'linear system solve did not converge'); % check the solve
h = S.histguess();
semilogy(h(:,2),h(:,1),'.-'); xlabel('time (s)'); ylabel('relres');

%% Construct a real problem

% these paramters make the problem more realistic, but take a while to
% solve
d = 4;
p_order = 5;
pname = sprintf('cgikl-d-%i-p-%i',d,p_order);

P = ntw_cgikl(d); % create a 5d NTW problem

tol=1e-6;
maxiter=10000;

%% Check that the problem runs
X = spectral_galerkin(P.A,P.b,P.s,p_order,'verbose',1,'lowmem',1,...
    'solver',@(A,b) minres(A,b,1e-2,10),'parallel_rhs',1);
    
%% Compute our preconditioning study
% We want to evaluate the following points
% 1. No preconditioning
% 2. Random points in the space (25)
% 3. The midpoint of the space
% 4. The mean matrix

all_results = [];

%%
% 1. No preconditioning
fprintf('No preconditioning trial\n');
S = instrumented_solver('minres','maxiter',maxiter,'tol',tol);
t0 = tic;
X = spectral_galerkin(P.A,P.b,P.s,p_order,'verbose',1,'lowmem',1,...
    'solver',S.solver,'parallel_rhs',1);
dt = toc(t0);    
[flag resvec] = S.stats();


results = [];
results.X = X;
results.dt = dt;
results.flag = flag;
results.resvec = resvec;
results.hist = S.histguess();
results.matvectimes = S.matvectimes();

all_results.nopre = results;
save(sprintf('partial.mat'),'all_results');


%%
% 2. Random points
npts = 25;
pts = zeros(npts,d);
for i=1:npts
    for j=1:d
        pts(i,j) = P.s(j).l + (P.s(j).r-P.s(j).l)*rand;
    end
end

results = [];
for i=1:npts
    fprintf('Random point %i\n', i);
    s0 = pts(i,:);
    results(i).pt = s0;
    
    t0 = tic;
    A0 = P.A(s0);
    R = chol(A0);
    results(i).dt_precond = toc(t0);

    pcon1=@(x) galerkin_preconditioner(R',x);
    pcon2=@(x) galerkin_preconditioner(R,x);

    S = instrumented_solver('minres','maxiter',maxiter,'tol',tol,...
            'precond',{pcon1,pcon2});
    t0 = tic;
    X = spectral_galerkin(P.A,P.b,P.s,p_order,'verbose',1,'lowmem',1,...
        'solver',S.solver,'parallel_rhs',1);
    dt = toc(t0);
    [flag resvec] = S.stats();


    results(i).X = X;
    results(i).dt = dt;
    results(i).flag = flag;
    results(i).resvec = resvec;
    results(i).hist = S.histguess();
    results(i).matvectimes = S.matvectimes();        
    
    all_results.random = results;
    save(sprintf('partial.mat'),'all_results');
end


%% 
% 3. Midpoint
fprintf('Midpoint trial\n');
results = [];

s0 = midpoint(P.s);
results.pt = s0;

t0 = tic;
A0 = P.A(s0);
R = chol(A0);
results.dt_precond = toc(t0);

pcon1=@(x) galerkin_preconditioner(R',x);
pcon2=@(x) galerkin_preconditioner(R,x);

S = instrumented_solver('minres','maxiter',maxiter,'tol',tol,...
    'precond',{pcon1,pcon2});

t0 = tic;
X = spectral_galerkin(P.A,P.b,P.s,p_order,'verbose',1,'lowmem',1,...
    'solver',S.solver,'parallel_rhs',1);
dt = toc(t0);    

[flag resvec] = S.stats();
results.X = X;
results.dt = dt;
results.flag = flag;
results.resvec = resvec;
results.hist = S.histguess();
results.matvectimes = S.matvectimes();

all_results.midpoint = results;
save(sprintf('partial.mat'),'all_results');

%%
% 4. Mean matrix
norders = 3;
results = [];
for i=2:(2+norders)
    fprintf('Mean order %i\n', i);
    results(i).order = i;
    
    t0 = tic;
    A0 = mean_matrix(P.A,P.s,i);
    R = chol(A0);
    results(i).dt_precond = toc(t0);

    pcon1=@(x) galerkin_preconditioner(R',x);
    pcon2=@(x) galerkin_preconditioner(R,x);

    S = instrumented_solver('minres','maxiter',maxiter,'tol',tol,...
            'precond',{pcon1,pcon2});
    t0 = tic;
    X = spectral_galerkin(P.A,P.b,P.s,p_order,'verbose',1,'lowmem',1,...
        'solver',S.solver,'parallel_rhs',1);
    dt = toc(t0);
    [flag resvec] = S.stats();


    results(i).X = X;
    results(i).dt = dt;
    results(i).flag = flag;
    results(i).resvec = resvec;
    results(i).hist = S.histguess();
    results(i).matvectimes = S.matvectimes();        
    
    all_results.mean = results;
    save(sprintf('partial.mat'),'all_results');
end

%%
% 5. Max eigenvalue (fast approx)

fprintf('Max-eig trial (Gauss Fast)\n');
results = [];

t0 = tic;
[s0,A0] = max_eig_matrix(P.A,P.s,max(3,ceil(floor(p_order)/2)));
results.pt = s0;
R = chol(A0);
results.dt_precond = toc(t0);

pcon1=@(x) galerkin_preconditioner(R',x);
pcon2=@(x) galerkin_preconditioner(R,x);

S = instrumented_solver('minres','maxiter',maxiter,'tol',tol,...
    'precond',{pcon1,pcon2});

t0 = tic;
X = spectral_galerkin(P.A,P.b,P.s,p_order,'verbose',1,'lowmem',1,...
    'solver',S.solver,'parallel_rhs',1);
dt = toc(t0);    

[flag resvec] = S.stats();
results.X = X;
results.dt = dt;
results.flag = flag;
results.resvec = resvec;
results.hist = S.histguess();
results.matvectimes = S.matvectimes();

all_results.maxeig_fast = results;
save(sprintf('partial.mat'),'all_results');

%%
% 5. Max eigenvalue (exact for galerkin)

fprintf('Max-eig trial (exact)\n');
results = [];

t0 = tic;
[s0,A0] = max_eig_matrix(P.A,P.s,2*(p_order+1));
results.pt = s0;
R = chol(A0);
results.dt_precond = toc(t0);

pcon1=@(x) galerkin_preconditioner(R',x);
pcon2=@(x) galerkin_preconditioner(R,x);

S = instrumented_solver('minres','maxiter',maxiter,'tol',tol,...
    'precond',{pcon1,pcon2});

t0 = tic;
X = spectral_galerkin(P.A,P.b,P.s,p_order,'verbose',1,'lowmem',1,...
    'solver',S.solver,'parallel_rhs',1);
dt = toc(t0);    

[flag resvec] = S.stats();
results.X = X;
results.dt = dt;
results.flag = flag;
results.resvec = resvec;
results.hist = S.histguess();
results.matvectimes = S.matvectimes();

all_results.maxeig_exact = results;
save(sprintf('partial.mat'),'all_results');

%%
% 5. Max eigenvalue (fmincon)

fprintf('Max-eig trial (fmincon)\n');
results = [];

t0 = tic;
f = @(s) -max(abs(eigs(P.A(s))));
s0 = fmincon(f,rand(1,d),[],[],[],[],[P.s.l],[P.s.r],[],...
    optimset('Algorithm','active-set'));
results.pt = s0;
A0 = P.A(s0);
R = chol(A0);
results.dt_precond = toc(t0);

pcon1=@(x) galerkin_preconditioner(R',x);
pcon2=@(x) galerkin_preconditioner(R,x);

S = instrumented_solver('minres','maxiter',maxiter,'tol',tol,...
    'precond',{pcon1,pcon2});

t0 = tic;
X = spectral_galerkin(P.A,P.b,P.s,p_order,'verbose',1,'lowmem',1,...
    'solver',S.solver,'parallel_rhs',1);
dt = toc(t0);    

[flag resvec] = S.stats();
results.X = X;
results.dt = dt;
results.flag = flag;
results.resvec = resvec;
results.hist = S.histguess();
results.matvectimes = S.matvectimes();

all_results.maxeig_fmincon = results;
save(sprintf('partial.mat'),'all_results');

%%
% 6. Min eigenvalue (fmincon)

fprintf('Min-eig trial (fmincon)\n');
results = [];

t0 = tic;
f = @(s) min(abs(eigs(P.A(s),6,'SA')));
s0 = fmincon(f,rand(1,d),[],[],[],[],[P.s.l],[P.s.r],[],...
    optimset('Algorithm','active-set'));
results.pt = s0;
A0 = P.A(s0);
R = chol(A0);
results.dt_precond = toc(t0);

pcon1=@(x) galerkin_preconditioner(R',x);
pcon2=@(x) galerkin_preconditioner(R,x);

S = instrumented_solver('minres','maxiter',maxiter,'tol',tol,...
    'precond',{pcon1,pcon2});

t0 = tic;
X = spectral_galerkin(P.A,P.b,P.s,p_order,'verbose',1,'lowmem',1,...
    'solver',S.solver,'parallel_rhs',1);
dt = toc(t0);    

[flag resvec] = S.stats();
results.X = X;
results.dt = dt;
results.flag = flag;
results.resvec = resvec;
results.hist = S.histguess();
results.matvectimes = S.matvectimes();

all_results.mineig_fmincon = results;
save(sprintf('partial.mat'),'all_results');


%%
% Save results
save(sprintf('%s.mat',pname),'all_results');
