primitive.others.cc
/*
Authors
Martin Schlather, schlather@math.uni-mannheim.de
Definition of other simple function
Copyright (C) 2017 -- 2018 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 3
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.
*/
#include <Rmath.h>
#include <R_ext/Lapack.h>
#include <R_ext/Linpack.h>
#include "questions.h"
#include "RF.h"
#include "primitive.others.h"
#include "Coordinate_systems.h"
#include "operator.h"
#include "QMath.h"
//#include "AutoRandomFields.h"
//#include "shape.h"
//#include "Processes.h"
// Sigma(x) = diag>0 + A'xx'A
void kappa_Angle(int i, model *cov, int *nr, int *nc){
int dim = OWNXDIM(0);
*nr = i == ANGLE_DIAG ? dim : 1;
*nc = i <= ANGLE_DIAG && (dim==3 || i != ANGLE_LATANGLE) &&
(dim==2 || i != ANGLE_RATIO) ? 1 : -1;
}
void AngleMatrix(model *cov, double *A) {
double c, s,
*diag = P(ANGLE_DIAG);
if (GLOBAL.coords.anglemode == radians) {
c = COS(P0(ANGLE_ANGLE));
s = SIN(P0(ANGLE_ANGLE));
} else {
assert(GLOBAL.coords.anglemode == degree);
c = COS(P0(ANGLE_ANGLE) * piD180);
s = SIN(P0(ANGLE_ANGLE) * piD180);
}
int
dim = OWNXDIM(0) ;
assert(dim ==2 || dim == 3);
if (dim == 2) {
A[0] = c;
A[1] = s;
A[2] = -s;
A[3] = c;
} else {
double pc, ps;
if (GLOBAL.coords.anglemode == radians) {
pc = COS(P0(ANGLE_LATANGLE));
ps = SIN(P0(ANGLE_LATANGLE));
} else {
pc = COS(P0(ANGLE_LATANGLE) * piD180);
ps = SIN(P0(ANGLE_LATANGLE) * piD180);
}
/*
c -s 0 pc 0 -ps c*pc -s -c*ps
s c 0 * 0 1 0 = s*pc c -s*ps
0 0 1 ps 0 pc ps 0 pc
*/
A[0] = c * pc;
A[1] = s * pc;
A[2] = ps;
A[3] = -s;
A[4] = c;
A[5] = 0.0;
A[6] = -c * ps;
A[7] = -s * ps;
A[8] = pc;
}
if (diag!= NULL) {
int k,d,j;
for (k=d=0; d<dim; d++) {
for (j=0; j<dim; j++) {
A[k++] *= diag[j];
}
}
} else {
double ratio = P0(ANGLE_RATIO);
A[1] /= ratio;
A[3] /= ratio;
}
}
void Angle(double *x, model *cov, double *v) { /// to do: extend to 3D!
double A[9];
int
dim = OWNXDIM(0) ;
AngleMatrix(cov, A);
Ax(A, x, dim, dim, v);
}
void invAngle(double *x, model *cov, double *v) { /// to do: extend to 3D!
double A[9],
c = COS(P0(ANGLE_ANGLE)),
s = SIN(P0(ANGLE_ANGLE)),
*diag = P(ANGLE_DIAG);
int d,
dim = OWNXDIM(0) ;
bool inf = x[0] == RF_INF;
for (d=1; d<dim; d++) inf &= x[d] == RF_INF;
if (inf) {
for (d=0; d<dim; v[d++] = RF_INF);
return;
}
bool neginf = x[0] == RF_NEGINF;
for (d=1; d<dim; d++) neginf &= x[d] == RF_NEGINF;
if (neginf) {
for (d=0; d<dim; v[d++] = RF_NEGINF);
return;
}
if (dim == 2) {
A[0] = c;
A[1] = -s;
A[2] = s;
A[3] = c;
} else {
double pc = COS(P0(ANGLE_LATANGLE)),
ps = SIN(P0(ANGLE_LATANGLE));
A[0] = c * pc;
A[1] = -s;
A[2] = -c * ps;
A[3] = s * pc;
A[4] = c;
A[5] = -s * ps;
A[6] = ps;
A[7] = 0.0;
A[8] = pc;
}
if (diag!= NULL) {
int j, k;
for (k=d=0; d<dim; d++) {
for (j=0; j<dim; j++) {
A[k++] /= diag[d];
}
}
} else {
double ratio = P0(ANGLE_RATIO);
A[2] *= ratio;
A[3] *= ratio;
}
Ax(A, x, dim, dim, v);
}
int checkAngle(model *cov){
int dim = OWNXDIM(0);
if (dim != 2 && dim != 3)
SERR1("'%.50s' only works for 2 and 3 dimensions", NICK(cov));
if (PisNULL(ANGLE_DIAG)) {
if (PisNULL(ANGLE_RATIO)) {
SERR2("either '%.50s' or '%.50s' must be given",
KNAME(ANGLE_RATIO), KNAME(ANGLE_DIAG))
} else if (dim != 2) {
SERR1("'%.50s' may be given only if dim=2", KNAME(ANGLE_RATIO))
}
} else {
if (!PisNULL(ANGLE_RATIO))
SERR2("'%.50s' and '%.50s' may not given at the same time",
KNAME(ANGLE_RATIO), KNAME(ANGLE_DIAG));
}
VDIM0 = dim;
VDIM1 = 1;
cov->mpp.maxheights[0] = RF_NA;
cov->matrix_indep_of_x = true;
RETURN_NOERROR;
}
void rangeAngle(model *cov, range_type* range){
model *calling = cov->calling;
assert(calling != NULL);
range->min[ANGLE_ANGLE] = 0.0;
range->max[ANGLE_ANGLE] =
calling->vdim[0] == SCALAR && isDollar(calling) &&
calling->sub[0] != cov && // so a parameter of dollar
SYSTEMS(CALLING) == 1
&& isSymmetric(ISO(CALLING, 0)) ? PI : TWOPI;
range->pmin[ANGLE_ANGLE] = 0.0;
range->pmax[ANGLE_ANGLE] = range->max[ANGLE_ANGLE];
range->openmin[ANGLE_ANGLE] = false;
range->openmax[ANGLE_ANGLE] = true;
range->min[ANGLE_LATANGLE] = 0.0;
range->max[ANGLE_LATANGLE] = PI;
range->pmin[ANGLE_LATANGLE] = 0.0;
range->pmax[ANGLE_LATANGLE] = PI;
range->openmin[ANGLE_LATANGLE] = false;
range->openmax[ANGLE_LATANGLE] = true;
range->min[ANGLE_RATIO] = 0;
range->max[ANGLE_RATIO] = RF_INF;
range->pmin[ANGLE_RATIO] = 1e-5;
range->pmax[ANGLE_RATIO] = 1e5;
range->openmin[ANGLE_RATIO] = false;
range->openmax[ANGLE_RATIO] = true;
range->min[ANGLE_DIAG] = 0;
range->max[ANGLE_DIAG] = RF_INF;
range->pmin[ANGLE_DIAG] = 1e-5;
range->pmax[ANGLE_DIAG] = 1e5;
range->openmin[ANGLE_DIAG] = false;
range->openmax[ANGLE_DIAG] = true;
}
void idcoord(double *x, model *cov, double *v){
int d,
vdim = VDIM0;
for (d=0; d<vdim; d++) v[d]=x[d];
}
int checkidcoord(model *cov){
if (PREVISO(0) != OWNISO(0)) SERR("unequal iso's");
VDIM0 = OWNTOTALXDIM;
VDIM1 = 1;
RETURN_NOERROR;
}
// obsolete??!!
#define NULL_TYPE 0
void NullModel(double VARIABLE_IS_NOT_USED *x,
model VARIABLE_IS_NOT_USED *cov,
double VARIABLE_IS_NOT_USED *v){ return; }
void logNullModel(double VARIABLE_IS_NOT_USED *x,
model VARIABLE_IS_NOT_USED *cov,
double VARIABLE_IS_NOT_USED *v,
int VARIABLE_IS_NOT_USED *Sign){ return; }
Types TypeNullModel(Types required, model *cov,
isotropy_type VARIABLE_IS_NOT_USED i) {
return TypeConsistency(isManifold(required) ? PREVTYPE(0) : required,
(Types) P0INT(NULL_TYPE));
}
void rangeNullModel(model VARIABLE_IS_NOT_USED *cov, range_type *range){
range->min[NULL_TYPE] = TcfType;
range->max[NULL_TYPE] = OtherType;
range->pmin[NULL_TYPE] = TcfType;
range->pmax[NULL_TYPE] = OtherType;
range->openmin[NULL_TYPE] = false;
range->openmax[NULL_TYPE] = false;
}
int checkMath(model *cov){
int i,
variables = DefList[COVNR].kappas,
err = NOERROR;
if (variables > 2) kdefault(cov, variables-1, 1.0);
if (isEarth(OWNISO(0))) {
defn *C = DefList + COVNR;
if ((C->cov==MathCos || C->cov==MathSin || C->cov==MathTan)) {
// printI(cov);
// printf("%.50s\n", TYPE_NAMES[DEFTYPE(0)]);
// APMI0(cov);
SERR1("only radians as angular system allowed for '%.50s'.", NICK(cov));
}
}
for (i=0; i<variables; i++) {
model *sub = cov->kappasub[i];
if (sub != NULL) {
if (i >= 2) SERR("only numerics allowed");
bool plus = DefList[SUBNR].cov == Mathplus ||
DefList[SUBNR].check == checkplus ||
DefList[SUBNR].cov == Mathminus
;
if ((err = CHECK_PASSTF(sub, plus ? OWNTYPE(0) : ShapeType,
// auch falls cov = TrendType ist
1, cov->frame)) != NOERROR){
// if ((err = CHECK(sub, cov->tsdim, cov->xdimown,
/// plus ? cov->typus : ShapeType,
// // auch falls cov = TrendType ist
// OWNDOM(0), OWNISO(0),
// 1, cov->frame)) != NOERROR){
RETURN_ERR(err);
}
if (sub->vdim[0] != 1 || sub->vdim[1] != 1)
SERR("only scalar functions are allowed");
setbackward(cov, sub);
} else if (PisNULL(i)) {
if (i==0 || (DefList[COVNR].cov!=Mathplus &&
DefList[COVNR].cov!=Mathminus &&
DefList[COVNR].cov!=Mathbind
)) SERR("not enough arguments given")
else break;
}
}
cov->pref[Trendproc] = 5;
cov->pref[Direct] = 1;
RETURN_NOERROR;
}
void rangeMath(model VARIABLE_IS_NOT_USED *cov, range_type *range){
int i,
variables = DefList[COVNR].kappas;
set_maxdim(OWN, 0, OWNLOGDIM(0));
assert(MAXDIM(OWN, 0) > 0);
for (i=0; i<variables; i++) {
range->min[i] = RF_NEGINF;
range->max[i] = RF_INF;
range->pmin[i] = -1e5;
range->pmax[i] = 1e5;
range->openmin[i] = true;
range->openmax[i] = true;
}
}
void Mathminus(double *x, model *cov, double *v){
MATH_DEFAULT;
double f = P0(MATH_FACTOR);
if (ISNAN(f)) f = 1.0;
*v = ( (PisNULL(1) && cov->kappasub[1]==NULL) ? -w[0] : w[0]- w[1]) * f;
}
void Mathplus(double *x, model *cov, double *v){
MATH_DEFAULT;
double f = P0(MATH_FACTOR);
if (ISNAN(f)) f = 1.0;
*v = ( (PisNULL(1) && cov->kappasub[1]==NULL)? w[0] : w[0] + w[1]) * f;
}
void Mathdiv(double *x, model *cov, double *v){
MATH_DEFAULT;
double f = P0(MATH_FACTOR);
if (ISNAN(f)) f = 1.0;
*v = (w[0] / w[1]) * f;
}
void Mathmult(double *x, model *cov, double *v){
MATH_DEFAULT;
double f = P0(MATH_FACTOR);
if (ISNAN(f)) f = 1.0;
*v = (w[0] * w[1]) * f;
}
#define IS_IS 1
void Mathis(double *x, model *cov, double *v){
int i,
variables = DefList[COVNR].kappas;
double w[3];
for (i=0; i<variables; i++) {
model *sub = cov->kappasub[i];
if (sub != NULL) {
COV(x, sub, w + i);
} else {
w[i] = i == IS_IS ? P0INT(i) : P0(i);
}
}
switch((int) w[IS_IS]) {
case 0 : *v = (double) (FABS(w[0] - w[2]) <= GLOBAL.nugget.tol); break;
case 1 : *v = (double) (FABS(w[0] - w[2]) > GLOBAL.nugget.tol); break;
case 2 : *v = (double) (w[2] + GLOBAL.nugget.tol >= w[0]); break;
case 3 : *v = (double) (w[2] + GLOBAL.nugget.tol > w[0]); break;
case 4 : *v = (double) (w[0] + GLOBAL.nugget.tol >= w[2]); break;
case 5 : *v = (double) (w[0] + GLOBAL.nugget.tol > w[2]); break;
default : BUG;
}
}
void rangeMathIs(model *cov, range_type *range){
rangeMath(cov, range);
range->min[IS_IS] = 0;
range->max[IS_IS] = 5;
range->pmin[IS_IS] = 0;
range->pmax[IS_IS] = 5;
range->openmin[IS_IS] = false;
range->openmax[IS_IS] = false;
}
void Mathbind(double *x, model *cov, double *v){
int vdim = VDIM0;
MATH_DEFAULT_0(vdim);
double f = P0(DefList[COVNR].kappas - 1);
f = ISNAN(f) ? 1.0 : f;
for (i=0; i<vdim; i++) v[i] = w[i] * f;
}
int check_bind(model *cov) { // R.c
int err;
if ((err = checkMath(cov)) != NOERROR) RETURN_ERR(err);
kdefault(cov, BIND_NCOL, 1);
int
ncol = P0INT(BIND_NCOL),
variables = BIND_VARIABLES;
while (variables > 0 && cov->nrow[variables - 1] == 0 &&
cov->kappasub[variables - 1] == NULL) variables--;
VDIM0 = variables / ncol;
VDIM1 = ncol;
if (VDIM0 * VDIM1 != variables)
SERR1("'%.50s' does not fit the number of components given", KNAME(BIND_NCOL));
cov->ptwise_definite = pt_mismatch;
RETURN_NOERROR;
}
bool allowedDbind(model *cov) {
defn *C = DefList + COVNR;
int i,
kappas = C->kappas;
for (i=0; i<kappas; i++) if (cov->kappasub[i] != NULL) break;
if (i<kappas) {
bool allowed = true,
*D = cov->allowedD;
for (int j=(int) FIRST_DOMAIN; j<=(int) LAST_DOMAINUSER; D[j++] = false);
for (; i<kappas; i++) {
model *k = cov->kappasub[i];
if (k != NULL) {
allowed &= allowedD(k);
for (int j=(int) FIRST_DOMAIN; j<=(int) LAST_DOMAINUSER; j++)
D[j] |= k->allowedD[j];
}
}
//printD(cov);printf("allowd.c=%d\n", allowed);
return allowed;
}
//printf("Allowd.c=%d\n", allowedItrue(cov));
return allowedItrue(cov);
}
bool allowedIbind(model *cov) {
defn *C = DefList + COVNR;
int i,
kappas = C->kappas;
for (i=0; i<kappas; i++) if (cov->kappasub[i] != NULL) break;
if (i<kappas) {
bool allowed = true,
*I = cov->allowedI;
for (int j=(int) FIRST_ISOUSER; j<=(int) LAST_ISOUSER; I[j++] = false);
for (; i<kappas; i++) {
model *k = cov->kappasub[i];
if (k != NULL) {
allowed |= allowedI(k);
for (int j=(int) FIRST_ISOUSER; j<=(int) LAST_ISOUSER; j++)
I[j] &= k->allowedI[j];
}
}
// printI(cov);printf("allowdI.c=%d\n", allowed);
return allowed;
}
return allowedItrue(cov);
}
Types Typebind(Types required, model VARIABLE_IS_NOT_USED *cov,
isotropy_type VARIABLE_IS_NOT_USED iso) {
return required;
/*
defn *C = DefList + COVNR;
int i,
kappas = C->kappas;
for (i=0; i<kappas; i++) if (cov->kappasub[i] != NULL) break;
if (i<kappas) {
// if (is(required))
// printf("bind: requ=%.50s %.50s\n", TYPE_NAMES[required], ISO_NAMES[iso]);
P MI0(cov);
BUG;
return BadType;
} else {
return required;
}
*/
}
void Mathc(double VARIABLE_IS_NOT_USED *x, model *cov, double *v) {
double f = P0(CONST_C);
*v = ISNAN(f) ? 1.0 : f;
}
int check_c(model *cov){
bool negdef = isnowNegDef(cov);
bool tcf = isnowTcf(cov) || equalsSphericalIsotropic(OWNISO(0));
if (negdef) {
if (cov->calling == NULL) BUG;
model *pp = cov->calling->calling;
// the following defines a (positive) constant on the level of a trend
// to be a trend, and not a positive definite function
// For spatially constant covariance functions, see RMconstant
if (pp == NULL ||
(CALLINGNR == PLUS && !isnowNegDef(pp) && !isnowTrend(pp))){
RETURN_ERR(ERRORFAILED); // by definition,
}
}
if (cov->kappasub[0] != NULL) SERR("only numerics allowed");
cov->ptwise_definite =
P0(CONST_C) > 0 ? pt_posdef : P0(CONST_C) < 0 ? pt_negdef : pt_zero;
if (tcf) MEMCOPY(cov->pref, PREF_ALL, sizeof(pref_shorttype));
GLOBAL.internal.warn_mathdef = GLOBAL.internal.warn_mathdef == False ? False
: isNegDef(PREVTYPE(0)) ? True : Nan;
RETURN_NOERROR;
}
void rangec(model *cov, range_type *range){
rangeconstant(cov, range);
}
bool allowedIp(model *cov) {
bool *I = cov->allowedI;
for (int i=(int) FIRST_ISOUSER; i<=(int) LAST_ISOUSER; I[i++] = false);
if (!PisNULL(PROJ_ISO)) {
isotropy_type iso = (isotropy_type) P0INT(PROJ_ISO);
I[iso] = true;
switch(iso) {
case ISOTROPIC: case EARTH_ISOTROPIC: case SPHERICAL_ISOTROPIC :
I[iso] = true;
break;
case DOUBLEISOTROPIC :
case VECTORISOTROPIC :
ERR("'VECTORISOTROPY' not programmed yet");
case SYMMETRIC : case EARTH_SYMMETRIC: case SPHERICAL_SYMMETRIC :
ERR2("Use '%.50s' within arbitrarty mathematical definitions (i.e. in '%.50s') or use the argument 'proj' within definite functions)",
CoordinateSystemOf(iso), NICK(cov));
break;
case CARTESIAN_COORD: case SPHERICAL_COORD: case EARTH_COORD:
I[iso] = true;
break;
case GNOMONIC_PROJ : case ORTHOGRAPHIC_PROJ :
ERR("Do not use projection in 'R,p', but use 'RMtrafo' instead.");
break;
case UNREDUCED :
I[CARTESIAN_COORD] = I[SPHERICAL_COORD] = I[EARTH_COORD] = true;
break;
default : ERR2("'%.50s' not allowed for '%.50s'", ISO_NAMES[iso],KNAME(PROJ_ISO));
}
} else {
I[CARTESIAN_COORD] = I[SPHERICAL_COORD] = I[EARTH_COORD] = true;
}
return false;
}
void proj(double *x, model *cov, double *v){
double f = P0(PROJ_FACTOR);
if (ISNAN(f)) f = 1.0;
int p = P0INT(PROJ_PROJ);
if (p >= 1) {
*v = x[p - 1] * f;
// PMIR(cov); BUG;
// printf("nr=%d f=%10g p=%d x=%10g v=%10g\n", cov->zaehler, f, p, x[p-1], *v);
return;
}
if (p == PROJ_TIME) *v = x[OWNTOTALXDIM - 1] * f;
else if (p == PROJ_SPACE) {
double dummy = 0.0;
int spdim = OWNTOTALXDIM - 1;
for (int d = 0; d<spdim; d++) dummy += x[d] * x[d];
*v = SQRT(dummy) * f;
} else BUG;
}
bool setproj(model *cov){
isotropy_type
iso = PisNULL(PROJ_ISO) ? PREVISO(0) : (isotropy_type) P0INT(PROJ_ISO);
if (!isFixed(iso)) return false; // siehe unten -- nicht immer gebraucht!!
domain_type dom = PREVDOM(0);
Types type = PREVTYPE(0);
bool fixed = isFixed(dom);
if (equalsSpaceIsotropic(iso)) {
NotProgrammedYet("");
if (fixed) {
#if MAXSYSTEMS == 1
set_system(OWN, 0, PREVLOGDIM(0), MAXINT, 2, type, dom, iso);
#else
BUG ?? return false;
set_system(OWN, 0, totlogicaldim - 1, MAXINT, 1, type, dom, ISOTROPIC);
set_system(OWN, 1, 1, 1, 1, OWNTYPE(1), OWNDOM(1), ISOTROPIC);
#endif
} else {
#if MAXSYSTEMS == 1
set_iso(OWN, 0, iso); // rest from the definition
#else
set_iso(OWN, 0, iso);
not programmed yet
#endif
}
} else if (isAnySpherical(iso)) {
int s = 0;
set_system(OWN, s, PREVLOGDIM(s), MAXINT, PREVLOGDIM(s), PREVTYPE(s),
PREVDOM(s), CoordinateSystemOf(PREVISO(s)));
} else if (equalsUnreduced(iso)) {
if (!fixed) return false;
int last = PREVLASTSYSTEM;
for (int s=0; s<=last; s++)
set_system(OWN, s, PREVLOGDIM(s), MAXINT, PREVLOGDIM(s), PREVTYPE(s),
PREVDOM(s), CoordinateSystemOf(PREVISO(s)));
} else {
if (fixed) set_system(OWN, 0, PREVLOGDIM(0), MAXINT, 1, type, dom, iso);
else set_iso(OWN, 0, iso);
}
return true;
}
int checkproj(model *cov){
kdefault(cov, PROJ_FACTOR, 1.0);
kdefault(cov, PROJ_PROJ, 1);
// if (isoown != iso && (iso != UNREDUCED || !equalsCoordinateSystem(isoown))) {
// SERR2("Offered system ('%.50s') does not match the required one ('%.50s')",
// ISO_NAMES[isoown], ISO_NAMES[iso]);
// }
if (P0INT(PROJ_PROJ) < 0 && !GetTime(cov))
SERR2("'%.50s' or '%.50s' used in a context that is not spatio-temporal.",
PROJECTION_NAMES[0], PROJECTION_NAMES[1]);
RETURN_NOERROR;
}
Types Typeproj(Types required, model *cov,
isotropy_type VARIABLE_IS_NOT_USED required_iso) {
if (isBad(TypeConsistency(required, ShapeType)) && // ex-mathdeftype
isBad(TypeConsistency(required, TrendType))) return BadType;
return atleastSpecialised(OWNISO(0), required_iso) ? required : BadType;
}
void rangeproj(model VARIABLE_IS_NOT_USED *cov, range_type *range){
range->min[PROJ_PROJ] = -PROJECTIONS;
range->max[PROJ_PROJ] = OWNXDIM(0);
range->pmin[PROJ_PROJ] = 1;
range->pmax[PROJ_PROJ] = OWNXDIM(0);
range->openmin[PROJ_PROJ] = false;
range->openmax[PROJ_PROJ] = false;
range->min[PROJ_ISO] = FIRST_ISOUSER;
range->max[PROJ_ISO] = LAST_ISOUSER;
range->pmin[PROJ_ISO] = FIRST_ISOUSER;
range->pmax[PROJ_ISO] = LAST_ISOUSER;
range->openmin[PROJ_ISO] = false;
range->openmax[PROJ_ISO] = false;
range->min[PROJ_FACTOR] = RF_NEGINF;
range->max[PROJ_FACTOR] = RF_INF;
range->pmin[PROJ_FACTOR] = - 1e5;
range->pmax[PROJ_FACTOR] = 1e5;
range->openmin[PROJ_FACTOR] = true;
range->openmax[PROJ_FACTOR] = true;
}
void declarefct(double VARIABLE_IS_NOT_USED *x, model *cov, double *v){
int vdimSq = VDIM0 * VDIM1;
for (int i=0; i<vdimSq; v[i++] = 0.0);
}
void rangedeclare(model *cov, range_type *range){
int k = DefList[COVNR].kappas;
for (int i=0; i<k; i++) {
range->min[i] = RF_NEGINF;
range->max[i] = RF_INF;
range->pmin[i] = RF_NEGINF;
range->pmax[i] = RF_INF;
range->openmin[i] = true;
range->openmax[i] = true;
}
}
int checkdeclare(model *cov) {
// defn *C = DefList + COVNR;
int v = cov->calling->vdim[0];
// kappas = C->kappas;
// for (int i=0; i<kappas; i++) if (!PisNULL(i)) P(i)[0]=RF_NA;
if (v <= 0) v=1;
VDIM0 = VDIM1 = v;
RETURN_NOERROR;
}