Skip to content

Commit

Permalink
Merge branch 'harmonic-notrend' into develop
Browse files Browse the repository at this point in the history
  • Loading branch information
@davidfrantz
davidfrantz committed Jan 9, 2024
2 parents dcb6173 + 5e37f5d commit 6d57afb
Show file tree
Hide file tree
Showing 5 changed files with 120 additions and 32 deletions.
3 changes: 3 additions & 0 deletions docs/source/history/vdev.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,6 +50,9 @@ Develop version
In this case, FORCE will only be able to filter nodata values, but no other quality flags as defined with ``SCREEN_QAI``.
Feature requested by Marcel Schwieder and Felix Lobert.

- In the TSA submodule, it is now possible to disable fitting a monotonic trend within the harmonic model.
The new parameter ``HARMONIC_TREND = TRUE/FALSE`` was added. Default behaviour (as before) is ``TRUE``.

- **FORCE AUX**

- new auxilliary program `force-init`.
Expand Down
6 changes: 6 additions & 0 deletions src/aux-level/param-aux.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1482,6 +1482,12 @@ void write_par_hl_tsi(FILE *fp, bool verbose){
}
fprintf(fp, "RBF_CUTOFF = 0.95\n");

if (verbose){
fprintf(fp, "# Should a monotonic trend be considered in the harmonic interpolation?\n");
fprintf(fp, "# Type: Logical. Valid values: {TRUE,FALSE}\n");
}
fprintf(fp, "HARMONIC_TREND = TRUE\n");

if (verbose){
fprintf(fp, "# Definition of how many modes per season are used for harmonic interpolation,\n");
fprintf(fp, "# i.e. uni-modal (1), bi-modal (2), or tri-modal (3).\n");
Expand Down
141 changes: 109 additions & 32 deletions src/higher-level/interpolate-hl.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -471,26 +471,61 @@ int s;
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++**/
int interpolate_harmonic(tsa_t *ts, small *mask_, int nc, int nt, int nr, int ni, short nodata, par_tsi_t *tsi){
int t, i, k, nk, p;
int ncoef;
int n_coef, i_coef;
double y_pred, y_err;
gsl_matrix *x = NULL, *cov = NULL;
gsl_vector *y = NULL, *c = NULL, *x_pred = NULL;
enum { _TERM_OFF_, _TERM_TREND_, _TERM_UNI_, _TERM_BI_, _TERM_TRI_, _TERM_LEN_ };
bool model_terms[_TERM_LEN_];


if ((ncoef = tsi->harm_nmodes*2 + 2) < 4){
n_coef = 0;
memset(model_terms, 0, _TERM_LEN_);


// add offset
model_terms[_TERM_OFF_] = true;
n_coef++;

// add trend
if (tsi->harm_trend){
model_terms[_TERM_TREND_] = true;
n_coef++;
}

// add uni-modal frequency
if (tsi->harm_nmodes >= 1){
model_terms[_TERM_UNI_] = true;
n_coef += 2;
}

// add bi-modal frequency
if (tsi->harm_nmodes >= 2){
model_terms[_TERM_BI_] = true;
n_coef += 2;
}

// add tri-modal frequency
if (tsi->harm_nmodes >= 3){
model_terms[_TERM_TRI_] = true;
n_coef += 2;
}

// at least have offset and uni-modal frequency
if (n_coef < 3){
printf("not enough coefficients for harmonic fitting. "
"Something is wrong. Abort.\n");
return FAILURE;
}

gsl_set_error_handler_off();

#pragma omp parallel private(i,t,k,nk,x,cov,y,c,x_pred,y_pred,y_err) shared(mask_,ts,nc,nt,nr,ni,tsi,nodata,ncoef) default(none)
#pragma omp parallel private(i,t,k,nk,i_coef,x,cov,y,c,x_pred,y_pred,y_err) shared(mask_,ts,nc,nt,nr,ni,tsi,nodata,n_coef,model_terms) default(none)
{

c = gsl_vector_alloc (ncoef);
cov = gsl_matrix_alloc (ncoef, ncoef);
x_pred = gsl_vector_alloc(ncoef);
c = gsl_vector_alloc (n_coef);
cov = gsl_matrix_alloc (n_coef, n_coef);
x_pred = gsl_vector_alloc(n_coef);

#pragma omp for
for (p=0; p<nc; p++){
Expand All @@ -513,13 +548,13 @@ gsl_vector *y = NULL, *c = NULL, *x_pred = NULL;
}

//if (nk == 0){
if (nk < ncoef){
if (nk < n_coef){
for (i=0; i<ni; i++) ts->tsi_[i][p] = nodata;
continue;
}

//printf("nk: %d\n", nk);
x = gsl_matrix_alloc(nk, ncoef);
x = gsl_matrix_alloc(nk, n_coef);
y = gsl_vector_alloc(nk);


Expand All @@ -528,18 +563,36 @@ gsl_vector *y = NULL, *c = NULL, *x_pred = NULL;
if (ts->tss_[t][p] == nodata ||
ts->d_tss[t].ce < tsi->harm_fit_range[_MIN_].ce ||
ts->d_tss[t].ce > tsi->harm_fit_range[_MAX_].ce){

continue;

} else {
gsl_matrix_set(x, k, 0, 1.0);
gsl_matrix_set(x, k, 1, ts->d_tss[t].ce);
gsl_matrix_set(x, k, 2, cos(2 * M_PI / 365 * ts->d_tss[t].ce));
gsl_matrix_set(x, k, 3, sin(2 * M_PI / 365 * ts->d_tss[t].ce));
if (ncoef > 4) gsl_matrix_set(x, k, 4, cos(4 * M_PI / 365 * ts->d_tss[t].ce));
if (ncoef > 5) gsl_matrix_set(x, k, 5, sin(4 * M_PI / 365 * ts->d_tss[t].ce));
if (ncoef > 6) gsl_matrix_set(x, k, 6, cos(6 * M_PI / 365 * ts->d_tss[t].ce));
if (ncoef > 7) gsl_matrix_set(x, k, 7, sin(6 * M_PI / 365 * ts->d_tss[t].ce));

i_coef = 0;

if (model_terms[_TERM_OFF_]) gsl_matrix_set(x, k, i_coef++, 1.0);

if (model_terms[_TERM_TREND_]) gsl_matrix_set(x, k, i_coef++, ts->d_tss[t].ce);

if (model_terms[_TERM_UNI_]){
gsl_matrix_set(x, k, i_coef++, cos(2 * M_PI / 365 * ts->d_tss[t].ce));
gsl_matrix_set(x, k, i_coef++, sin(2 * M_PI / 365 * ts->d_tss[t].ce));
}

if (model_terms[_TERM_BI_]){
gsl_matrix_set(x, k, i_coef++, cos(4 * M_PI / 365 * ts->d_tss[t].ce));
gsl_matrix_set(x, k, i_coef++, sin(4 * M_PI / 365 * ts->d_tss[t].ce));
}

if (model_terms[_TERM_TRI_]){
gsl_matrix_set(x, k, i_coef++, cos(6 * M_PI / 365 * ts->d_tss[t].ce));
gsl_matrix_set(x, k, i_coef++, sin(6 * M_PI / 365 * ts->d_tss[t].ce));
}

gsl_vector_set(y, k, ts->tss_[t][p]);

k++;

}


Expand All @@ -551,14 +604,26 @@ gsl_vector *y = NULL, *c = NULL, *x_pred = NULL;
// interpolate for each equidistant timestep
for (i=0; i<ni; i++){

gsl_vector_set(x_pred, 0, 1.0);
gsl_vector_set(x_pred, 1, ts->d_tsi[i].ce);
gsl_vector_set(x_pred, 2, cos(2 * M_PI / 365 * ts->d_tsi[i].ce));
gsl_vector_set(x_pred, 3, sin(2 * M_PI / 365 * ts->d_tsi[i].ce));
if (ncoef > 4) gsl_vector_set(x_pred, 4, cos(4 * M_PI / 365 * ts->d_tsi[i].ce));
if (ncoef > 5) gsl_vector_set(x_pred, 5, sin(4 * M_PI / 365 * ts->d_tsi[i].ce));
if (ncoef > 6) gsl_vector_set(x_pred, 6, cos(6 * M_PI / 365 * ts->d_tsi[i].ce));
if (ncoef > 7) gsl_vector_set(x_pred, 7, sin(6 * M_PI / 365 * ts->d_tsi[i].ce));
i_coef = 0;

if (model_terms[_TERM_OFF_]) gsl_vector_set(x_pred, i_coef++, 1.0);

if (model_terms[_TERM_TREND_]) gsl_vector_set(x_pred, i_coef++, ts->d_tsi[i].ce);

if (model_terms[_TERM_UNI_]){
gsl_vector_set(x_pred, i_coef++, cos(2 * M_PI / 365 * ts->d_tsi[i].ce));
gsl_vector_set(x_pred, i_coef++, sin(2 * M_PI / 365 * ts->d_tsi[i].ce));
}

if (model_terms[_TERM_BI_]){
gsl_vector_set(x_pred, i_coef++, cos(4 * M_PI / 365 * ts->d_tsi[i].ce));
gsl_vector_set(x_pred, i_coef++, sin(4 * M_PI / 365 * ts->d_tsi[i].ce));
}

if (model_terms[_TERM_TRI_]){
gsl_vector_set(x_pred, i_coef++, cos(6 * M_PI / 365 * ts->d_tsi[i].ce));
gsl_vector_set(x_pred, i_coef++, sin(6 * M_PI / 365 * ts->d_tsi[i].ce));
}

gsl_multifit_robust_est(x_pred, c, cov, &y_pred, &y_err);
ts->tsi_[i][p] = (short)y_pred;
Expand All @@ -578,14 +643,26 @@ gsl_vector *y = NULL, *c = NULL, *x_pred = NULL;

} else {

gsl_vector_set(x_pred, 0, 1.0);
gsl_vector_set(x_pred, 1, ts->d_tss[t].ce);
gsl_vector_set(x_pred, 2, cos(2 * M_PI / 365 * ts->d_tss[t].ce));
gsl_vector_set(x_pred, 3, sin(2 * M_PI / 365 * ts->d_tss[t].ce));
if (ncoef > 4) gsl_vector_set(x_pred, 4, cos(4 * M_PI / 365 * ts->d_tss[t].ce));
if (ncoef > 5) gsl_vector_set(x_pred, 5, sin(4 * M_PI / 365 * ts->d_tss[t].ce));
if (ncoef > 6) gsl_vector_set(x_pred, 6, cos(6 * M_PI / 365 * ts->d_tss[t].ce));
if (ncoef > 7) gsl_vector_set(x_pred, 7, sin(6 * M_PI / 365 * ts->d_tss[t].ce));
i_coef = 0;

if (model_terms[_TERM_OFF_]) gsl_vector_set(x_pred, i_coef++, 1.0);

if (model_terms[_TERM_TREND_]) gsl_vector_set(x_pred, i_coef++, ts->d_tss[t].ce);

if (model_terms[_TERM_UNI_]){
gsl_vector_set(x_pred, i_coef++, cos(2 * M_PI / 365 * ts->d_tss[t].ce));
gsl_vector_set(x_pred, i_coef++, sin(2 * M_PI / 365 * ts->d_tss[t].ce));
}

if (model_terms[_TERM_BI_]){
gsl_vector_set(x_pred, i_coef++, cos(4 * M_PI / 365 * ts->d_tss[t].ce));
gsl_vector_set(x_pred, i_coef++, sin(4 * M_PI / 365 * ts->d_tss[t].ce));
}

if (model_terms[_TERM_TRI_]){
gsl_vector_set(x_pred, i_coef++, cos(6 * M_PI / 365 * ts->d_tss[t].ce));
gsl_vector_set(x_pred, i_coef++, sin(6 * M_PI / 365 * ts->d_tss[t].ce));
}

gsl_multifit_robust_est(x_pred, c, cov, &y_pred, &y_err);
ts->nrt_[k][p] = (short)(ts->tss_[t][p] - y_pred);
Expand Down
1 change: 1 addition & 0 deletions src/higher-level/param-hl.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -199,6 +199,7 @@ void register_tsa(params_t *params, par_hl_t *phl){
register_intvec_par(params, "RBF_SIGMA", 1, 365, &phl->tsa.tsi.rbf_sigma, &phl->tsa.tsi.rbf_nk);
register_float_par(params, "RBF_CUTOFF", 0, 1, &phl->tsa.tsi.rbf_cutoff);
register_int_par(params, "HARMONIC_MODES", 1, 3, &phl->tsa.tsi.harm_nmodes);
register_bool_par(params, "HARMONIC_TREND", &phl->tsa.tsi.harm_trend);
register_datevec_par(params, "HARMONIC_FIT_RANGE", "1900-01-01", "2099-12-31", &phl->tsa.tsi.harm_fit_range, &phl->tsa.tsi.harm_fit_nrange);
register_int_par(params, "INT_DAY", 1, INT_MAX, &phl->tsa.tsi.step);
register_enum_par(params, "STANDARDIZE_TSI", _TAGGED_ENUM_STD_, _STD_LENGTH_, &phl->tsa.tsi.standard);
Expand Down
1 change: 1 addition & 0 deletions src/higher-level/param-hl.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -183,6 +183,7 @@ typedef struct {
int rbf_nk; // number of kernels for RBF fit
int *rbf_sigma; // sigmas for RBF fit
float rbf_cutoff; // cutoff for RBF fit
int harm_trend; // trend in harmonic model?
int harm_nmodes; // number of modes for harmonic
date_t *harm_fit_range; // date range for fitting harmonic
int harm_fit_nrange; // number of dates for fitting harmonic
Expand Down

0 comments on commit 6d57afb

Please sign in to comment.