#include "expr.h"
#include "exprP.h"
#include "error.h"
+#include <gsl/gsl_randist.h>
#include <math.h>
#include <errno.h>
#include <stdio.h>
+#include "case.h"
#include "data-in.h"
+#include "dictionary.h"
#include "error.h"
#include "julcal/julcal.h"
#include "magic.h"
#include "misc.h"
+#include "moments.h"
#include "pool.h"
-#include "random.h"
-#include "stats.h"
+#include "settings.h"
#include "str.h"
#include "var.h"
#include "vfm.h"
#include "vfmP.h"
double
-expr_evaluate (const struct expression *e, const struct ccase *c, int case_num,
+expr_evaluate (const struct expression *e, const struct ccase *c, int case_idx,
union value *v)
{
unsigned char *op = e->op;
case OP_CFVAR:
{
int n_args = *op++;
- int nv = 0;
- double sum[2] =
- {0.0, 0.0};
+ double weight, mean, variance;
sp -= n_args - 1;
- for (i = 0; i < n_args; i++)
- if (sp[i].f != SYSMIS)
- {
- nv++;
- sum[0] += sp[i].f;
- sum[1] += sp[i].f * sp[i].f;
- }
- if (nv < *op++)
+
+ moments_of_values (sp, n_args,
+ &weight, &mean, &variance, NULL, NULL);
+
+ if (weight < *op++ || mean == SYSMIS
+ || mean == 0 || variance == SYSMIS)
sp->f = SYSMIS;
else
- sp->f = calc_cfvar (sum, nv);
+ sp->f = sqrt (variance) / mean;
}
break;
case OP_MAX:
case OP_MEAN:
{
int n_args = *op++;
- int nv = 0;
- double sum[1] =
- {0.0};
+ double weight, mean;
sp -= n_args - 1;
- for (i = 0; i < n_args; i++)
- if (sp[i].f != SYSMIS)
- {
- nv++;
- sum[0] += sp[i].f;
- }
- if (nv < *op++)
- sp->f = SYSMIS;
- else
- sp->f = calc_mean (sum, nv);
+
+ moments_of_values (sp, n_args,
+ &weight, &mean, NULL, NULL, NULL);
+ sp->f = weight < *op++ ? SYSMIS : mean;
}
break;
case OP_MIN:
case OP_SD:
{
int n_args = *op++;
- int nv = 0;
- double sum[2];
-
- sum[0] = sum[1] = 0.0;
+ double weight, variance;
sp -= n_args - 1;
- for (i = 0; i < n_args; i++)
- if (sp[i].f != SYSMIS)
- {
- nv++;
- sum[0] += sp[i].f;
- sum[1] += sp[i].f * sp[i].f;
- }
- if (nv < *op++)
- sp->f = SYSMIS;
- else
- sp->f = calc_stddev (calc_variance (sum, nv));
+ moments_of_values (sp, n_args,
+ &weight, NULL, &variance, NULL, NULL);
+ sp->f = weight < *op++ ? SYSMIS : sqrt (variance);
}
break;
case OP_SUM:
case OP_VARIANCE:
{
int n_args = *op++;
- int nv = 0;
- double sum[2];
-
- sum[0] = sum[1] = 0.0;
+ double weight, variance;
sp -= n_args - 1;
- for (i = 0; i < n_args; i++)
- if (sp[i].f != SYSMIS)
- {
- nv++;
- sum[0] += sp[i].f;
- sum[1] += sp[i].f * sp[i].f;
- }
- if (nv < *op++)
- sp->f = SYSMIS;
- else
- sp->f = calc_variance (sum, nv);
+ moments_of_values (sp, n_args,
+ &weight, NULL, &variance, NULL, NULL);
+ sp->f = weight < *op++ ? SYSMIS : variance;
}
break;
break;
case OP_XDATE_JDAY:
if (sp->f != SYSMIS)
- sp->f = 86400. * julian_to_jday (sp->f / 86400.);
+ sp->f = julian_to_jday (sp->f / 86400.);
break;
case OP_XDATE_MDAY:
if (sp->f != SYSMIS)
break;
case OP_NORMAL:
if (sp->f != SYSMIS)
- sp->f *= rng_get_double_normal (pspp_rng ());
+ sp->f = gsl_ran_gaussian (get_rng (), sp->f);
break;
case OP_UNIFORM:
if (sp->f != SYSMIS)
- sp->f *= rng_get_double (pspp_rng ());
+ sp->f *= gsl_rng_uniform (get_rng ());
break;
case OP_SYSMIS:
sp->f = sp->f == SYSMIS || !finite (sp->f);
break;
}
assert (c != NULL);
- sp->f = c->data[v->var[rindx - 1]->fv].f;
+ sp->f = case_num (c, v->var[rindx - 1]->fv);
}
break;
case OP_VEC_ELEM_STR:
sp->c = pool_alloc (e->pool, v->width + 1);
sp->c[0] = v->width;
assert (c != NULL);
- memcpy (&sp->c[1], c->data[v->fv].s, v->width);
+ memcpy (&sp->c[1], case_str (c, v->fv), v->width);
}
break;
case OP_NUM_VAR:
sp++;
assert (c != NULL);
- sp->f = c->data[(*vars)->fv].f;
+ sp->f = case_num (c, (*vars)->fv);
if (is_num_user_missing (sp->f, *vars))
sp->f = SYSMIS;
vars++;
sp->c = pool_alloc (e->pool, width + 1);
sp->c[0] = width;
assert (c != NULL);
- memcpy (&sp->c[1], &c->data[(*vars)->fv], width);
+ memcpy (&sp->c[1], case_str (c, (*vars)->fv), width);
vars++;
}
break;
sp->f = SYSMIS;
else
{
- sp->f = c->data[(*vars)->fv].f;
+ sp->f = case_num (c, (*vars)->fv);
if (is_num_user_missing (sp->f, *vars))
sp->f = SYSMIS;
}
if (c == NULL)
memset (sp->c, ' ', width);
else
- memcpy (&sp->c[1], &c->data[(*vars)->fv], width);
+ memcpy (&sp->c[1], case_str (c, (*vars)->fv), width);
vars++;
}
case OP_NUM_SYS:
sp++;
assert (c != NULL);
- sp->f = c->data[*op++].f == SYSMIS;
+ sp->f = case_num (c, *op++) == SYSMIS;
break;
case OP_NUM_VAL:
sp++;
assert (c != NULL);
- sp->f = c->data[*op++].f;
+ sp->f = case_num (c, *op++);
break;
case OP_CASENUM:
sp++;
- sp->f = case_num;
+ sp->f = case_idx;
break;
case OP_SENTINEL: