function [xx,yy,mo, dp, actrms] = OccamMapSP(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 of self potential
% targrms = target misfit, root mean square weighted 
%                  (if empty, OccamMap will choose something sensible)
% x = vector of data x coords
% y = vector of data y coords
% d = vector of data values - for SP east and north gradients
% s = vector of errors in d
%         s can be empty, in which case 1% erros will be assigned
% x0 = low x value of model
% x1 = high x value of model
% nx = number of x values in model
% y0 = low y value of model
% y1 = high y value of model
% 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, November 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

% minimize the value of the model in the lower left corner (i.e. make the SP zero here)
k = k+1;
Si(k) = 2;
Sj(k) = 2;
Ss(k) = 1;

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

% find G

% start with 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);

% get the size of the cells
dx = xx(2) - xx(1);
dy = yy(2) - yy(1);
    wtx = 1/(2*dx);
    wty = 1/(2*dy);

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

G = zeros(M, N);
for i = 1:M/2
     % 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 
    
 % this is the x part of the data vector
    G(i, j1) =  -(1-yfrac(i))/dx;
    G(i, j2) = (1-yfrac(i))/dx;
    G(i, j3) =  -yfrac(i)/dx;
    G(i, j4) = yfrac(i)/dx;
% this is the y part of the data vector
    G(i+M/2, j1) =  -(1-xfrac(i))/dy;
    G(i+M/2, j2) =  -xfrac(i)/dy;
    G(i+M/2, j3) = (1-xfrac(i))/dy;
    G(i+M/2, j4) = xfrac(i)/dy;
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 the target rms is empty, use half the maximum misfit.  Largely untested.
if isempty(targrms)
    %targrms = full((log10(max(r)) + log10(min(r)))/2)    
    targrms = full(log10(max(r)));
    targrms = (10^targrms)/2;
end

% if the target rms is negative, use it's absolute value as a fraction of 
% maximum misfit.  I think this is totally untested.
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



  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;
  dp = G*m;
  dp = full(dp);
  dp = dp.*s;
  mo = full(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)
              mo(j,i) = NaN;
          end
      end
  end
  
  
%loglog(mmu, actrms , 'r*' )

figure
pcolor(xx,yy,mo);
shading flat;
axis equal
%axis tight
%caxis([8 20]);

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 
