% Copyright (c) 2010 Kristopher L. Kuhlman (klkuhlm at sandia dot gov)
% 
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, including without limitation the rights
% to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
% copies of the Software, and to permit persons to whom the Software is
% furnished to do so, subject to the following conditions:
% 
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software.
% 
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
% AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
% OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
% THE SOFTWARE.
%
% this function is part of the set of MATLAB scripts 
% comprising the implementation of a universal kriging
% program that computes estimates of potential via
% cokriging, taking into account boundary conditions.
%
% Kuhlman, K.L., and E. Pardo-Igúzquiza, 2010. Universal
% cokriging of hydraulic heads accounting for boundary 
% conditions, Journal of Hydrology, 384(1–2), 14–25.
% http://dx.doi.org/10.1016/j.jhydrol.2010.01.002

clear;

%% Hard-wired parameters to set 
% <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

% small displacement to simulate deriv
m.rho = 10.0; 

% are boundary conditions being used? (1=true, 0=false)
bc_flag = 0;

% type of trend to incorporate 
% (0=constant (aka O.K.), 2=linear, 5=quadratic,)
L = 2;

m.type = 3;  % model 1=Gaussian, 2=generalized Cauchy, 3=Matern

% covariance model parameters
m.a     = 2048.215;   % range, 
m.sigsq = 641.1818;   % sill
if m.type == 2
    m.alpha = 1.5;
elseif m.type == 3
    m.nu = 2.2;
end

m.nugh = 0.7035;  % head direct-covariance nugget
m.nugc = 0.0;   % head gradient cross-covariance nugget (always zero?)
m.nugg = 0.0;  % gradient direct-covariance nugget (zero helps enforce BC)

% output mesh related things
% nx = 56;  ny = 77;
% xmin = 601750.0;  ymin = 3566550.0;
% dx = 400.0; dy = 400.0;
% nx = 40; ny = 40;
% xmin = 0.0; ymin = 0.01;
% dx = 0.2499; dy = 0.2499;

dx = 500.0; dy = 500.0;
xmin = 4000.0; ymin = 5000.0;
xmax = 31000.0; ymax = 26000.0;

%% uniform grid of calculation points
[X,Y] = meshgrid(xmin:dx:xmax,ymin:dy:ymax);             
D = numel(X);

nx = size(X,2);
ny = size(X,1);

% observation points in a long x,y vector
d.estlocs(1:D,1:2) = [reshape(X,D,1),reshape(Y,D,1)];

% the square-edged model domain
bound = load('eulogio_case/domain_boundary.dat');

% only krig at locations inside domain
IN = inpolygon(X,Y,bound(:,1),bound(:,2));

% >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
% ((((((((((((((((((((((((((((((((((((((((((((((((((
%% read / prepare head observations

z = load('eulogio_case/MAY1970.DAT');
N = size(z,1);
d.zlocs = z(1:N,1:2);
d.zdata = z(1:N,3);
% ))))))))))))))))))))))))))))))))))))))))))))))))))

% [[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[
%%  read / prepare BC "observations"

if bc_flag == 1
%     nfclocs = load('test01_nslocs.dat');   % no-flow boundary data
%     chclocs = load('test02_ewlocs.dat'); % constant head boundary data
%     chclocs = [];
%     wclocs = [nfclocs; chclocs];
    wclocs = load('eulogio_case/locs.dat');
    M = size(wclocs,1);

    % unit vector pointing normal to no-flow boundary
    % or parallel to specified head boundary
%     nfdata = load('test01_nfvecs.dat');   % no-flow vectors
%     chdata = load('test02_chvecs.dat');  % constant head vectors
%     %chdata = [];
%     d.wdata = [nfdata;chdata];
    d.wdata = load('eulogio_case/vecs.dat');

    % compute pairs of points straddling boundary,
    % separated by 2*rho, with orientation wdata
    d.wlocs(1:M,1:2,1:2) = cat(3, wclocs(1:M,1:2) + d.wdata(1:M,1:2)*m.rho, ...
                                  wclocs(1:M,1:2) - d.wdata(1:M,1:2)*m.rho);
elseif bc_flag == 0
    M=0;
else
    error('invalid choice for bc_flag; no bc=0, with bc=1')
end
% ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]

% test -- include both data at one point
% zlocs(21,1:2) = wclocs(20,1:2);
% % zlocs(22,1:2) = wclocs(11,1:2);
% N=21;
% % zdata(21:22) = sum(zlocs(21:22,1:2),2);
% zdata(21) = sum(zlocs(21,1:2),2);

[z0,MSE] = solve_fd_gradient_system(N,M,L,D,d,m,bc_flag);

% make output easier to plot
z0 = reshape(z0,ny,nx);
MSE = reshape(MSE,ny,nx);

% to blank out regions outside boundary
z0(~IN) = NaN;
MSE(~IN) = NaN;

%% plot
figure;
% subplot(1,2,1)
[C,h]=contour(X,Y,z0,[500:5:600,610:10:700]);
colormap([0,0,0])
clabel(C,h)
% axis([-0.1,10.1,-0.1,10.1])
axis([4000,31000,6000,25500])
set(gca,'YTick',6000:6000:24000);
set(gca,'YTickLabel',{'6000','12000','18000','24000'});
set(gca,'XTick',4000:6000:28000);
set(gca,'XTickLabel',{'4000','10000','16000','22000','28000'});
% cmap = colormap;
% title('estimated head');
hold on;
plot(d.zlocs(:,1),d.zlocs(:,2),'r*');
plot(bound(:,1),bound(:,2),'k-')
if bc_flag == 1
    quiver(d.wlocs(:,1),d.wlocs(:,2),d.wdata(:,1),d.wdata(:,2),0.1,'LineWidth',4);
end
% colorbar;
daspect([1 1 1]);
xlabel('Easting (m)','fontsize',12)
ylabel('Northing (m)','fontsize',12)
hold off;

figure;
% subplot(1,2,2)
[C,h]=contour(X,Y,MSE);
clabel(C,h)
colormap([0,0,0])
% axis([-0.1,10.1,-0.1,10.1])
axis([4000,31000,5000,25500])
title('MSE');
hold on;
plot(d.zlocs(:,1),d.zlocs(:,2),'r*');
if bc_flag == 1
    quiver(d.wlocs(:,1),d.wlocs(:,2),d.wdata(:,1),d.wdata(:,2),0.1,'LineWidth',4);
end
% colorbar;
daspect([1 1 1]);
xlabel('Easting (m)')
ylabel('Northing (m)')
hold off;

