https://github.com/cran/gstat
Tip revision: c5921c572113e39c365daf8757d324c5cfa49790 authored by Edzer J. Pebesma on 12 October 2006, 00:00:00 UTC
version 0.9-35
version 0.9-35
Tip revision: c5921c5
block.c
/*
Gstat, a program for geostatistical modelling, prediction and simulation
Copyright 1992, 2003 (C) Edzer J. Pebesma
Edzer J. Pebesma, e.pebesma@geo.uu.nl
Department of physical geography, Utrecht University
P.O. Box 80.115, 3508 TC Utrecht, The Netherlands
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version. As a special exception, linking
this program with the Qt library is permitted.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
(read also the files COPYING and Copyright)
*/
/*
* block.c: calculate block discretization at a prediction location
*/
#include <stdio.h>
#include <math.h>
#include "defs.h"
#include "userio.h"
#include "debug.h"
#include "data.h"
#include "utils.h"
#include "glvars.h"
#include "block.h"
/* Gaussian quadrature adapted from Carr & Palmer, MG 25(5), p. 514:
gauss[4]: -0.8611363116, -0.3399810436, 0.3399810436, 0.8611363116
gauss_w[4]: 0.3478548452, 0.6521451548, 0.6521451548, 0.3478548452
both divided by 2
*/
static double gauss[4] = {
-0.4305681558,
-0.1699905218,
0.1699905218,
0.4305681558
};
double gauss_w[4] = {
0.1739274226,
0.3260725774,
0.3260725774,
0.1739274226
};
int reset_anyway = 0;
void reset_block_discr(void) {
reset_anyway = 1;
}
DATA *block_discr(DATA *d, const DPOINT *block, const DPOINT *where) {
/*
* Purpose : get locations that discretize the estimation area
* Created by : Edzer J. Pebesma
* Date : april 6th, 1992
* Modified : jan 13, 1993, june 3, 1993, june 23 1993
* Prerequisites : where;
* the contents of block (if not NULL) may not change
* during the run of the program: only once the full
* discretisation is calculated.
* Returns : DATA list with 1 point (x,y,z) or a list;
* Side effects : dynamic memory (re)allocation for d;
* uses IS_BLOCK(where) to determine whether point or block is needed;
* discretisations is rectangular using dimensions of *block;
*
* only the x, y and z field are used for locations,
* weights are put in the weight union.
*/
int i, k, l, m, ndim;
int ndiscr = 0, restart = 0;
/*
* if restart != 0, block descretizations are re-calculated.
*/
/*
* ndim: number of dimensions (x, y, z, ..?) in data
* ndiscr: total number of points discretizing the block in this run
* lastndiscr: number of discretizing points in last call
* gl_nblockdiscr: number of points in each dimension in discretizing grid
*/
static double *dx = NULL, *dy = NULL, *dz = NULL;
static double *weight = NULL;
static int max = 0;
DATA *area = NULL, **data;
/*
* d has the points of the discretization of the estimation area;
* dx, dy, dz and weight are remembered in case blockdiscretization has not
* changed
*/
if (IS_POINT(where) && d != NULL) { /* point shortcut */
d->list[0]->x = where->x;
d->list[0]->y = where->y;
d->list[0]->z = where->z;
d->list[0]->u.weight = 1.0; /* weight */
d->n_list = 1;
if (DEBUG_BLOCK) {
printlog("block discretization (dist is weight):\n");
d->mode = X_BIT_SET | Y_BIT_SET | Z_BIT_SET;
print_data_list(d);
}
return d;
}
/*
* calculate the neccesary number of block discrimination points:
*/
ndim = 0;
ndiscr = 1;
if (IS_BLOCK(where)) {
if ((area = get_data_area()) == NULL) {
if (block->x > 0.0) ndim++;
if (block->y > 0.0) ndim++;
if (block->z > 0.0) ndim++;
for (i = 0; i < ndim; i++)
ndiscr *= gl_nblockdiscr;
} else
ndiscr = area->n_list;
}
if (area == NULL && ndiscr > max) {
dx = (double *) erealloc(dx, ndiscr * sizeof(double));
dy = (double *) erealloc(dy, ndiscr * sizeof(double));
dz = (double *) erealloc(dz, ndiscr * sizeof(double));
weight = (double *) erealloc(weight, ndiscr * sizeof(double));
max = ndiscr;
}
/*
* (re)allocate the memory neccesary:
*/
if (d == NULL) { /* allocate ndiscr: */
d = (DATA *) emalloc (sizeof(DATA));
init_one_data(d);
assert(get_n_vars() > 0);
data = get_gstat_data();
d->pp_norm2 = data[0]->pp_norm2; /* assign distance function */
d->what_is_u = U_ISWEIGHT;
d->n_X = 0;
d->colnx = d->colny = d->colnz = 1;
d->list = (DPOINT **) emalloc(ndiscr * sizeof(DPOINT *));
for (i = 0; i < ndiscr; i++)
d->list[i] = (DPOINT *) emalloc (sizeof(DPOINT));
d->n_max = d->n_list = ndiscr;
restart = 1;
} else if (ndiscr > d->n_max) { /* resize if n_max < ndiscr: */
d->list = (DPOINT **) erealloc(d->list, ndiscr * sizeof(DPOINT *));
for (i = d->n_max; i < ndiscr; i++)
d->list[i] = (DPOINT *) emalloc (sizeof(DPOINT));
d->n_max = d->n_list = ndiscr;
restart = 1;
} else if (reset_anyway) {
reset_anyway = 0;
restart = 1;
}
if (restart && ndim > 0 && area == NULL) {
/* set up block regular or Gaussian block discretisation */
i = 0;
switch (ndim) {
case 1:
if (block->y > 0 || block->z > 0)
ErrMsg(ER_IMPOSVAL,
"block_discr(): block y and z dimensions must be 0");
if (gl_gauss) {
for (k = 0; k < 4; k++) {
dx[i] = block->x * gauss[k];
dy[i] = dz[i] = 0.0;
weight[i] = gauss_w[k];
i++;
}
} else {
for (k = 0; k < gl_nblockdiscr; k++) {
dx[i] = block->x *
(-0.5 + (1.0/gl_nblockdiscr) * (0.5 + k));
dy[i] = dz[i] = 0.0;
weight[i] = 1.0 / (1.0 * ndiscr);
i++;
}
}
break;
case 2:
if (block->z > 0)
ErrMsg(ER_IMPOSVAL, "block_discr(): 2D block->z must be zero");
if (gl_gauss) {
for (k = 0; k < 4; k++) {
for (l = 0; l < 4; l++) {
dx[i] = block->x * gauss[k];
dy[i] = block->y * gauss[l];
weight[i] = gauss_w[k] * gauss_w[l];
dz[i] = 0.0;
i++;
} /* for l */
} /* for k */
} else {
for (k = 0; k < gl_nblockdiscr; k++) {
for (l = 0; l < gl_nblockdiscr; l++) {
dx[i] = block->x *
(-0.5 + (1.0/gl_nblockdiscr) * (0.5 + k));
dy[i] = block->y *
(-0.5 + (1.0/gl_nblockdiscr) * (0.5 + l));
dz[i] = 0.0;
weight[i] = 1.0 / (1.0 * ndiscr);
i++;
} /* for l */
} /* for k */
}
break;
case 3:
if (gl_gauss) {
for (k = 0; k < 4; k++) {
for (l = 0; l < 4; l++) {
for (m = 0; m < 4; m++) {
dx[i] = block->x * gauss[k];
dy[i] = block->y * gauss[l];
dz[i] = block->z * gauss[m];
weight[i] = gauss_w[k] *
gauss_w[l] * gauss_w[m];
i++;
} /* for m */
} /* for l */
} /* for k */
} else {
for (k = 0; k < gl_nblockdiscr; k++) {
for (l = 0; l < gl_nblockdiscr; l++) {
for (m = 0; m < gl_nblockdiscr; m++) {
dx[i] = block->x *
(-0.5+(1.0/gl_nblockdiscr) * (0.5 + k));
dy[i] = block->y *
(-0.5+(1.0/gl_nblockdiscr) * (0.5 + l));
dz[i] = block->z *
(-0.5+(1.0/gl_nblockdiscr) * (0.5 + m));
weight[i] = 1.0 / (1.0 * ndiscr);
i++;
} /* for m */
} /* for l */
} /* for k */
}
break;
default:
ErrMsg(ER_IMPOSVAL, "block_discr()");
} /* switch */
if (i != ndiscr) /* silly check ! */
ErrMsg(ER_IMPOSVAL, "block_discr()");
} /* end if restart */
if (IS_BLOCK(where)) {
if (area != NULL) {
for (i = 0; i < ndiscr; i++) {
d->list[i]->x = area->list[i]->x + where->x;
d->list[i]->y = area->list[i]->y + where->y;
d->list[i]->z = area->list[i]->z + where->z;
if (area->colnvariance)
d->list[i]->u.weight = area->list[i]->variance;
else
d->list[i]->u.weight = 1.0 / ndiscr;
}
} else {
for (i = 0; i < ndiscr; i++) {
d->list[i]->x = where->x - dx[i];
d->list[i]->y = where->y - dy[i];
d->list[i]->z = where->z - dz[i];
d->list[i]->u.weight = weight[i];
}
}
} else { /* where is point */
d->list[0]->x = where->x;
d->list[0]->y = where->y;
d->list[0]->z = where->z;
d->list[0]->u.weight = 1.0; /* weight */
}
d->n_list = ndiscr;
if (DEBUG_BLOCK) {
printlog("block discretization (dist is weight):\n");
d->mode = X_BIT_SET | Y_BIT_SET | Z_BIT_SET;
print_data_list(d);
}
return d;
}