/*
* C Source code
* This program performs ordinary block kriging for a 20 X 20 X 20
* three-dimensional grid around a set of 3-d data.
* Output is coordinates/estimates.
* Tail of output file contains estimation variance values computed
* using the methods described by Journel and Huijbregts (1978) and
* Isaaks and Srivastava (1989).
*
* References:
*
* Carr, J.R., and Palmer, J.A. 1993. Revisiting the accurate
* calculation of block-sample covariances using Gauss quadrature.
* Mathematical Geology 25, no. 5: 507-524.
*
* Deutsch, C.V. 1996. Correcting for negative weights in ordinary kriging.
* Computers & Geosciences 22, no. 7: 765-773.
*
* Isaaks, E.H., and Srivastava, R.M. 1989. Applied Geostatistics.
* Oxford University Press Inc., New York, NY.
*
* Journel, H.G., and Huijbregts, C.J. 1978. Mining Geostatistics.
* Academic Press, London. 500 p.
*
* Press, W.H., Teukolsky, S.A., Vetterling, W.T., and Flannery, B.P.
* 1996. Numerical Recipes in C, Second Edition. Cambridge University
* Press. 994 p.
*
* Written by Bart Faulkner, 1998
* RS Kerr Environmental Research Center
* PO Box 1198
* Ada, OK 74820
*
*/
#include "krig.h"
int ndata, numnodes_x, numnodes_y, numnodes_z;
int *indx;
double *east, *north, *elev, *bchem, *estdist, *estcovar, *vv, *b, var;
double **smpdist, **smpcovar, **a;
double ***massest;
/* Determine the size of the file and allocate external arrays */
int main(void) {
int i, j, entno;
double ea, no, el, bc;
FILE *ifp;
numnodes_x=10;
numnodes_y=10;
numnodes_z=10;
ifp=fopen("bchemfile.u","r");
for (i=0; fscanf(ifp, "%d",&entno) != EOF; ++i) {
fscanf(ifp,"%f %f %f %f\n",&ea,&no,&el,&bc);
ndata=i;
}
ndata += 1;
fclose(ifp);
indx = calloc(ndata+1, sizeof(int));
east = calloc(ndata, sizeof(double));
north = calloc(ndata, sizeof(double));
elev = calloc(ndata, sizeof(double));
bchem = calloc(ndata, sizeof(double));
estdist = calloc(ndata+1, sizeof(double));
estcovar = calloc(ndata+1, sizeof(double));
vv=calloc(ndata+1, sizeof(double));
b=calloc(ndata+1, sizeof(double));
smpdist = calloc(ndata+1, sizeof(double *));
smpcovar = calloc(ndata+1, sizeof(double *));
a = calloc(ndata+1, sizeof(double *));
for (i=0; i<ndata+1; ++i) {
smpdist[i] = calloc(ndata+1, sizeof(double));
smpcovar[i] = calloc(ndata+1, sizeof(double));
a[i] = calloc(ndata+1, sizeof(double));
}
massest = calloc(numnodes_x, sizeof(double **));
for (i=0; i<numnodes_x; ++i) {
massest[i] = calloc(numnodes_y, sizeof(double *));
for (j=0; j<numnodes_y; ++j) {
massest[i][j] = calloc(numnodes_z, sizeof(double));
}
}
krig();
return 0;
}
/* Start Kriging functions */
int krig(void) {
int i, j, entno, ll;
int ac_count=0, o, p, q, m, n;
int chk=1, chkdcmp=1;
double maxx, maxy, maxz, minx, miny, minz, largest=0., inc, ii, jj, kk;
double xinc, yinc, zinc, dist;
double *xgrid, *ygrid, *zgrid;
/* double var[10][10][10];
double upper95[10][10][10];
double lower95[10][10][10]; */
double *ac_wts, *ac_est;
double sumwts=0., gl_covar=0., gvar=0., t_smpcovar=0., glovar, var_y;
double svar[3], disp=0.;
double dx, dy, dz, xcen, ycen, zcen;
FILE *ifp, *ofp;
svar[0]=624.383;
svar[1]=1088.323;
svar[2]=10.0618;
maxx=maxy=maxz=0.;
minx=miny=minz=(double)1.0e12;
ifp=fopen("bchemfile.u","r");
ofp=fopen("test.out","w");
xgrid = calloc(numnodes_x, sizeof(double));
ygrid = calloc(numnodes_y, sizeof(double));
zgrid = calloc(numnodes_z, sizeof(double));
ac_wts = calloc(ndata, sizeof(double));
ac_est = calloc(ndata, sizeof(double));
for (i=0; fscanf(ifp, "%d",&entno) != EOF; ++i) {
fscanf(ifp,"%lf %lf %lf %lf\n",&east[i],&north[i],&elev[i],&bchem[i]);
if (east[i] > maxx) maxx=east[i];
if (east[i] < minx) minx=east[i];
if (north[i] > maxy) maxy=north[i];
if (north[i] < miny) miny=north[i];
if (elev[i] > maxz) maxz=elev[i];
if (elev[i] < minz) minz=elev[i];
}
largest=sqrt( pow(maxx-minx,2.) + pow(maxy-miny,2.) + pow(maxz-minz,2.) );
inc = largest/10000.;
var=svar[1];
for (i=0; i<10000; ++i) {
if (-(svar[1]-sinehole((double)i*inc,svar,3)) > disp)
disp = -(svar[1]-sinehole((double)i*inc,svar,3));
}
var += disp;
printf("%lf\n",var);
xinc=(double)(maxx-minx)/(numnodes_x-1.);
yinc=(double)(maxy-miny)/(numnodes_y-1.);
zinc=(double)(maxz-minz)/(numnodes_z-1.);
xgrid[0]=minx;
ygrid[0]=miny;
zgrid[0]=minz;
for (i=1; i<numnodes_x; ++i) xgrid[i]=xgrid[i-1]+xinc;
for (i=1; i<numnodes_y; ++i) ygrid[i]=ygrid[i-1]+yinc;
for (i=1; i<numnodes_z; ++i) zgrid[i]=zgrid[i-1]+zinc;
dx=fabs(xgrid[2]-xgrid[1]);
dy=fabs(ygrid[2]-ygrid[1]);
dz=fabs(zgrid[2]-zgrid[1]);
/* Calculate sample covariance matrix */
for (o=0; o<numnodes_x; ++o) {
for (p=0; p<numnodes_y; ++p) {
for (q=0; q<numnodes_z; ++q) {
if (chkdcmp == 1) {
for (m=0; m<ndata; ++m) {
for (n=m+1; n<ndata; ++n) {
ii=fabs(east[m]-east[n]);
jj=fabs(north[m]-north[n]);
kk=fabs(elev[m]-elev[n]);
smpdist[m][n]=sqrt(pow(ii,2.)+pow(jj,2.)+pow(kk,2.));
smpcovar[m][n]=var-sinehole(smpdist[m][n],svar,3);
a[m][n]=var-sinehole(smpdist[m][n],svar,3);
if (smpcovar[m][n] < 0.)
printf("COV not bounded by (VAR =%lf) at smpcovar\n",var);
smpcovar[n][m]=smpcovar[m][n];
}
}
/* To guarantee unbiasedness in the estimate, the kriging wts are
constrained such that the sum of the wts is unity. This
constraint is built into the linear system: */
for (n=0; n<ndata; ++n) {
smpcovar[ndata][n]=1.;
smpcovar[n][ndata]=1.;
a[ndata][n]=1.;
a[n][ndata]=1.;
/* The model of Journel and Huijbregts (1978) requires that
smpcovar[0]=sill. It allows exact interpolation, should
an estimate node fall exactly on a sample location: */
smpcovar[n][n]=var;
a[n][n]=var;
}
smpcovar[ndata][ndata]=0.;
a[ndata][ndata]=0.;
/* Replace the sample covariance matrix with it's L-U permutation
so it can be used in the LUBKSB solver */
chkdcmp=ludcmp();
}
/* Create the vector, estcovar, that quantifies the spatial dependence
between estimation block and the sample points that will be weighted */
xcen=xgrid[o];
ycen=ygrid[p];
zcen=zgrid[q];
for (i=0; i<ndata; ++i) {
estcovar[i]=g_quad(east[i], north[i], elev[i], xcen, ycen,
zcen, dx, dy, dz, svar);
ac_est[i]+=estcovar[i];
}
/* Insert 1 in the last position to ensure unbiasedness */
estcovar[ndata]=1.;
/* The linear system has the form [smpcovar]<wts>=<estcovar>.
Solve for the kriging weights using ludcmp() and lubksb() */
b=estcovar;
chk=lubksb();
chk=negwts();
var_y=0.;
for (i=0; i<ndata; ++i) {
for (j=0; j<ndata; ++j) {
var_y+=b[i]*b[j]*smpcovar[i][j];
if (var_y < 0.) printf("POSITIVE DEFINITE CONDITION VIOLATED\n");
}
}
for (i=0; i<ndata; ++i) {
massest[o][p][q]+=b[i]*bchem[i];
ac_wts[i]+=b[i];
}
printf("%d %d %d %lf\n",o,p,q,massest[o][p][q]);
}
}
}
/* Carry out the global uncertainty computations */
for (o=0; o<numnodes_x; ++o) {
for (p=0; p<numnodes_y; ++p) {
for (q=0; q<numnodes_z; ++q) {
for (i=o+1; i<numnodes_x-1; ++i) {
for (j=p+1; j<numnodes_y-1; ++j) {
for (m=q+1; m<numnodes_z-1; ++m) {
ii=fabs(xgrid[o]-xgrid[i]);
jj=fabs(ygrid[p]-ygrid[j]);
kk=fabs(zgrid[q]-zgrid[m]);
dist=sqrt(pow(ii,2.)+pow(jj,2.)+pow(kk,2.));
gl_covar+=var-sinehole(dist,svar,3);
ac_count+=1;
}
}
}
}
}
}
gl_covar /= (double)ac_count;
for (ll=0; ll<ndata; ll++) {
/* fprintf(ofp,"%lf %lf %lf %lf\n",east[ll],north[ll],elev[ll],bchem[ll]);*/
ac_wts[ll] /= ( ((double)numnodes_x)*((double)numnodes_y)*((double)numnodes_z));
sumwts += ac_wts[ll];
ac_est[ll] /= ( ((double)numnodes_x)*((double)numnodes_y)*((double)numnodes_z));
gvar += ac_wts[ll]*ac_est[ll];
fprintf(ofp,"%d %lf %lf\n", ll, ac_wts[ll], ac_est[ll]);
}
for (i=1; i<ndata; ++i) {
for (j=1; j<ndata; ++j) {
t_smpcovar += ac_wts[i]*ac_wts[j]*smpcovar[i][j];
fprintf(ofp,"y%lf %lf\n",smpcovar[i][j],t_smpcovar);
}
}
glovar = gl_covar + t_smpcovar - 2.*gvar;
printf("sumwts : %lf\n\n",sumwts);
printf("gl_covar : %lf\n",gl_covar);
printf("t_smpcovar: %lf\n",t_smpcovar);
printf("gvar : %lf\n",gvar);
printf("----------------------\n");
printf("glovar : %lf\n",glovar);
fclose(ifp);
fclose(ofp);
/* free(indx);
free(east);
free(north);
free(elev);
free(bchem);
free(estdist);
free(estcovar);
free(vv);
free(smpdist);
free(smpcovar);
free(a);*/
return 0;
}
double sinehole(double h, double svar[3], int nparam)
{
double gam;
if (nparam != 3) printf("incorrect number of parameters for sinehole");
if (h != 0) gam=svar[0]+(svar[1]-svar[0])*(1.-sin(4.4934*h/svar[2])/(4.4934*h/svar[2]));
else gam=svar[0];
return gam;
}
/* This function is a modified version of the Fortran program
presented in Carr and Palmer (1993) */
double g_quad(double x, double y, double z,
double xcen, double ycen, double zcen,
double dx, double dy, double dz, double svar[3])
{
int i, j, k, n;
double tj[4], wj[4], tk[4], wk[4], tm[4], wm[4], t[4], w[4], somg, dist;
n=4;
t[0] = -.8611363116;
t[3] = -t[0];
t[1] = -.3399810436;
t[2] = -t[1];
w[0] = w[3] = .3478548452;
w[1] = w[2] = .6521451548;
for (i=0; i<n; ++i) {
tj[i]=ycen+t[i]*dy/2.;
tk[i]=xcen+t[i]*dx/2.;
tm[i]=zcen+t[i]*dz/2.;
wj[i]=w[i]*dy/2.;
wk[i]=w[i]*dx/2.;
wm[i]=w[i]*dz/2.;
}
somg=0.;
for (i=0; i<n; ++i) {
for (j=0; j<n; ++j) {
for (k=0; k<n; ++k) {
dist=sqrt( pow(x-tk[i],2.) + pow(y-tj[j],2.) + pow(z-tm[k],2.) );
somg=somg+(var-sinehole(dist,svar,3))*wk[i]*wj[j]*wm[k];
}
}
}
return somg/(dx*dy*dz);
}
/* This function is a modified version of the program
presented in Press et al. (1996) */
#define tiny 1.0e-20
ludcmp()
{
int i, j, k, imax=0, n, np;
double big, dum, sum, tmp;
n=np=ndata+1;
for (i=0; i<n; i++) {
big=0.;
for (j=0; j<n; j++) {
tmp=fabs(a[i][j]);
if (tmp > big) big=tmp;
}
if (big == 0.) printf("\n\nSingular matrix at ludcmp()\n");
vv[i]=1./big;
}
for (j=0; j<n; j++) {
for (i=0; i<j; i++) {
sum=a[i][j];
for (k=0; k<i; k++) sum -= a[i][k]*a[k][j];
a[i][j]=sum;
}
big=0.;
for (i=j; i<n; i++) {
sum=a[i][j];
for (k=0; k<j; k++)
sum -= a[i][k]*a[k][j];
a[i][j]=sum;
dum=vv[i]*fabs(sum);
if (dum >= big) {
imax=i;
big=dum;
}
}
if (j != imax) {
for (k=0; k<n; k++) {
dum=a[imax][k];
a[imax][k]=a[j][k];
a[j][k]=dum;
}
vv[imax]=vv[j];
}
indx[j]=imax;
if (a[j][j] == 0.) a[j][j]=tiny;
if (j != n-1) {
dum=1./a[j][j];
for (i=j+1; i<n; i++) a[i][j] *= dum;
}
}
/* for (i=0; i<n; ++i) {
for (j=0; j<n; ++j) {
printf("%lf\n",a[i][j]);
}
}*/
return 0;
}
/* This function is a modified version of the program
presented in Press et al. (1996) */
lubksb()
{
int i,ii=0,ll,j,n,id=0;
double sum;
n=ndata+1;
for (i=0; i<n; i++) {
ll=indx[i];
sum=b[ll];
b[ll]=b[i];
if (id) {
for (j=ii; j<=i-1; j++) sum -= (a[i][j])*(b[j]);
}
else if (sum) {
ii=i;
id=1;
}
b[i]=sum;
}
for (i=n-1; i>=0; i--) {
sum=b[i];
for (j=i+1;j<n;j++) sum -= a[i][j]*b[j];
b[i]=sum/a[i][i];
/*if (a[i][i] < 1.e-20 && a[i][i] > -1.e-20) printf("%lf\n",b[i]);*/
}
return 0;
}
/* This function is a modified version of the Fortran program
presented in Deutsch (1996) */
negwts()
{
int i,nneg=0;
double avgneg=0.,avgcov=0.,sumwts;
for (i=0; i<ndata; ++i) {
if (b[i] < 0.) {
nneg+=1;
avgneg-=b[i];
avgcov+=estcovar[i];
}
}
if (nneg > 1) {
avgneg/=(double)nneg;
avgcov/=(double)nneg;
sumwts=0.;
for (i=0; i<ndata; ++i) {
if (b[i] < avgneg && estcovar[i] < avgcov) b[i]=0.;
if (b[i] < 0.) b[i]=0;
sumwts+=b[i];
}
if (sumwts > 0.) {
for (i=0; i<ndata; ++i) b[i]/=sumwts;
}
}
return 0;
}