clear;

% default values
n = 500;                      % number of interior nodes
E = 1;                        % Number of experiments
k = 2;                        % Number of measurements
c = 3;                        % Initial guess for s(1)
method = 2;                   % 1: gradient, 2: GN

%%%%%%%%% asking for input or not %%%%%%%%%%%%%%%%%%%
%     E k c n method
IN = [0 1 1 0 1];
if IN(1); n = input('No. of interior nodes = '); end
if IN(2); E = input('No. of experiments  = '); end
if IN(3); k = input('No. of measurements = '); end
if IN(4); c = input('initial s(1) = '); end
if IN(5); method = input('method = '); end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

h = 1/(n + 1);                % step size
bvsE = [1:E; E:-1:1];         % experiment boundary values
B = round((1:k)/(k+1)*n);     % uniformly distributed

% generate data
ss = [1/10; 10];              % true solution
[As,bE] = zgetAb(ss,bvsE,n);  % A(ss) and b(bvs)
uE = As \ bE; cE = uE(B,:);   % measured data
s0 = [c; 1];                  % initial guess

tol = 1.e-6; maxit = 400;
if method == 1; maxit = 800; end

t0 = cputime;     
[sz,iterz,nfz,fhz,ghz] = ...
    zsolver(@zinvprJ,s0,tol,maxit,method,n,bvsE,cE,B);
tz = cputime - t0;
rel_err = norm(sz-ss)/norm(ss);
fprintf('zsolver: [n E k c method] = [%i %i %i %i %i], s* = [%g %g]\n',...
    n,E,k,c,method,ss); 
fprintf('sz: [%8.6f %8.6f], rel_err: %6.2e, cput: %gs\n\n',sz,rel_err,tz);

% call your solver 
if exist('mysolver.m','file');
t0 = cputime;
[sm,iterm,nfm,fhm,ghm] = ...
    mysolver(@myinvprJ,s0,tol,maxit,method,n,bvsE,cE,B);
tm = cputime - t0;
rel_err = norm(sm-ss)/norm(ss);
fprintf('mysolver: [n E k c method] = [%i %i %i %i %i]\n',n,E,k,c,method); 
fprintf('sm: [%8.6f %8.6f], rel_err: %6.2e, cput: %gs\n\n',sm,rel_err,tm);
end