Commit 8a6f1af9 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

cleanup

parent 4bc2e132
......@@ -56,12 +56,12 @@ void detrend(long nts, double missval1, double *array1, double *array2)
n++;
}
work1 = DIV(SUB(sumjx, DIV(MUL(sumx, sumj), n) ),
SUB(sumjj, DIV(MUL(sumj, sumj), n)) );
work2 = SUB(DIV(sumx, n), MUL(work1, DIV(sumj, n)));
work1 = DIVMN( SUBMN(sumjx, DIVMN( MULMN(sumx, sumj), n) ),
SUBMN(sumjj, DIVMN( MULMN(sumj, sumj), n)) );
work2 = SUBMN( DIVMN(sumx, n), MULMN(work1, DIVMN(sumj, n)));
for ( j = 0; j < nts; j++ )
array2[j] = SUB(array1[j], ADD(work2, MUL(j, work1)));
array2[j] = SUBMN(array1[j], ADDMN(work2, MULMN(j, work1)));
}
......
......@@ -56,8 +56,8 @@ double correlation_s(const double * restrict in0, const double * restrict in1,
}
out = IS_NOT_EQUAL(wsum0, 0) ?
DIV((sum01 * wsum0 - sum0 * sum1),
SQRT((sum00 * wsum0 - sum0 * sum0) *
DIVMN((sum01 * wsum0 - sum0 * sum1),
SQRTMN((sum00 * wsum0 - sum0 * sum0) *
(sum11 * wsum0 - sum1 * sum1))) : missval1;
return (out);
......
......@@ -680,27 +680,24 @@ void ctl_options(FILE *gdp, int yrev, int zrev, int sequential, int bigendian, i
static
void ctl_undef(FILE *gdp, int vlistID)
{
double missval;
missval = vlistInqVarMissval(vlistID, 0);
double missval = vlistInqVarMissval(vlistID, 0);
fprintf(gdp, "UNDEF %g\n", missval);
}
static
void ctl_vars(FILE *gdp, int filetype, int vlistID, int nvarsout, int *vars)
{
int varID, nvars;
int ltype, code;
int zaxisID, nlev;
int i, j;
int len;
char varname[CDI_MAX_NAME], varlongname[CDI_MAX_NAME], varunits[CDI_MAX_NAME];
nvars = vlistNvars(vlistID);
int nvars = vlistNvars(vlistID);
fprintf(gdp, "VARS %d\n", nvarsout);
for ( varID = 0; varID < nvars; varID++ )
for ( int varID = 0; varID < nvars; varID++ )
{
if ( vars[varID] == TRUE )
{
......
......@@ -145,7 +145,7 @@ void *Math(void *argument)
break;
case SQRT:
for ( i = 0; i < gridsize; i++ )
array2[i] = DBL_IS_EQUAL(array1[i], missval1) ? missval1 : SQRT(array1[i]);
array2[i] = DBL_IS_EQUAL(array1[i], missval1) ? missval1 : SQRTMN(array1[i]);
break;
case EXP:
for ( i = 0; i < gridsize; i++ )
......
......@@ -139,14 +139,14 @@ void *Regres(void *argument)
for ( i = 0; i < gridsize; i++ )
{
temp1 = SUB(work[2][varID][levelID].ptr[i],
DIV(MUL(work[0][varID][levelID].ptr[i], work[3][varID][levelID].ptr[i]), work[4][varID][levelID].ptr[i]));
temp2 = SUB(work[1][varID][levelID].ptr[i],
DIV(MUL(work[0][varID][levelID].ptr[i], work[0][varID][levelID].ptr[i]), work[4][varID][levelID].ptr[i]));
field2.ptr[i] = DIV(temp1, temp2);
field1.ptr[i] = SUB(DIV(work[3][varID][levelID].ptr[i], work[4][varID][levelID].ptr[i]),
MUL(DIV(work[0][varID][levelID].ptr[i], work[4][varID][levelID].ptr[i]), field2.ptr[i]));
temp1 = SUBMN(work[2][varID][levelID].ptr[i],
DIVMN( MULMN(work[0][varID][levelID].ptr[i], work[3][varID][levelID].ptr[i]), work[4][varID][levelID].ptr[i]));
temp2 = SUBMN(work[1][varID][levelID].ptr[i],
DIVMN( MULMN(work[0][varID][levelID].ptr[i], work[0][varID][levelID].ptr[i]), work[4][varID][levelID].ptr[i]));
field2.ptr[i] = DIVMN(temp1, temp2);
field1.ptr[i] = SUBMN( DIVMN(work[3][varID][levelID].ptr[i], work[4][varID][levelID].ptr[i]),
MULMN( DIVMN(work[0][varID][levelID].ptr[i], work[4][varID][levelID].ptr[i]), field2.ptr[i]));
}
/*
......
......@@ -82,7 +82,7 @@ void *Runstat(void *argument)
int lmean = operfunc == func_mean || operfunc == func_avg;
int lstd = operfunc == func_std || operfunc == func_std1;
int lvarstd = operfunc == func_std || operfunc == func_var || operfunc == func_std1 || operfunc == func_var1;
double divisor = operfunc == func_std1 || operfunc == func_var1;
int divisor = operfunc == func_std1 || operfunc == func_var1;
int streamID1 = streamOpenRead(cdoStreamName(0));
......
......@@ -46,7 +46,7 @@ void *Seasstat(void *argument)
int vdate1 = 0, vtime1 = 0;
int nrecs;
int varID, levelID, recID;
long nsets;
int nsets;
int i;
int year, month, day, seas, seas0 = 0;
int nmiss;
......@@ -76,7 +76,7 @@ void *Seasstat(void *argument)
int lmean = operfunc == func_mean || operfunc == func_avg;
int lstd = operfunc == func_std || operfunc == func_std1;
int lvarstd = operfunc == func_std || operfunc == func_var || operfunc == func_std1 || operfunc == func_var1;
double divisor = operfunc == func_std1 || operfunc == func_var1;
int divisor = operfunc == func_std1 || operfunc == func_var1;
int streamID1 = streamOpenRead(cdoStreamName(0));
......
......@@ -113,7 +113,7 @@ void *Subtrend(void *argument)
missval1 = missval;
missval2 = missval;
for ( i = 0; i < gridsize; i++ )
field4.ptr[i] = SUB(field1.ptr[i], ADD(vars2[varID][levelID].ptr[i], MUL(vars3[varID][levelID].ptr[i], tsID)));
field4.ptr[i] = SUBMN(field1.ptr[i], ADDMN(vars2[varID][levelID].ptr[i], MULMN(vars3[varID][levelID].ptr[i], tsID)));
nmiss = 0;
for ( i = 0; i < gridsize; i++ )
......
......@@ -75,7 +75,7 @@ void *Timselstat(void *argument)
int lmean = operfunc == func_mean || operfunc == func_avg;
int lstd = operfunc == func_std || operfunc == func_std1;
int lvarstd = operfunc == func_std || operfunc == func_var || operfunc == func_std1 || operfunc == func_var1;
double divisor = operfunc == func_std1 || operfunc == func_var1;
int divisor = operfunc == func_std1 || operfunc == func_var1;
int streamID1 = streamOpenRead(cdoStreamName(0));
......
......@@ -72,17 +72,22 @@
#include "pstream.h"
typedef struct {
short varID;
short levelID;
} recinfo_t;
void *Timstat(void *argument)
{
enum {HOUR_LEN=4, DAY_LEN=6, MON_LEN=8, YEAR_LEN=10};
int timestat_date = TIMESTAT_MEAN;
int vdate0 = 0, vtime0 = 0;
int nrecs;
int varID, levelID, recID;
int varID, levelID;
int streamID3 = -1;
int vlistID3, taxisID3 = -1;
int nmiss;
int nlevel;
int lvfrac = FALSE;
int nwpv; // number of words per value; real:1 complex:2
char indate1[DATE_LEN+1], indate2[DATE_LEN+1];
......@@ -144,7 +149,7 @@ void *Timstat(void *argument)
int lmean = operfunc == func_mean || operfunc == func_avg;
int lstd = operfunc == func_std || operfunc == func_std1;
int lvarstd = operfunc == func_std || operfunc == func_var || operfunc == func_std1 || operfunc == func_var1;
double divisor = operfunc == func_std1 || operfunc == func_var1;
int divisor = operfunc == func_std1 || operfunc == func_var1;
if ( operfunc == func_mean )
{
......@@ -176,8 +181,8 @@ void *Timstat(void *argument)
if ( taxisInqType(taxisID2) == TAXIS_FORECAST ) taxisDefType(taxisID2, TAXIS_RELATIVE);
vlistDefTaxis(vlistID2, taxisID2);
int nvars = vlistNvars(vlistID1);
int nrecords = vlistNrecs(vlistID1);
int nvars = vlistNvars(vlistID1);
int maxrecs = vlistNrecs(vlistID1);
if ( cmplen == 0 && CDO_Reduce_Dim )
for ( varID = 0; varID < nvars; ++varID )
......@@ -221,8 +226,7 @@ void *Timstat(void *argument)
streamDefVlist(streamID3, vlistID3);
}
int *recVarID = (int*) Malloc(nrecords*sizeof(int));
int *recLevelID = (int*) Malloc(nrecords*sizeof(int));
recinfo_t *recinfo = (recinfo_t *) malloc(maxrecs*sizeof(recinfo_t));
dtlist_type *dtlist = dtlist_new();
dtlist_set_stat(dtlist, timestat_date);
......@@ -231,9 +235,16 @@ void *Timstat(void *argument)
int gridsize = vlistGridsizeMax(vlistID1);
if ( vlistNumber(vlistID1) != CDI_REAL ) gridsize *= 2;
int FIELD_MEMTYPE = 0;
if ( CDO_Memtype == MEMTYPE_FLOAT ) FIELD_MEMTYPE = MEMTYPE_FLOAT;
field_t field;
field_init(&field);
field.ptr = (double*) Malloc(gridsize*sizeof(double));
field.memtype = FIELD_MEMTYPE;
if ( FIELD_MEMTYPE == MEMTYPE_FLOAT )
field.ptrf = (float*) Malloc(gridsize*sizeof(float));
else
field.ptr = (double*) Malloc(gridsize*sizeof(double));
field_t **vars1 = field_malloc(vlistID1, FIELD_PTR);
field_t **samp1 = field_malloc(vlistID1, FIELD_NONE);
......@@ -244,7 +255,7 @@ void *Timstat(void *argument)
int otsID = 0;
while ( TRUE )
{
long nsets = 0;
int nsets = 0;
while ( (nrecs = streamInqTimestep(streamID1, tsID)) )
{
dtlist_taxisInqTimestep(dtlist, taxisID1, nsets);
......@@ -256,30 +267,32 @@ void *Timstat(void *argument)
if ( DATE_IS_NEQ(indate1, indate2, cmplen) ) break;
for ( recID = 0; recID < nrecs; recID++ )
for ( int recID = 0; recID < nrecs; recID++ )
{
streamInqRecord(streamID1, &varID, &levelID);
if ( tsID == 0 )
{
recVarID[recID] = varID;
recLevelID[recID] = levelID;
recinfo[recID].varID = varID;
recinfo[recID].levelID = levelID;
}
nwpv = vars1[varID][levelID].nwpv;
gridsize = gridInqSize(vars1[varID][levelID].grid);
field_t *pvar1 = &vars1[varID][levelID];
nwpv = pvar1->nwpv;
gridsize = gridInqSize(pvar1->grid);
if ( nsets == 0 )
{
streamReadRecord(streamID1, vars1[varID][levelID].ptr, &nmiss);
vars1[varID][levelID].nmiss = nmiss;
streamReadRecord(streamID1, pvar1->ptr, &nmiss);
pvar1->nmiss = nmiss;
if ( nmiss > 0 || samp1[varID][levelID].ptr )
{
if ( samp1[varID][levelID].ptr == NULL )
samp1[varID][levelID].ptr = (double*) Malloc(nwpv*gridsize*sizeof(double));
for ( int i = 0; i < nwpv*gridsize; i++ )
if ( DBL_IS_EQUAL(vars1[varID][levelID].ptr[i], vars1[varID][levelID].missval) )
if ( DBL_IS_EQUAL(pvar1->ptr[i], pvar1->missval) )
samp1[varID][levelID].ptr[i] = 0;
else
samp1[varID][levelID].ptr[i] = 1;
......@@ -287,10 +300,13 @@ void *Timstat(void *argument)
}
else
{
streamReadRecord(streamID1, field.ptr, &field.nmiss);
if ( CDO_Memtype == MEMTYPE_FLOAT )
streamReadRecordF(streamID1, field.ptrf, &field.nmiss);
else
streamReadRecord(streamID1, field.ptr, &field.nmiss);
field.size = gridsize;
field.grid = vars1[varID][levelID].grid;
field.missval = vars1[varID][levelID].missval;
field.grid = pvar1->grid;
field.missval = pvar1->missval;
if ( field.nmiss > 0 || samp1[varID][levelID].ptr )
{
if ( samp1[varID][levelID].ptr == NULL )
......@@ -301,29 +317,34 @@ void *Timstat(void *argument)
}
for ( int i = 0; i < nwpv*gridsize; i++ )
if ( !DBL_IS_EQUAL(field.ptr[i], vars1[varID][levelID].missval) )
if ( !DBL_IS_EQUAL(field.ptr[i], pvar1->missval) )
samp1[varID][levelID].ptr[i]++;
}
if ( lvarstd )
{
farsumq(&vars2[varID][levelID], field);
farsum(&vars1[varID][levelID], field);
field_t *pvar2 = &vars2[varID][levelID];
farsumq(pvar2, field);
farsum(pvar1, field);
}
else
{
farfun(&vars1[varID][levelID], field, operfunc);
farfun(pvar1, field, operfunc);
}
}
}
if ( nsets == 0 && lvarstd )
for ( varID = 0; varID < nvars; varID++ )
{
for ( int recID = 0; recID < maxrecs; recID++ )
{
int varID = recinfo[recID].varID;
int levelID = recinfo[recID].levelID;
field_t *pvar1 = &vars1[varID][levelID];
field_t *pvar2 = &vars2[varID][levelID];
if ( vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;
nlevel = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID));
for ( levelID = 0; levelID < nlevel; levelID++ )
farmoq(&vars2[varID][levelID], vars1[varID][levelID]);
farmoq(pvar2, *pvar1);
}
vdate0 = vdate;
......@@ -335,42 +356,38 @@ void *Timstat(void *argument)
if ( nrecs == 0 && nsets == 0 ) break;
if ( lmean )
for ( varID = 0; varID < nvars; varID++ )
{
for ( int recID = 0; recID < maxrecs; recID++ )
{
int varID = recinfo[recID].varID;
int levelID = recinfo[recID].levelID;
field_t *pvar1 = &vars1[varID][levelID];
if ( vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;
nlevel = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID));
for ( levelID = 0; levelID < nlevel; levelID++ )
{
if ( samp1[varID][levelID].ptr == NULL )
farcdiv(&vars1[varID][levelID], (double)nsets);
else
{
// farround(&samp1[varID][levelID]); not necessary
fardiv(&vars1[varID][levelID], samp1[varID][levelID]);
}
}
if ( samp1[varID][levelID].ptr == NULL )
farcdiv(pvar1, (double)nsets);
else
fardiv(pvar1, samp1[varID][levelID]);
}
else if ( lvarstd )
for ( varID = 0; varID < nvars; varID++ )
{
if ( vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;
nlevel = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID));
for ( levelID = 0; levelID < nlevel; levelID++ )
{
if ( samp1[varID][levelID].ptr == NULL )
{
if ( lstd )
farcstd(&vars1[varID][levelID], vars2[varID][levelID], nsets, divisor);
else
farcvar(&vars1[varID][levelID], vars2[varID][levelID], nsets, divisor);
}
else
{
if ( lstd )
farstd(&vars1[varID][levelID], vars2[varID][levelID], samp1[varID][levelID], divisor);
else
farvar(&vars1[varID][levelID], vars2[varID][levelID], samp1[varID][levelID], divisor);
}
for ( int recID = 0; recID < maxrecs; recID++ )
{
int varID = recinfo[recID].varID;
int levelID = recinfo[recID].levelID;
field_t *pvar1 = &vars1[varID][levelID];
field_t *pvar2 = &vars2[varID][levelID];
if ( vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;
if ( samp1[varID][levelID].ptr == NULL )
{
if ( lstd ) farcstd(pvar1, *pvar2, nsets, divisor);
else farcvar(pvar1, *pvar2, nsets, divisor);
}
else
{
if ( lstd ) farstd(pvar1, *pvar2, samp1[varID][levelID], divisor);
else farvar(pvar1, *pvar2, samp1[varID][levelID], divisor);
}
}
......@@ -383,35 +400,36 @@ void *Timstat(void *argument)
}
if ( lvfrac && operfunc == func_mean )
for ( varID = 0; varID < nvars; varID++ )
{
for ( int recID = 0; recID < maxrecs; recID++ )
{
int varID = recinfo[recID].varID;
int levelID = recinfo[recID].levelID;
field_t *pvar1 = &vars1[varID][levelID];
if ( vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;
nlevel = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID));
for ( levelID = 0; levelID < nlevel; levelID++ )
{
nwpv = vars1[varID][levelID].nwpv;
gridsize = gridInqSize(vars1[varID][levelID].grid);
missval = vars1[varID][levelID].missval;
if ( samp1[varID][levelID].ptr )
{
int irun = 0;
for ( int i = 0; i < nwpv*gridsize; ++i )
{
if ( (samp1[varID][levelID].ptr[i] / nsets) < vfrac )
{
vars1[varID][levelID].ptr[i] = missval;
irun++;
}
}
if ( irun )
{
nmiss = 0;
for ( int i = 0; i < nwpv*gridsize; ++i )
if ( DBL_IS_EQUAL(vars1[varID][levelID].ptr[i], missval) ) nmiss++;
vars1[varID][levelID].nmiss = nmiss;
}
}
nwpv = pvar1->nwpv;
gridsize = gridInqSize(pvar1->grid);
missval = pvar1->missval;
if ( samp1[varID][levelID].ptr )
{
int irun = 0;
for ( int i = 0; i < nwpv*gridsize; ++i )
{
if ( (samp1[varID][levelID].ptr[i] / nsets) < vfrac )
{
pvar1->ptr[i] = missval;
irun++;
}
}
if ( irun )
{
nmiss = 0;
for ( int i = 0; i < nwpv*gridsize; ++i )
if ( DBL_IS_EQUAL(pvar1->ptr[i], missval) ) nmiss++;
pvar1->nmiss = nmiss;
}
}
}
......@@ -424,15 +442,17 @@ void *Timstat(void *argument)
streamDefTimestep(streamID3, otsID);
}
for ( recID = 0; recID < nrecords; recID++ )
for ( int recID = 0; recID < maxrecs; recID++ )
{
varID = recVarID[recID];
levelID = recLevelID[recID];
int varID = recinfo[recID].varID;
int levelID = recinfo[recID].levelID;
field_t *pvar1 = &vars1[varID][levelID];
if ( otsID && vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;
streamDefRecord(streamID2, varID, levelID);
streamWriteRecord(streamID2, vars1[varID][levelID].ptr, vars1[varID][levelID].nmiss);
streamWriteRecord(streamID2, pvar1->ptr, pvar1->nmiss);
if ( cdoDiag )
{
if ( samp1[varID][levelID].ptr )
......@@ -459,9 +479,7 @@ void *Timstat(void *argument)
streamClose(streamID1);
if ( field.ptr ) Free(field.ptr);
if ( recVarID ) Free(recVarID);
if ( recLevelID ) Free(recLevelID);
Free(recinfo);
cdoFinish();
......
......@@ -42,15 +42,15 @@ int correlation_t(long gridsize, double missval1, double missval2, int *nofvals,
if ( nvals > 0 )
{
temp0 = MUL(work0[i], work1[i]);
temp1 = SUB(work4[i], DIV(temp0, nvals));
temp2 = MUL(work0[i], work0[i]);
temp3 = MUL(work1[i], work1[i]);
temp4 = SUB(work2[i], DIV(temp2, nvals));
temp5 = SUB(work3[i], DIV(temp3, nvals));
temp6 = MUL(temp4, temp5);
temp0 = MULMN(work0[i], work1[i]);
temp1 = SUBMN(work4[i], DIVMN(temp0, nvals));
temp2 = MULMN(work0[i], work0[i]);
temp3 = MULMN(work1[i], work1[i]);
temp4 = SUBMN(work2[i], DIVMN(temp2, nvals));
temp5 = SUBMN(work3[i], DIVMN(temp3, nvals));
temp6 = MULMN(temp4, temp5);
cor = DIV(temp1, SQRT(temp6));
cor = DIVMN(temp1, SQRTMN(temp6));
if ( cor < -1 ) cor = -1;
else if ( cor > 1 ) cor = 1;
......@@ -87,8 +87,8 @@ int covariance_t(long gridsize, double missval1, double missval2, int *nofvals,
if ( nvals > 0 )
{
temp = DIV(MUL(work0[i], work1[i]), dnvals*dnvals);
covar = SUB(DIV(work2[i], dnvals), temp);
temp = DIVMN( MULMN(work0[i], work1[i]), dnvals*dnvals);
covar = SUBMN( DIVMN(work2[i], dnvals), temp);
if ( DBL_IS_EQUAL(covar, missval1) ) nmiss++;
}
......
......@@ -235,11 +235,11 @@ void *Timstat3(void *argument)
fnvals0 = iwork[0][varID][levelID][i];
fnvals1 = iwork[1][varID][levelID][i];
temp0 = DIV(MUL(fwork[0][varID][levelID].ptr[i], fwork[0][varID][levelID].ptr[i]), fnvals0);
temp1 = DIV(MUL(fwork[2][varID][levelID].ptr[i], fwork[2][varID][levelID].ptr[i]), fnvals1);
temp2 = SUB(fwork[1][varID][levelID].ptr[i], temp0);
temp3 = SUB(fwork[3][varID][levelID].ptr[i], temp1);
statistic = DIV(temp2, ADD(temp2, MUL(rconst, temp3)));
temp0 = DIVMN( MULMN(fwork[0][varID][levelID].ptr[i], fwork[0][varID][levelID].ptr[i]), fnvals0);
temp1 = DIVMN( MULMN(fwork[2][varID][levelID].ptr[i], fwork[2][varID][levelID].ptr[i]), fnvals1);
temp2 = SUBMN(fwork[1][varID][levelID].ptr[i], temp0);
temp3 = SUBMN(fwork[3][varID][levelID].ptr[i], temp1);
statistic = DIVMN(temp2, ADDMN(temp2, MULMN(rconst, temp3)));
if ( fnvals0 <= 1 || fnvals1 <= 1 )
fractil_1 = fractil_2 = missval1;
......@@ -271,35 +271,35 @@ void *Timstat3(void *argument)
for ( j = 0; j < n_in; j++ )
{
fnvals = iwork[j][varID][levelID][i];
tmp = DIV(MUL(fwork[2*j][varID][levelID].ptr[i], fwork[2*j][varID][levelID].ptr[i]), fnvals);
temp0 = ADD(temp0, DIV(SUB(fwork[2*j+1][varID][levelID].ptr[i], tmp), var_factor[j]));
deg_of_freedom = ADD(deg_of_freedom, fnvals);
tmp = DIVMN( MULMN(fwork[2*j][varID][levelID].ptr[i], fwork[2*j][varID][levelID].ptr[i]), fnvals);
temp0 = ADDMN(temp0, DIVMN( SUBMN(fwork[2*j+1][varID][levelID].ptr[i], tmp), var_factor[j]));
deg_of_freedom = ADDMN(deg_of_freedom, fnvals);
}
if ( !DBL_IS_EQUAL(temp0, missval1) && temp0 < 0 ) /* This is possible because */
temp0 = 0; /* of rounding errors */
stddev_estimator = SQRT(DIV(temp0, deg_of_freedom));
stddev_estimator = SQRTMN( DIVMN(temp0, deg_of_freedom));
mean_estimator = -rconst;
for ( j = 0; j < n_in; j++ )
{
fnvals = iwork[j][varID][levelID][i];
mean_estimator = ADD(mean_estimator,
MUL(mean_factor[j],
DIV(fwork[2*j][varID][levelID].ptr[i], fnvals)));
mean_estimator = ADDMN(mean_estimator,
MULMN(mean_factor[j],
DIVMN(fwork[2*j][varID][levelID].ptr[i], fnvals)));
}
temp1 = 0;
for ( j = 0; j < n_in; j++ )
{
fnvals = iwork[j][varID][levelID][i];
temp1 = ADD(temp1, DIV(MUL(MUL(mean_factor[j], mean_factor[j]),
temp1 = ADDMN(temp1, DIVMN( MULMN( MULMN(mean_factor[j], mean_factor[j]),
var_factor[j]), fnvals));
}
norm = SQRT(temp1);
norm = SQRTMN(temp1);
temp2 = DIV(DIV(mean_estimator, norm), stddev_estimator);
temp2 = DIVMN( DIVMN(mean_estimator, norm), stddev_estimator);
fractil = deg_of_freedom < 1 ? missval1 :
student_t_inv (deg_of_freedom, 1 - risk/2, __func__);
......
......@@ -145,14 +145,14 @@ void *Trend(void *argument)