function [xx, yy, m, dp, actrms] = OccamMap(targrms,x,y,d,s, x0, x1, nx, y0, y1, ny, masksize);
%
% [m, actrms] = OccamMap(targrms,x,y,d,s, x0, x1, nx, y0, y1, ny)
% Occam's inversion of for a smooth 2D map
%
% targrms = target misfit, normalized root mean square  
%                  (if empty, OccamMap will try to choose something sensible)
% x = vector of data x coords
% y = vector of data y coords
% d = vector of data values
% s = vector of errors in d
%         s can be empty, in which case 1% erros will be assigned
% x0 = low x value of model (left edge)
% x1 = high x value of model (right edge)
% nx = number of x values in model 
% y0 = low y value of model (bottom edge)
% y1 = high y value of model (top edge)
% ny = number of y values in model
% masksize = radius for filling plot with NaNs to mask off areas with no data.
%
% on output
% xx = x axis values
% yy = y axis values
% m = model 
% dp = predicted data
% actrms = actual RMS of model response

% Steven Constable, Scripps Institution of Oceanography
% version 2.0, July 2018

t0 = clock;

% test to see if data are all on map. (Code crashes if data are off the map.)
if min(x) < x0
  disp('Sorry, x0 not low enough for data set.  Please try again.') 
  return
end
if max(x) > x1
  disp('Sorry, x1 not large enough for data set.  Please try again.') 
  return
end
if min(y) < y0
  disp('Sorry, y0 not low enough for data set.  Please try again.') 
  return
end
if max(y) > y1
  disp('Sorry, y1 not large enough for data set.  Please try again.') 
  return
end

% make the x and y coords of the map using the values passed 
xx = linspace(x0,x1,nx);
yy = linspace(y0,y1,ny);  % note that y goes UP in space, DOWN in the vector

% assign 1 percent errors if errors are not supplied
if (isempty(s)) 
  s = 0.01*d;  
end

M = length(d);      % the length of the data
N = ny*nx;          % the length of the model

% make the roughening matrix, R (but here called S for sparse)
k = 0;          % this is an entry counter for the S matrix
l = 0;          % this is a row-counter for the S matrix
ssize = 2*(nx*(ny-1) + (nx-1)*ny);
Si = zeros(ssize,1);  % row index of non-zero elements of R (roughening matrix) 
Sj = zeros(ssize,1);  % column index
Ss = zeros(ssize,1);  % value of R, usually -1 or 1, but could be scaled by dx and dy

for i = 0:nx-1  % vertical roughness measure
    for j = 1:ny-1
        k = k+1;
        l = l+1;
        Si(k) = l;
        Sj(k) = i*ny + j;
        Ss(k) = -1;
        k = k+1;
        Si(k) = l;
        Sj(k) = i*ny + j + 1;
        Ss(k) = 1;
    end
end

for i = 0:nx-2  % horizontal roughness measure
     for j = 0:ny-1
         k = k+1;
         l = l+1;
         Si(k) = l;
         Sj(k) = i*ny+1 + j;
         Ss(k) = -1;
         k = k+1;
         Si(k) = l;
         Sj(k) = (i+1)*ny+1 + j;
         Ss(k) = 1;
     end
 end

 del = sparse(Si, Sj, Ss); 
 delsq = del'*del;
 clear del

% make G, the forward funtion matrix

sqrt2 = sqrt(2);

% start by finding the data positions in map coordinates
idyf = interp1(yy, 1:ny, y);  % use linear interpolation, the default
idxf = interp1(xx, 1:nx, x);  % note that the units are not meters, but indices 

% Find the lower left corner of the box containing the data 
idy = floor(idyf);
idx = floor(idxf);

% and the distance into the box in the two directions
xfrac = idxf - idx;
yfrac = idyf - idy;

G = zeros(M, N);

for i = 1:M
% the indices of the corners of the box containing the datum
    j1 = ny*(idx(i)) + idy(i); % bottom left  (space coords)
    j2 = ny*(idx(i)+1) + idy(i);   % bottom right 
    j3 = ny*(idx(i)) + idy(i)+1;  % top  left
    j4 = ny*(idx(i)+1) + idy(i)+1;   %  top right 

% distances from the corners    
    wt1 = sqrt2-sqrt(xfrac(i)^2 + yfrac(i)^2); 
    wt2 = sqrt2-sqrt((1-xfrac(i))^2 + yfrac(i)^2);
    wt3 = sqrt2-sqrt(xfrac(i)^2 + (1-yfrac(i))^2);
    wt4 = sqrt2-sqrt((1-xfrac(i))^2 + (1-yfrac(i))^2);
    
    summ = wt1 + wt2+ wt3+ wt4; % because we're going to normalized the weights
    G(i, j1) = wt1/summ;
    G(i, j2) = wt2/summ;
    G(i, j3) = wt3/summ;
    G(i, j4) = wt4/summ;
end

G = sparse(G);

W = diag(1./s);
W = sparse(W);
G = W*G;
ds = sparse(d);
ds = W*ds;

logmu = -5:1:7;
mu = 10.^logmu;
k = length(mu);

m = inv(mu(1)*delsq + G'*G)*(G'*ds);
r(1) = rms(ds - G*m);
m = reshape (m, ny, nx);

for i = 1:k
       m = inv(mu(i)*delsq + G'*G)*(G'*ds);
       r(i) = rms(ds - G*m);
       m = reshape (m, ny, nx);
end

%figure; loglog(mu, r , '*' );hold on
disp( ['Limits of misfit: ' num2str(r(1)) '  ' num2str(r(end))]) 
 
if isempty(targrms)
    %targrms = full((log10(max(r)) + log10(min(r)))/2)    
    targrms = full(log10(max(r)));
    targrms = (10^targrms)/2;
end
if (targrms < 0)
    fact = abs(targrms);
    %targrms = full((log10(max(r)) + log10(min(r)))/2)    
    maxrms = full(log10(max(r)));
    targrms = fact*(10^maxrms);
end

if (r(1) > targrms) 
    disp(' Appears that requested RMS is too small')
    targrms = full(log10(max(r))/2);
    targrms = 10^targrms;
end

if (r(end) < targrms) 
    disp(' Appears that requested RMS is too large')
    targrms = full(log10(max(r))/2);
    targrms = 10^targrms;
end


% interpolate to find the mu that will produce the required RMS misfit
  mmu = interp1( full(r), logmu, targrms,'pchip');
  mmu = 10.^mmu;
  m = inv(mmu*delsq + G'*G)*(G'*ds);
  actrms = full(rms(ds - G*m));
  Ws = diag(s);
  Ws = sparse(Ws);   
  G = Ws*G;   % undo the weighting to generate a forward model
  dp = G*m;
  dp = full(dp);
  m = reshape (m, ny, nx);

 % now NaN out parts of the map too far from actual data

for i = 1:nx
    for j = 1:ny
        distance = sqrt((x - xx(i)).^2 + (y - yy(j)).^2);
        if (min(distance) > masksize)
              m(j,i) = NaN;
        end
     end
end
  
  
%loglog(mmu, actrms , 'r*' )

% plot a figure of the model

figure
pcolor(xx,yy,m);
shading flat;
axis equal
axis tight

hold on

colorbar

% add the data locations

plot (x, y, 'w.')

disp(['time taken = ' num2str(etime(clock,t0))]);

end



function r = rms(v)
n = length(v);
chi = sum(v.^2);
r = sqrt(chi/n);
return
end 
