https://github.com/cran/RandomFields
Tip revision: 25ed5157dd19e7ca108a1d5ed2581a6a98c376d5 authored by Martin Schlather on 23 February 2011, 00:00:00 UTC
version 2.0.43
version 2.0.43
Tip revision: 25ed515
markov.cc
/*
Authors
Martin Schlather, martin.schlather@math.uni-goettingen.de
wrapper for Havard Rue's GMRF library
Copyright (C) 2001 -- 2011 Martin Schlather,
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.
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., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
*/
// was passiert wenn geo-coef nicht verfuegbar?
// wie wird *problem geloescht ??
#include <math.h>
#include <stdio.h>
#include <assert.h>
#include "RF.h"
#ifdef INCLUDEMARKOV
#include "GMRFLib/GMRFLib.h"
#endif
#define MARKOV_DEBUG false
// #define MARKOV_DEBUG true
typedef struct markov_storage {
// typedef cannot be moved to RF.h since
// it requires GMRFLib/GMRFLib.h and the inclusion
// in RF.h leads to conflicts with blas.h definitions
// used in circulant.cc
int rowcol[2];
#ifdef INCLUDEMARKOV
GMRFLib_problem_tp *problem;
#endif
} markov_storage;
void Markov_destruct(void **S)
{
if (*S!=NULL) {
markov_storage *s;
s = *((markov_storage**)S);
#ifdef INCLUDEMARKOV
GMRFLib_free_problem(s->problem);
#endif
free(*S);
*S = NULL;
}
}
void SetParamMarkov(int *action,int *neighbours, double *precision,
int *cyclic, int *maxmem) {
markov_param *lp = &(GLOBAL.markov);
if (*action) {
lp->neighbours = *neighbours;
if (lp->neighbours < 2) {
if (PL>0) PRINTF("minimal neighbourhood is 2");
lp->neighbours = 2;
}
if (lp->neighbours > 3) {
if (PL>0) PRINTF("maximal neighbourhood is 3");
lp->neighbours = 3;
}
lp->precision = *precision;
lp->cyclic = *cyclic;
lp->maxmem = *maxmem;
} else {
*neighbours = lp->neighbours;
*precision = lp->precision;
*cyclic = lp->cyclic;
*maxmem = lp->maxmem;
}
}
int init_markov(method_type *meth){
SET_DESTRUCT(Markov_destruct);
int err;
#ifdef INCLUDEMARKOV
int name, endfor, v, d;
double param, range, gridstep=0.0, scale=0.0;
markov_storage* s;
GMRFLib_geo_problem_tp *geo_problem = NULL;
int reduceddim reduceddim = meth->reduceddim;
cov_model *cov=meth->cov;
cov_fct *cf = &(CovList[cov->nr]);
if (MARKOV_DEBUG) {
FILE *fp;
if ((fp = fopen("available-geo-coefs.dat", "w")))
{
PRINTF("print available geo-coefs to file\n");
GMRFLib_print_geo_coefs(fp);
fclose(fp);
}
}
SET_DESTRUCT(Markov_destruct);
if ((meth->S = malloc(sizeof(markov_storage)))==0){
err=ERRORMEMORYALLOCATION; goto ErrorHandling;
}
s = (markov_storage*) meth->S;
s->problem = NULL;
s->rowcol[1] = 1; // falls nur 1-dimensionales Feld
//////////////////////////////////////////////////////////////////////
err = ERRRORCURRENTLY;
goto ErrorHandling;
//////////////////////////////////////////////////////////////////////
if (reduceddim > 2) {err=ERRORWRONGDIM; goto ErrorHandling;}
if (loc->totalpoints > lp->maxmem) {
sprintf(ERRORSTRING_OK, "%d", lp->maxmem);
sprintf(ERRORSTRING_WRONG,"%ld", loc->totalpoints);
err = ERRORMARKOVMAXMEMORY;
goto ErrorHandling;
}
err = NOERROR;
if (cf->implemented[Markov] != IMPLEMENTED) err=ERRORNOTDEFINED;
if (!cov->diag) err=ERRORONLYISOTROPIC;
if (err!=NOERROR) goto ErrorHandling;
// current restriction:
gridstep = meth->xorig[0][XSTEP];
assert(gridstep != 0.0);
for (d=1; d<cov->tsdim; d++) {
if (meth->xorig[d][XSTEP] != gridstep) {
err=ERRORONLYREGULARGRID; goto ErrorHandling;
}
}
double *aniso;
aniso = getAnisoMatrix(meth);
if (meth->type == TypeIso && aniso[0] != 0.0) {
scale = 1.0 / aniso[0];
free(aniso);
} else {
free(aniso);
err=ERRORONLYISOTROPIC; goto ErrorHandling;
}
param = 1.0; // dummy
name = GMRFLib_CORTP_EXP; // exponential
range = 3.0;
if (kc->nr == GAUSS) {
name = GMRFLib_CORTP_GAUSS;
range = sqrt(3.0);
} else if (kc->nr == WHITTLE) {
#define ONEMINUSEPS 1.0 - 1e-10
if (gp->general.naturalscaling) {
GetNaturalScaling(cov, &range, &err);
error("markov not programmed anymore");
if (err != NOERROR) goto ErrorHandling;
range = 1.0 / range;
} else {
double xmax, x, xmin, value;
name = GMRFLib_CORTP, MATERN;
xmax = 1.0;
WM(&xmax, cov->p[0][0], &value);
while (value > 0.05) {
xmax *= 2.0;
WM(&xmax, cov->p[0][0], &value);
}
xmin = (xmax > 1.0) ? 0.5 * xmax : 0.0;
while (ONEMINUSEPS > xmin / xmax) {
x = 0.5 * (xmin + xmax);
WM(&x, cov->p[0][0], &value);
if (value > 0.05) xmin = x;
else xmax = x;
}
range = 1.0 / x;
}
} else assert(kc->nr == EXPONENTIAL);
range *= scale / gridstep;
v = 0;
for (d=0; d<cov->tsdim; d++) {
assert(v <= 1);
s->rowcol[v++] = loc->length[d];
}
GMRFLib_set_error_handler_off();
if (GMRFLib_is_geo_coefs(name, lp->neighbours, param, range)
!= GMRFLib_SUCCESS) {
sprintf(ERRORSTRING_WRONG,
"GMRFLib_is_geo_coefs: err=#%d; param=%f, range=%f",
err, param, range);
err = ERRORMARKOVPARAMETER;
goto ErrorHandling;
}
err = GMRFLib_init_geo_problem(&geo_problem, name,
lp->neighbours, param, range,
s->rowcol[0], s->rowcol[1],
&lp->precision, lp->cyclic);
if (err != GMRFLib_SUCCESS) {
sprintf(ERRORSTRING_WRONG,
"GMRFLib_init_geo_problem: err=#%d; param=%f, range=%f",
err, param, range);
err = ERRORMARKOVPARAMETER;
goto ErrorHandling;
}
GetRNGstate();
GMRFLib_uniform_init((unsigned int) (4294967296.0 * UNIFORM_RANDOM));
PutRNGstate();
err = GMRFLib_init_problem(&(s->problem), NULL, NULL, NULL, NULL,
geo_problem->graph,
geo_problem->Qfunc, geo_problem->Qfunc_arg,
NULL, NULL, GMRFLib_NEW_PROBLEM);
if (err != GMRFLib_SUCCESS) {
sprintf(ERRORSTRING_WRONG,
"GMRFLib_init_problem: err=#%d; param=%f, range=%f",
err, param, range);
err = ERRORMARKOVPARAMETER;
goto ErrorHandling;
}
err = NOERROR;
if (MARKOV_DEBUG) {
FILE *fp;
int i, j;
/*
print the Q-matrix
*/
fp = fopen("Q.dat", "w");
PRINTF("write [Q.dat]\n");
for(i = 0; i < geo_problem->graph->n; i++)
{
fprintf(fp, "%d %d %.16f\n", i, i,
geo_problem->Qfunc(i, i, geo_problem->Qfunc_arg));
for(j = 0; j < geo_problem->graph->nnbs[i]; j++)
{
fprintf(fp, "%d %d %.16f\n", i,
geo_problem->graph->nbs[i][j],
geo_problem -> Qfunc(i,
geo_problem->graph->nbs[i][j],
geo_problem->Qfunc_arg));
}
}
fclose(fp);
}
// if (noerror_repeat) return NOERROR_REPEAT;
// else return err;
if (geo_problem != NULL) GMRFLib_free_geo_problem(geo_problem);
return NOERROR_REPEAT;
#else
err = ERRORMARKOVNOTINCLUDED; goto ErrorHandling;
#endif
ErrorHandling:
#ifdef INCLUDEMARKOV
if (geo_problem != NULL) GMRFLib_free_geo_problem(geo_problem);
#endif
return err;
}
void do_markov(method_type *meth, res_type *res ) {
#ifdef INCLUDEMARKOV
markov_storage *S;
long k, i, j;
int ncol, nrow, idx;
GMRFLib_problem_tp *problem;
S = (markov_storage*) meth->S;
problem = S->problem;
nrow = S->rowcol[0];
ncol = S->rowcol[1];
GMRFLib_sample(problem);
k = 0;
for(j=0; j<ncol; j++) {
for(i=0; i<nrow; i++, k++) {
GMRFLib_lattice2node(&idx, i, j, nrow, ncol);
res[k] += (res_type) problem->sample[idx];
}
}
if (MARKOV_DEBUG) {//////////////////////
FILE *fp;
double *marg_var;
fp = fopen("sample.dat", "w");
PRINTF("write [sample.dat]\n");
for(i=0; i<nrow; i++)
{
for(j=0; j<ncol; j++)
{
GMRFLib_lattice2node(&idx, i, j, nrow, ncol);
fprintf(fp, " %lf", problem->sample[idx]);
}
fprintf(fp, "\n");
}
fclose(fp);
GMRFLib_Qinv(problem, GMRFLib_QINV_ALL);
fp = fopen("variances.dat", "w");
PRINTF("write [variances.dat]\n");
for(i=0; i<nrow; i++)
{
for(j=0; j<ncol; j++)
{
GMRFLib_lattice2node(&idx, i, j, nrow, ncol);
marg_var = GMRFLib_Qinv_get(problem, idx, idx);
fprintf(fp, " %lf", *marg_var);
}
fprintf(fp, "\n");
}
fclose(fp);
} //////////////////////////////
#else
assert(false);
#endif
}
// ruwe ??