% Script: TestInvSteady
% Programmer: T. Shores
% Date: 7/02/02
% Description:  Routine for solving the inverse problem of
% determining the dispersion coefficient D(x) below, given
% a set of (noisy) datapoints of the soluton c(x). Following
% BVP is solved by a pseudo Sinc-Galerkin method.
% -(D(x)*c')' + (v(x)*c)' + lam(x)*c = f(x), 0 < x < infty.
%  A*c(0) + B*c'(0) = G
%  u(infty) = 0
% A variational approach is based on this form of the problem:
% For all test functions u(x)
% int(D*c'*u'-v*c*u'-lam*c*u,x=0..infty) = int(f*u,x=0..infty)
%
% You must run the file sincfcn or have all these fcns available.
% You must define a test problem by COMPILING TestInvProbx.m

% GLOBAL PARAMETERS

global tndx = 0;
global hfactor = 1.0; % should be between 1 and 2, 1 is default
% sinc setup parameters for forward problem
global alpha;     % left exponential rate factor
global bta;       % right exponential rate factor
global d;         % width of strip
global mm;        % starting left index for nodes
global nn;        % ending right index for nodes
global hx;        % sinc stepsize
global nodes;     % sinc nodes
global I1;        % sinc matrix I^{1}
global vnodes;    % = arrayf('v',nodes);
global phinodes   % = arrayf('phi',nodes);
global phipnodes  % = arrayf('phip',nodes);
global gamnodes   % = arrayf('gam',nodes);
global gampnodes  % = arrayf('gamp',nodes);
global wnodes     % = arrayf('w',nodes);
global wpnodes    % = arrayf('wp',nodes);
global lamnodes   % = arrayf('lam',nodes);
global pnodes     % = arrayf('p',nodes);
global ppnodes    % = arrayf('pp',nodes);
global fnodes     % = arrayf('f',nodes);
global q0nodes    % = arrayf('q0',nodes);
global q0pnodes   % = arrayf('q0p',nodes);
global q1nodes    % = arrayf('q1',nodes);
global q1pnodes   % = arrayf('q1p',nodes);
global q2nodes    % = arrayf('q2',nodes);
global q2pnodes   % = arrayf('q2p',nodes);

% regularization setup parameters for inverse problem
global alp;       % Tikhonov regularization parameter
global dtnodes; % data nodes
global cdtnodes;% data vector (values at dtnodes)
global N          % length of data vector
global H          % data interpolation matrix

% DEFINE THE VARIABLE PARAMETER(S)

function retval = D(x)
% usage: D(x)
% description: this computes the fcn D(x)

retval =  1 - exp(-2*x)/2;

endfunction

% DEFINE THE WEIGHT FUNCTIONS

function retval = gam(x)
% usage: gam(x)
% description: this computes the weight fcn gam(x)
% for deriv approx's

%mu = 0.0;
%retval = x^(1+mu)/(10+x)^2;
%retval = x*sech(x);
retval = x/(5+x)^2;

endfunction

%

function retval = gamp(x)
% usage: gamp(x)
% description: this computes the fcn gam'(x)

%mu = 0.0;
%retval = x^mu*((10+x)*(1+mu)-2*x)/(10+x)^3;
%retval = (1-x*tanh(x))*sech(x);
retval = (5-x)/(5+x)^3;

endfunction

%

function retval = w(x)
% usage: w(x)
% description:  this is the integration weight
% Warning: If you change this, you must also change wp

retval = gam(x);

endfunction

%

function retval = wp(x)
% usage: wp(x)
% description: this is the derivative of the integration weight

retval = gamp(x);

endfunction

% DEFINE THE TRANSFORMATION FUNCTIONS

function retval = q0(x)
% usage: q0(x)
% description: this computes the fcn q0(x)

retval = p(x)*x/(x+1);

endfunction

%

function retval = q0p(x)
% usage: q0p(x)
% description: this computes the fcn q0'(x)

retval = pp(x)*x/(x+1) + p(x)/(x+1)^2;

endfunction

%

function retval = q1(x)
% usage: q1(x)
% description: this computes the fcn q1(x)

retval = p(x)*(-(pp(0)-1)*x+1)/(x+1);

endfunction

%

function retval = q1p(x)
% usage: q1p(x)
% description: this computes the fcn q1'(x)

retval = pp(x)*(-(pp(0)-1)*x+1)/(x+1) - p(x)*pp(0)/(x+1)^2;

endfunction


%

function retval = q2(x)
% usage: q2(x)
% description: this computes the fcn q2(x)

retval = 1 - 10/(10+x*x*exp(x));

endfunction

%

function retval = q2p(x)
% usage: q2p(x)
% description: this computes the fcn q2'(x)

retval = 10*(2*x+x*x)*exp(x)/(10+x*x*exp(x))^2;

endfunction

% DEFINE THE INTEGRANDS FOR BILINEAR FORM T

function retval = T0(x)
% usage: T0(x)
% description: this computes the integrand fcn of T(q0,w*Skh(phi(x)))
global hx;
global tndx;

retval = (D(x)*q0p(x) - v(x)*q0(x))* ...
(wp(x)*Skh(tndx,hx,phi(x)) + w(x)*phip(x)*Skhp(tndx,hx,phi(x))) + ...
lam(x)*q0(x)*Skh(tndx,hx,phi(x));

endfunction

%

function retval = T1(x)
% usage: T1(x)
% description: this computes the integrand fcn of T(q1,Skh)
global hx;
global tndx;

retval = (D(x)*q1p(x) - v(x)*q1(x))* ... 
(wp(x)*Skh(tndx,hx,phi(x)) + w(x)*phip(x)*Skhp(tndx,hx,phi(x))) + ...
lam(x)*q1(x)*Skh(tndx,hx,phi(x));

endfunction


% DEFINE THE SHIFTED UNKNOWN FUNCTION

function retval = csup0(x)
% usage: csup0(x)
% description: this computes the fcn csup0(x)

retval = csoln(x) - q0(x) - csoln(0)*q1(x);

endfunction

% DEFINE THE CONFORMAL MAPPINGS

function retval = xphi(x)
% usage: phi(x)
% description: this computes the fcn phi(x)

retval = log(sinh(x));

endfunction

%

function retval = xphip(x)
% usage: phip(x)
% description: this computes the fcn phi'(x)

retval = 1/tanh(x);

endfunction

%

function retval = xpsi(x)
% usage: psi(x)
% description: this computes the fcn psi(x)

retval = log(exp(x)+sqrt(exp(2*x)+1));

endfunction


function retval = phi(x)
% usage: phi(x)
% description: this computes the fcn phi(x)

retval = log(x);

endfunction

%

function retval = phip(x)
% usage: phip(x)
% description: this computes the fcn phi'(x)

retval = 1/x;

endfunction

%

function retval = psi(x)
% usage: psi(x)
% description: this computes the fcn psi(x)

retval = exp(x);

endfunction

% BUILD AND SOLVE THE SYSTEM

function c = bvp(pDnodes)
% usage: c = bvp(pDnodes)
% description: this routine solves the system in question,
% using the functions defined in this file and return the
% the column vector [cpatzero;catzero;c0;0] whose first entry
% is the flux of solution at zero, second is value of
% solution at zero, last is value at infinity and remaining
% entries are the reduced soluton at sinc nodes. This enables
% the construction of exact solution c(x) in the form
%
%     c(x) = c0(x) + cpatzero*q0(x) + catzero*q1(x)
%
% where c0(x) is expressed as a weighted sinc sum.
% which is essential for interpolation purposes.
% It is assumed that the *input* pDnodes is also in this
% format.

% local variables
% nn: finishing index for nodes
% coef: coefficient matrix
% v0,v1: right hand side vectors
% w0,w1: right hand side parts of soln
% c0: proposed solution node vector for ctilde(x)
% errvec:  error vector of c0-ctilde at nodes

% global variables:
global A;
global B;
global G;
global edecay;
global hfactor;
global alpha;
global bta;
global d;
global mm;
global nn;
global nodes;
global hx;
global I1;
global vnodes;    % = arrayf('v',nodes);
global phinodes   % = arrayf('phi',nodes);
global phipnodes  % = arrayf('phip',nodes);
global gamnodes   % = arrayf('gam',nodes);
global gampnodes  % = arrayf('gamp',nodes);
global wnodes     % = arrayf('w',nodes);
global wpnodes    % = arrayf('wp',nodes);
global lamnodes   % = arrayf('lam',nodes);
global pnodes     % = arrayf('p',nodes);
global ppnodes    % = arrayf('pp',nodes);
global fnodes     % = arrayf('f',nodes);
global q0nodes    % = arrayf('q0',nodes);
global q0pnodes   % = arrayf('q0p',nodes);
global q1nodes    % = arrayf('q1',nodes);
global q1pnodes   % = arrayf('q1p',nodes);
global q2nodes    % = arrayf('q2',nodes);


% first, upack the input
[xDnodes,Dp0,Dinf] = unpackc(pDnodes);
Dnodes = xDnodes(2:length(xDnodes));
% set up coefficient matrix
% This is the discretization of T(c^0,S(i,hx)(phi(x))) coef matrix
%coef1 = hx*diag(Dnodes.*gampnodes.*wpnodes./gamnodes./phipnodes) + ...
%I1*diag(Dnodes.*gampnodes.*wnodes./gamnodes) - ...
%diag(Dnodes.*gamnodes.*wpnodes)*I1*diag(1./gamnodes) - ...
%I1*diag(Dnodes.*gamnodes.*phipnodes.*wnodes)*1/hx*I1*diag(1./gamnodes);
%coef2 = hx*diag(vnodes.*wpnodes ./phipnodes) + I1*diag(vnodes.*wnodes);
%coef3 = hx*diag(lamnodes .*wnodes./phipnodes);
%coef = coef1 - coef2 + coef3;
% an alternate grouping...made no significant difference, but faster
coef = hx*diag(((Dnodes.*gampnodes./gamnodes-vnodes).*wpnodes + ...
lamnodes.*wnodes)./phipnodes) - ...
(diag(Dnodes.*gamnodes.*wpnodes)+I1*diag( ...
Dnodes.*gamnodes.*phipnodes.*wnodes)/hx)*I1*diag(1./gamnodes) + ...
I1*diag(Dnodes.*gampnodes.*wnodes./gamnodes-vnodes.*wnodes);

% This is the discretization of T(qi,w(x)S(i,hx)(phi(x)))
v0 = Dnodes.*q0pnodes - vnodes.*q0nodes;
Tq0 = hx*(v0.*wpnodes + lamnodes.*wnodes.*q0nodes)./phipnodes + ...
I1*(v0.*wnodes);
v0 = Dnodes.*q1pnodes - vnodes.*q1nodes;
Tq1 = hx*(v0.*wpnodes + lamnodes.*wnodes.*q1nodes)./phipnodes + ...
I1*(v0.*wnodes);

% Set up the rhs vector
v0 = hx*wnodes.*fnodes./phipnodes;
% solve system and assign values for the function to return
ws = coef\[v0,Tq0,Tq1]; % this is w_s, w_0, w_1
w0 = ws(:,2);
w1 = ws(:,3);
ws = ws(:,1);
% compute c(0) and c'(0)
A1 = hx*sum(lamnodes.*(q1nodes-w1)./phipnodes)-v(0);
B1 = hx*sum(lamnodes.*(q0nodes-w0)./phipnodes)+xDnodes(1);
G1 = hx*sum((fnodes-lamnodes.*ws)./phipnodes);
catzero = [A,B;A1,B1]\[G;G1];
cpatzero = catzero(2);
catzero = catzero(1);
% finish the calculation for c0 (hence c indirectly)
c0 = ws - cpatzero*w0 - catzero*w1;
c = [cpatzero;catzero;c0;0];

endfunction

% Construct the Tikhonov functional

function retval = TikR(pDnodes)
% usage: R = TikR(pDnodes)
% description: Input is packed value of D at nodes, output is
% vector R such that R^T*R is Tikhonov functional at pDnodes.

% local variables

% global variables
global alpha;
global bta;
global d;
global mm;
global nn;
global hx;
global nodes;
global I1;
global alp;
global dtnodes;
global cdtnodes;
global N;
global H          % data interpolation matrix
global vnodes;    % = arrayf('v',nodes);
global phinodes   % = arrayf('phi',nodes);
global phipnodes  % = arrayf('phip',nodes);
global gamnodes   % = arrayf('gam',nodes);
global gampnodes  % = arrayf('gamp',nodes);
global wnodes     % = arrayf('w',nodes);
global wpnodes    % = arrayf('wp',nodes);
global lamnodes   % = arrayf('lam',nodes);
global pnodes     % = arrayf('p',nodes);
global ppnodes    % = arrayf('pp',nodes);
global fnodes     % = arrayf('f',nodes);
global q0nodes    % = arrayf('q0',nodes);
global q0pnodes   % = arrayf('q0p',nodes);
global q1nodes    % = arrayf('q1',nodes);
global q1pnodes   % = arrayf('q1p',nodes);
global q2nodes    % = arrayf('q2',nodes);
global q2pnodes   % = arrayf('q2p',nodes);
%retval = [H*bvp(pDnodes)-cdtnodes;...
%          sqrt(alp)*(Dnodes-1)];
% the idea is that
% ||D||^2 = S_1*D'(0)^2 + S_2*D(0)^2 + S_3*D(inf)^2 +
%             S_4*||Dred||^2 + S_5*||Dred'||^2
% with coefficients positive
Sc = [0.001;0.001;0.001;1;1];
D0pbias = 0;
D0bias = 0;
Dinfbias = 0;
mD = length(pDnodes);
vtmp = sqrt(hx./phipnodes);
Dtmp = pDnodes(3:mD-1);
retval = [H*bvp(pDnodes)-cdtnodes;...
          sqrt(alp)*[Sc(1)*(pDnodes(1)-D0pbias);...
	  Sc(2)*(pDnodes(2)-D0bias);...
	  Sc(3)*(pDnodes(mD)-Dinfbias);...
	  Sc(4)*vtmp.*Dtmp;...
	  Sc(5)*vtmp.*(Dtmp.*gampnodes./gamnodes - (1/hx)*...
	  gamnodes.*phipnodes.*(I1*(Dtmp./gamnodes)))]];

endfunction

% Some utilities

function v = packc(c,cp0,cinf)
% usage: v = packc(c,cp0,cinf)
% description: Input expanded node vector c, c'(0), c(inf)
% and output the packed form of fcn c.

global mm;
global nn;
global nodes;
global hx;
global I1;
global q0nodes    % = arrayf('q0',nodes);
global q0pnodes   % = arrayf('q0p',nodes);
global q1nodes    % = arrayf('q1',nodes);
global q1pnodes   % = arrayf('q1p',nodes);
global q2nodes    % = arrayf('q2',nodes);
global q2pnodes   % = arrayf('q2p',nodes);
global q2pnodes   % = arrayf('q2p',nodes);

v = [cp0;c(1);c(2:mm+nn+2) - cp0*q0nodes - ...
           c(1)*q1nodes - cinf*q2nodes;cinf];

endfunction

%

function [c,cp0,cinf] = unpackc(v)
% usage: [c,cp0,cinf] = unpackc(v)
% description: Input packed form of vector v
% that is, v = [vp0;v0;vred;vinf], and output
% the expanded node vector c, c'(0), c(inf).

global q0nodes    % = arrayf('q0',nodes);
global q1nodes    % = arrayf('q1',nodes);
global q2nodes    % = arrayf('q2',nodes);

mv = length(v); % really mm+nn+4
cp0 = v(1);
cinf = v(mv);
c = v(1)*q0nodes + v(2)*q1nodes + ...
                v(3:mv-1) + cinf*q2nodes;
c = [v(2);c];

endfunction

%

% An experiment (run a testproblem file first):
% sinc method setup parameters
alpha = 1;
bta = 1;
d = pi/3;
mm = 40;
% regularization setup parameters
alp = .001;
nonlinearLS;
% compute remaining sinc setup globals
sincfcn;
hx = sqrt(hfactor*pi*d/(alpha*mm));
% if we have exponential decay, make use of the economization
if (edecay)
  nn = round(log(alpha*mm*hx/bta)/hx);
else
  nn = round(alpha/bta*mm);
end
% make sure even I^-1 used
if (rem(mm+nn+1,2)>0)
  nn=nn+1;
end;
% the essential nodes
nodes = arrayf('psi',[-mm:nn]*hx)';
xnodes = [0;nodes]; % the 'x' is for 'extended'
% construct I^{1} just once
m = mm + nn + 1;
I1 = Idn(1,m);
% in the interest of speed, these (bad) globals
% set up coefficient matrix
vnodes = arrayf('v',nodes);
phinodes = arrayf('phi',nodes);
phipnodes = arrayf('phip',nodes);
gamnodes = arrayf('gam',nodes);
gampnodes = arrayf('gamp',nodes);
wnodes = arrayf('w',nodes);
wpnodes = arrayf('wp',nodes);
lamnodes = arrayf('lam',nodes);
pnodes = arrayf('p',nodes);
ppnodes = arrayf('pp',nodes);
fnodes = arrayf('f',nodes);
q0nodes = arrayf('q0',nodes);
q0pnodes = arrayf('q0p',nodes);
q1nodes = arrayf('q1',nodes);
q1pnodes = arrayf('q1p',nodes);
q2nodes = arrayf('q2',nodes);
q2pnodes = arrayf('q2p',nodes);
%
% compute remaining inverse problem globals
dtnodes = linspace(0,12,m)';
%dtnodes = xnodes;
% calculate the point evaluation matrix H
% Of course, we assume data node points are
% nonnegative, distinct and in increasing order.
N = length(dtnodes);
if (dtnodes(1)==0) % avoid certain fcn evalns at 0
  startrow = 2;
  H = zeros(N-1,m);
else
  startrow = 1;
  H = zeros(N,m);
end
gamdtnodes = arrayf('gam',dtnodes(startrow:N));
phidtnodes = arrayf('phi',dtnodes(startrow:N));
for jj = 1:m
  H(:,jj) = gamdtnodes.*Skh(jj-1-mm,hx,phidtnodes)/...
     gam(nodes(jj));
  end
% now prepend 1st row if needed
if (startrow == 2)
  H = [zeros(1,m);H];
end
H = [arrayf('q0',dtnodes),arrayf('q1',dtnodes),H,...
     arrayf('q2',dtnodes)];
%
% "Measurements": exact solution at data nodes plus noise
cdtnodes = arrayf('csoln',dtnodes).*(1+0.02*rand(size(dtnodes)));
% OK, here goes the solution, which we pack
Dnodestrue = packc(arrayf('D',xnodes),D(0),Dp(0));
x0 = [2.5;2;zeros(m,1);2];
nu0 = 0.1;
iterlim = 70;
h = 0.0001;
global storage;
storage = zeros(length(x0),iterlim+1);
%[Dnodesx,resid,k] = levmart('TikR',x0,nu0,iterlim,h);
%[Dnodesx,resid,k] = levmar('TikR',x0,nu0,iterlim,h);
[Dnodesx,resid,kk,nr] = bfgs('TikR',x0,nu0,iterlim,h);
%[Dnodesx,resid,kk,nr] = quasinewton('TikR',x0,nu0,iterlim,h);
% storage contains packed approximants to D and true solution
storage(:,iterlim+1) = Dnodestrue;
% compute the corresponding solutions
cstorage = zeros(m+1,kk+1);
dstorage = zeros(m+1,kk+1);
for jj=1:kk
  cstorage(:,jj) = unpackc(bvp(storage(:,jj)));
  dstorage(:,jj) = unpackc(storage(:,jj));
end
cstorage(:,kk+1) = unpackc(bvp(Dnodestrue));
dstorage(:,kk+1) = unpackc(Dnodestrue);
%clearplot;grid;hold on
%plot(xnodes,dstorage(:,kk-1:kk))
%clearplot;grid;
%plot(xnodes,cstorage(:,kk-1:kk))
%plot(xnodes,cstorage(:,kk-1:kk))
ctrue = unpackc(bvp(Dnodestrue));
ctrue = ctrue(2:length(ctrue));
cprox = unpackc(bvp(Dnodesx));
cprox = cprox(2:length(cprox));
cerror = sqrt(sum(hx./phipnodes.*(ctrue-cprox).^2))