Skip to content

Commit ed736b6

Browse files
AllinCottrellAllinCottrell
committed
make matrix_aggregate somewhat more general
1 parent fce716b commit ed736b6

File tree

2 files changed

+102
-145
lines changed

2 files changed

+102
-145
lines changed

‎lib/src/geneval.c‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -13642,7 +13642,7 @@ static NODE *eval_3args_func (NODE *l, NODE *m, NODE *r,
1364213642
if (mat_args == 1 && data_args == 2) {
1364313643
/* matrix arguments: can't be mixed with series/list */
1364413644
p->err = E_TYPES;
13645-
} else if (mat_args == 2 && ym->cols == 1) {
13645+
} else if (mat_args == 2 /* && ym->cols == 1 */) {
1364613646
/* experimental */
1364713647
matrix_version = 1;
1364813648
} else if (mat_args > 0) {

‎lib/src/genfuncs.c‎

Lines changed: 101 additions & 144 deletions
Original file line numberDiff line numberDiff line change
@@ -7593,210 +7593,167 @@ static gretl_matrix *series_aggregate (const gretl_matrix *X,
75937593
return A;
75947594
}
75957595

7596-
#define AGGR_SORT_ALL 1
7596+
static int row_starts_block (int i, const double *y,
7597+
int ny, int n)
7598+
{
7599+
int j;
75977600

7598-
#if AGGR_SORT_ALL
7601+
if (i == 0) {
7602+
return 1;
7603+
} else {
7604+
for (j=0; j<ny; j++) {
7605+
if (y[i] > y[i-1]) {
7606+
return 1;
7607+
}
7608+
y += n;
7609+
}
7610+
}
7611+
7612+
return 0;
7613+
}
75997614

7600-
/* Design: start by forming (y ~ X), where y is a vector, and sorting
7601-
this matrix by its first column. This is an expensive operation,
7602-
but once it's done the rest of the calculation is very
7603-
straightforward.
7615+
/* Design: start by forming (Y ~ X) and sorting this matrix by the Y
7616+
columns. This is an expensive operation, but once it's done the
7617+
rest of the calculation is very straightforward.
76047618
*/
76057619

76067620
gretl_matrix *matrix_aggregate (const gretl_matrix *X,
7607-
const gretl_matrix *y,
7621+
const gretl_matrix *Y,
76087622
const char *func,
76097623
int *err)
76107624
{
76117625
double (*dbuiltin) (int, int, const double *) = NULL;
76127626
int (*ibuiltin) (int, int, const double *) = NULL;
7613-
gretl_matrix *Tmp;
76147627
gretl_matrix *M;
76157628
gretl_matrix *ret;
7616-
double *z, *zsave;
7629+
double *y;
76177630
double rkj;
7631+
int *cols = NULL;
76187632
int *counts = NULL;
7619-
int nx = X->cols;
7620-
int nby = 0;
7621-
int i, j, k, nk;
7633+
int just_count = 0;
7634+
int nx, ny;
7635+
int n, nby = 0;
7636+
int i, j, k;
76227637

7623-
if (X == NULL || y == NULL || X->rows != y->rows) {
7638+
if (func == NULL || !strcmp(func, "null")) {
7639+
just_count = 1;
7640+
} else if (X == NULL || Y == NULL || X->rows != Y->rows) {
76247641
*err = E_NONCONF;
7642+
return NULL;
76257643
} else if (!get_aggregator(func, &dbuiltin, &ibuiltin)) {
7626-
/* can't handle a user-defined function here */
7627-
return series_aggregate(X, y, func, err);
7644+
/* we can't handle a user-defined function here */
7645+
return series_aggregate(X, Y, func, err);
76287646
}
76297647

7630-
if (*err) {
7631-
return NULL;
7648+
nx = just_count ? 0 : X->cols;
7649+
ny = Y->cols;
7650+
n = Y->rows;
7651+
7652+
if (ny > 1) {
7653+
cols = gretl_consecutive_list_new(0, ny - 1);
76327654
}
76337655

7634-
/* allocate space for y ~ X */
7635-
Tmp = gretl_matrix_alloc(X->rows, 1 + X->cols);
7636-
if (Tmp == NULL) {
7637-
*err = E_ALLOC;
7638-
return NULL;
7656+
if (nx == 0) {
7657+
/* make M = sorted Y */
7658+
if (ny > 1) {
7659+
M = gretl_matrix_sort_by_columns(Y, cols, err);
7660+
} else {
7661+
M = gretl_matrix_sort_by_column(Y, 0, err);
7662+
}
7663+
} else {
7664+
/* make M = (Y ~ X), sorted by Y */
7665+
gretl_matrix *Tmp = gretl_matrix_alloc(n, ny + nx);
7666+
7667+
if (Tmp == NULL) {
7668+
*err = E_ALLOC;
7669+
return NULL;
7670+
}
7671+
memcpy(Tmp->val, Y->val, n * ny * sizeof(double));
7672+
memcpy(Tmp->val + n * ny, X->val, n * nx * sizeof(double));
7673+
if (ny > 1) {
7674+
M = gretl_matrix_sort_by_columns(Tmp, cols, err);
7675+
} else {
7676+
M = gretl_matrix_sort_by_column(Tmp, 0, err);
7677+
}
7678+
gretl_matrix_free(Tmp);
76397679
}
76407680

7641-
/* make M = (y ~ X), sorted by y */
7642-
memcpy(Tmp->val, y->val, X->rows * sizeof(double));
7643-
memcpy(Tmp->val + X->rows, X->val, X->rows * nx * sizeof(double));
7644-
M = gretl_matrix_sort_by_column(Tmp, 0, err);
7645-
gretl_matrix_free(Tmp);
76467681
if (*err) {
76477682
return NULL;
76487683
}
76497684

76507685
/* determine the number of 'by' values */
76517686
nby = 0;
7652-
for (i=0; i<M->rows; i++) {
7653-
if (i == 0 || M->val[i] > M->val[i-1]) {
7687+
for (i=0; i<n; i++) {
7688+
if (row_starts_block(i, M->val, ny, n)) {
76547689
nby++;
76557690
}
76567691
}
76577692

76587693
/* allocate the return matrix */
7659-
ret = gretl_zero_matrix_new(nby, 2 + nx);
7694+
ret = gretl_zero_matrix_new(nby, ny + 1 + nx);
76607695
if (ret == NULL) {
76617696
*err = E_ALLOC;
76627697
gretl_matrix_free(M);
7698+
free(cols);
76637699
return NULL;
76647700
}
76657701

7666-
/* get the per-value counts, filling the first two columns
7667-
of @ret as we go
7668-
*/
7702+
/* find and record the per-value counts */
76697703
counts = calloc(nby, sizeof *counts);
76707704
counts[0] = 1;
7671-
for (i=1, k=0; i<M->rows; i++) {
7672-
if (M->val[i] > M->val[i-1]) {
7673-
gretl_matrix_set(ret, k, 0, M->val[i-1]);
7674-
gretl_matrix_set(ret, k, 1, (double) counts[k]);
7675-
counts[++k] = 1;
7676-
if (i == M->rows - 1) {
7677-
gretl_matrix_set(ret, k, 0, M->val[i]);
7678-
gretl_matrix_set(ret, k, 1, 1.0);
7705+
for (i=1, k=0; i<n; i++) {
7706+
y = M->val;
7707+
if (row_starts_block(i, y, ny, n)) {
7708+
for (j=0; j<ny; j++) {
7709+
gretl_matrix_set(ret, k, j, y[i-1]);
7710+
y += n;
76797711
}
7712+
gretl_matrix_set(ret, k, j, (double) counts[k]);
7713+
counts[++k] = 1;
76807714
} else {
76817715
counts[k] += 1;
7682-
if (i == M->rows - 1) {
7683-
gretl_matrix_set(ret, k, 0, M->val[i]);
7684-
gretl_matrix_set(ret, k, 1, (double) counts[k]);
7716+
if (i == n - 1) {
7717+
for (j=0; j<ny; j++) {
7718+
gretl_matrix_set(ret, k, j, y[i]);
7719+
y += n;
7720+
}
7721+
gretl_matrix_set(ret, k, j, (double) counts[k]);
76857722
}
76867723
}
76877724
}
76887725

7689-
z = M->val + M->rows; /* the @X portion of @M */
7726+
if (nx > 0) {
7727+
/* do the @X aggregation */
7728+
double *z = M->val + ny * n; /* the @X portion of @M */
7729+
double *zsave;
7730+
int nk;
76907731

7691-
/* do the actual aggregation */
7692-
for (k=0; k<nby; k++) {
7693-
zsave = z;
7694-
nk = counts[k] - 1;
7695-
for (j=0; j<nx; j++) {
7696-
if (dbuiltin) {
7697-
rkj = (*dbuiltin)(0, nk, z);
7698-
} else {
7699-
rkj = (double) (*ibuiltin)(0, nk, z);
7732+
for (k=0; k<nby; k++) {
7733+
zsave = z;
7734+
nk = counts[k] - 1;
7735+
for (j=0; j<nx; j++) {
7736+
if (dbuiltin) {
7737+
rkj = (*dbuiltin)(0, nk, z);
7738+
} else {
7739+
rkj = (double) (*ibuiltin)(0, nk, z);
7740+
}
7741+
gretl_matrix_set(ret, k, j + ny + 1, rkj);
7742+
/* skip to the next column */
7743+
z += n;
77007744
}
7701-
gretl_matrix_set(ret, k, j+2, rkj);
7702-
/* skip to the next column */
7703-
z += M->rows;
7745+
/* move to the next block of rows */
7746+
z = zsave + counts[k];
77047747
}
7705-
/* move to the next block of rows */
7706-
z = zsave + counts[k];
77077748
}
77087749

77097750
gretl_matrix_free(M);
77107751
free(counts);
7752+
free(cols);
77117753

77127754
return ret;
77137755
}
77147756

7715-
#else /* not AGGR_SORT_ALL */
7716-
7717-
gretl_matrix *matrix_aggregate (const gretl_matrix *X,
7718-
const gretl_matrix *y,
7719-
const char *func,
7720-
int *err)
7721-
{
7722-
double (*dbuiltin) (int, int, const double *) = NULL;
7723-
int (*ibuiltin) (int, int, const double *) = NULL;
7724-
gretl_matrix *yvals = NULL;
7725-
gretl_matrix *sel = NULL;
7726-
gretl_matrix *xi = NULL;
7727-
gretl_matrix *ret = NULL;
7728-
double by, rkj;
7729-
double *z;
7730-
int n, ni, nx, nby;
7731-
int i, j, k;
7732-
7733-
if (X == NULL || y == NULL || X->rows != y->rows) {
7734-
*err = E_NONCONF;
7735-
} else if (!get_aggregator(func, &dbuiltin, &ibuiltin)) {
7736-
/* can't handle a user-defined function here */
7737-
return series_aggregate(X, y, func, err);
7738-
}
7739-
if (!*err) {
7740-
yvals = gretl_matrix_values(y->val, y->rows, OPT_S, err);
7741-
}
7742-
7743-
if (*err) {
7744-
return NULL;
7745-
}
7746-
7747-
n = y->rows;
7748-
nby = yvals->rows;
7749-
nx = X->cols;
7750-
7751-
ret = gretl_matrix_alloc(nby, 2 + nx);
7752-
sel = gretl_matrix_alloc(n, 1);
7753-
7754-
if (ret == NULL || sel == NULL) {
7755-
*err = E_ALLOC;
7756-
gretl_matrix_free(yvals);
7757-
gretl_matrix_free(ret);
7758-
gretl_matrix_free(sel);
7759-
return NULL;
7760-
}
7761-
7762-
for (i=0; i<nby && !*err; i++) {
7763-
by = yvals->val[i];
7764-
ni = 0;
7765-
for (k=0; k<n; k++) {
7766-
if (y->val[k] == by) {
7767-
sel->val[k] = 1.0;
7768-
ni++;
7769-
} else {
7770-
sel->val[k] = 0.0;
7771-
}
7772-
}
7773-
xi = gretl_matrix_bool_sel(X, sel, 1, err);
7774-
if (*err) {
7775-
break;
7776-
}
7777-
gretl_matrix_set(ret, i, 0, by);
7778-
gretl_matrix_set(ret, i, 1, (double) ni);
7779-
z = xi->val;
7780-
for (j=0; j<nx; j++) {
7781-
if (dbuiltin) {
7782-
rkj = (*dbuiltin)(0, ni-1, z);
7783-
} else {
7784-
rkj = (double) (*ibuiltin)(0, ni-1, z);
7785-
}
7786-
gretl_matrix_set(ret, i, j+2, rkj);
7787-
z += ni;
7788-
}
7789-
gretl_matrix_free(xi);
7790-
}
7791-
7792-
gretl_matrix_free(yvals);
7793-
gretl_matrix_free(sel);
7794-
7795-
return ret;
7796-
}
7797-
7798-
#endif /* AGGR_SORT_ALL or not */
7799-
78007757
static int monthly_or_quarterly_dates (const DATASET *dset,
78017758
int pd, double sd0, int T,
78027759
double *x)

0 commit comments

Comments
 (0)