From 0fbbf6e63bc8c43003d972ac9036e0ce2dc19b8f Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Tue, 4 Apr 2023 19:25:56 +0200 Subject: [PATCH 1/6] Allow -E+n for modelnorm output See forum post for the origin of this PR. --- src/greenspline.c | 27 ++++++++++++++++++++++----- 1 file changed, 22 insertions(+), 5 deletions(-) diff --git a/src/greenspline.c b/src/greenspline.c index c068f984e2e..0e702e3b993 100644 --- a/src/greenspline.c +++ b/src/greenspline.c @@ -82,8 +82,9 @@ struct GREENSPLINE_CTRL { bool active; char *information; } D; - struct GREENSPLINE_E { /* -E[] */ + struct GREENSPLINE_E { /* -E[][+n] */ bool active; + bool norm; unsigned int mode; char *file; } E; @@ -233,7 +234,7 @@ static int usage (struct GMTAPI_CTRL *API, int level) { const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE); if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR); GMT_Usage (API, 0, "usage: %s [] -G [-A+f] [-C[[n|r|v][%%]][+c][+f][+i][+n]] " - "[-D] [-E[]] [-I[/[/]]] [-L[t][r]] [-N] [-Q[|]] " + "[-D] [-E[][+n]] [-I[/[/]]] [-L[t][r]] [-N] [-Q[|]] " "[-R/[//[//]]] [-Sc|l|t|r|p|q[]] [-T] " "[%s] [-W[w]] [-Z] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s]%s[%s] [%s]\n", name, GMT_V_OPT,GMT_bi_OPT, GMT_d_OPT, GMT_e_OPT, GMT_f_OPT, GMT_g_OPT, GMT_h_OPT, GMT_i_OPT, @@ -276,10 +277,11 @@ static int usage (struct GMTAPI_CTRL *API, int level) { GMT_Usage (API, 3, "+n Stop execution after reporting the eigenvalues - no solution is computed."); GMT_Usage (API, -2, "Note: Without -C we use Gauss-Jordan elimination to solve the linear system."); gmt_grdcube_info_syntax (API->GMT, 'D'); - GMT_Usage (API, 1, "\n-E[]"); + GMT_Usage (API, 1, "\n-E[][+n]"); GMT_Usage (API, -2, "Evaluate solution at input locations and report misfit statistics. " "Append a filename to save all data with two extra columns for model and misfit. " "If -C+i|c are used then we instead report the history of model variance and rms misfit."); + GMT_Usage (API, 3, "n: Add two last output columns with model and data norms."); GMT_Usage (API, 1, "\n-I[/[/]]"); GMT_Usage (API, -2, "Specify a regular set of output locations. Give equidistant increment for each dimension. " "Requires -R for specifying the output domain."); @@ -570,8 +572,13 @@ static int parse (struct GMT_CTRL *GMT, struct GREENSPLINE_CTRL *Ctrl, struct GM case 'E': /* Evaluate misfit -E[]*/ n_errors += gmt_M_repeated_module_option (API, Ctrl->E.active); if (opt->arg) { + if ((c = strstr (opt->arg, "+n")) == NULL) { + Ctrl->E.norm = true; + c[0] = '\0'; /* Hide modifier */ + } Ctrl->E.file = strdup (opt->arg); Ctrl->E.mode = 1; + if (Ctrl->E.norm) c[0] = '+'; /* Restore modifier */ } break; case 'G': /* Output file */ @@ -2651,7 +2658,7 @@ EXTERN_MSC int GMT_greenspline (void *V_API, int mode, void *args) { previous = gmt_M_memory_aligned (GMT, NULL, Out->header->size, gmt_grdfloat); gmt_grd_init (GMT, Out->header, options, true); if (Ctrl->E.active) { /* Want to write out misfit as function of eigenvalue */ - uint64_t e_dim[GMT_DIM_SIZE] = {1, 1, n_use, 4+Ctrl->W.active}; + uint64_t e_dim[GMT_DIM_SIZE] = {1, 1, n_use, 4+Ctrl->W.active+Ctrl->E.norm}; eigen = gmt_sort_svd_values (GMT, ssave, nm); /* Get sorted eigenvalues */ if ((E = GMT_Create_Data (API, GMT_IS_DATASET, GMT_IS_NONE, 0, e_dim, NULL, NULL, 0, 0, NULL)) == NULL) { GMT_Report (API, GMT_MSG_ERROR, "Unable to create a data set for saving misfit estimates per eigenvector\n"); @@ -2667,7 +2674,16 @@ EXTERN_MSC int GMT_greenspline (void *V_API, int mode, void *args) { GMT_Report (API, GMT_MSG_INFORMATION, "Evaluate spline for eigenvalue # %d\n", (int)e); gmt_M_memcpy (s, ssave, nm, double); /* Restore original values before call */ (void)gmt_solve_svd (GMT, A, (unsigned int)nm, (unsigned int)nm, v, s, b, 1U, obs, (double)e, GMT_SVD_EIGEN_NUMBER_CUTOFF); + /* obs (hence alpha) now has the solution for the coefficients based on the first e eigenvalues */ + if (Ctrl->E.active) { /* Compute the history of model misfit */ + double modnorm = 0.0; + for (k = 0; k < nm; k++) { + modnorm += obs[k] * obs[k]; + S->data[4][e] = modnorm; /* Model norm at this point */ + } + } + if (Ctrl->Q.active) { /* Derivatives of solution */ #ifdef _OPENMP #pragma omp parallel for private(row,V,col,ij,p,wp,r,C,part) shared(Grid,yp,xp,nm,GMT,Ctrl,X,G,par,Lz,alpha,Out,normalize,norm) @@ -2774,11 +2790,12 @@ EXTERN_MSC int GMT_greenspline (void *V_API, int mode, void *args) { char header[GMT_LEN64] = {""}; sprintf (header, "# eigenno\teigenval\tvar_percent\trms"); if (Ctrl->W.active) strcat (header, "\tchi2"); + if (Ctrl->E.norm) strcat (header, "\tmnorm"); if (GMT_Set_Comment (API, GMT_IS_DATASET, GMT_COMMENT_IS_COLNAMES, header, E)) { Return (API->error); } gmt_set_tableheader (API->GMT, GMT_OUT, true); /* So header is written */ - for (k = 0; k < 5; k++) GMT->current.io.col_type[GMT_OUT][k] = GMT_IS_FLOAT; /* Set plain float column types */ + for (k = 0; k < 6; k++) GMT->current.io.col_type[GMT_OUT][k] = GMT_IS_FLOAT; /* Set plain float column types */ if (GMT_Write_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_NONE, GMT_WRITE_SET, NULL, Ctrl->E.file, E) != GMT_NOERROR) { Return (API->error); } From 93510db203b644fc2691a63e81d948ad68d203ed Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Tue, 4 Apr 2023 19:33:10 +0200 Subject: [PATCH 2/6] Update greenspline.c --- src/greenspline.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/greenspline.c b/src/greenspline.c index 0e702e3b993..c39553c7807 100644 --- a/src/greenspline.c +++ b/src/greenspline.c @@ -281,7 +281,7 @@ static int usage (struct GMTAPI_CTRL *API, int level) { GMT_Usage (API, -2, "Evaluate solution at input locations and report misfit statistics. " "Append a filename to save all data with two extra columns for model and misfit. " "If -C+i|c are used then we instead report the history of model variance and rms misfit."); - GMT_Usage (API, 3, "n: Add two last output columns with model and data norms."); + GMT_Usage (API, 3, "+n: Add two last output columns with model and data norms."); GMT_Usage (API, 1, "\n-I[/[/]]"); GMT_Usage (API, -2, "Specify a regular set of output locations. Give equidistant increment for each dimension. " "Requires -R for specifying the output domain."); From dfbb8ae1a132a778bd516a1e0f95f3c884fd98c3 Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Tue, 4 Apr 2023 21:10:25 +0200 Subject: [PATCH 3/6] Update greenspline.c --- src/greenspline.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/greenspline.c b/src/greenspline.c index c39553c7807..d2709c11f44 100644 --- a/src/greenspline.c +++ b/src/greenspline.c @@ -572,7 +572,7 @@ static int parse (struct GMT_CTRL *GMT, struct GREENSPLINE_CTRL *Ctrl, struct GM case 'E': /* Evaluate misfit -E[]*/ n_errors += gmt_M_repeated_module_option (API, Ctrl->E.active); if (opt->arg) { - if ((c = strstr (opt->arg, "+n")) == NULL) { + if ((c = strstr (opt->arg, "+n"))) { Ctrl->E.norm = true; c[0] = '\0'; /* Hide modifier */ } From c654a71bb9a22de1a877de257a18b58efe6ff4cc Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Tue, 4 Apr 2023 21:14:28 +0200 Subject: [PATCH 4/6] Update greenspline.c --- src/greenspline.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/greenspline.c b/src/greenspline.c index d2709c11f44..05bf1e9552d 100644 --- a/src/greenspline.c +++ b/src/greenspline.c @@ -569,7 +569,7 @@ static int parse (struct GMT_CTRL *GMT, struct GREENSPLINE_CTRL *Ctrl, struct GM Ctrl->D.information = strdup (opt->arg); } break; - case 'E': /* Evaluate misfit -E[]*/ + case 'E': /* Evaluate misfit -E[][+n] */ n_errors += gmt_M_repeated_module_option (API, Ctrl->E.active); if (opt->arg) { if ((c = strstr (opt->arg, "+n"))) { From 8f4d78c09d4f49a28ce2651c45fe01c8053cac58 Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Tue, 4 Apr 2023 21:49:07 +0200 Subject: [PATCH 5/6] Update greenspline.c --- src/greenspline.c | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/src/greenspline.c b/src/greenspline.c index 05bf1e9552d..1f97945b915 100644 --- a/src/greenspline.c +++ b/src/greenspline.c @@ -2646,7 +2646,7 @@ EXTERN_MSC int GMT_greenspline (void *V_API, int mode, void *args) { unsigned int width = urint (floor (log10 ((double)n_use))) + 1; /* Width of maximum integer needed */ uint64_t e, p; gmt_grdfloat *current = NULL, *previous = NULL; - double l2_sum_n = 0.0, l2_sum_e = 0.0, predicted; + double l2_sum_n = 0.0, l2_sum_e = 0.0, predicted, modnorm = 0.0; static char *mkind[3] = {"", "Incremental", "Cumulative"}; char file[PATH_MAX] = {""}; struct GMT_SINGULAR_VALUE *eigen = NULL; @@ -2676,12 +2676,10 @@ EXTERN_MSC int GMT_greenspline (void *V_API, int mode, void *args) { (void)gmt_solve_svd (GMT, A, (unsigned int)nm, (unsigned int)nm, v, s, b, 1U, obs, (double)e, GMT_SVD_EIGEN_NUMBER_CUTOFF); /* obs (hence alpha) now has the solution for the coefficients based on the first e eigenvalues */ - if (Ctrl->E.active) { /* Compute the history of model misfit */ - double modnorm = 0.0; - for (k = 0; k < nm; k++) { + if (Ctrl->E.norm) { /* Compute the history of model misfit */ + for (k = 0; k < nm; k++) modnorm += obs[k] * obs[k]; - S->data[4][e] = modnorm; /* Model norm at this point */ - } + S->data[4][e] = modnorm; /* Model norm at this point */ } if (Ctrl->Q.active) { /* Derivatives of solution */ From 31bcb960c5888bb682c87106a894450e1eac4816 Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Mon, 17 Apr 2023 20:23:33 +0200 Subject: [PATCH 6/6] Fix column and sqrt --- src/greenspline.c | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/greenspline.c b/src/greenspline.c index 1f97945b915..1d0d8edc43b 100644 --- a/src/greenspline.c +++ b/src/greenspline.c @@ -1595,7 +1595,7 @@ GMT_LOCAL void greenspline_set_filename (char *name, unsigned int k, unsigned in EXTERN_MSC int GMT_greenspline (void *V_API, int mode, void *args) { openmp_int col, row; uint64_t n_read, p, k, i, j, seg, m, n, nm, n_cr = 0, n_ok = 0, ij, ji, ii, n_duplicates = 0, n_skip = 0; - unsigned int dimension = 0, normalize, n_cols, n_layers = 1, w_col, L_Max = 0; + unsigned int dimension = 0, normalize, n_cols, n_layers = 1, w_col, e_col, L_Max = 0; int64_t *kolumn = NULL; size_t n_alloc; int error = GMT_NOERROR, out_ID, way, n_columns, n_use; @@ -2659,6 +2659,7 @@ EXTERN_MSC int GMT_greenspline (void *V_API, int mode, void *args) { gmt_grd_init (GMT, Out->header, options, true); if (Ctrl->E.active) { /* Want to write out misfit as function of eigenvalue */ uint64_t e_dim[GMT_DIM_SIZE] = {1, 1, n_use, 4+Ctrl->W.active+Ctrl->E.norm}; + e_col = 4 + Ctrl->W.active; eigen = gmt_sort_svd_values (GMT, ssave, nm); /* Get sorted eigenvalues */ if ((E = GMT_Create_Data (API, GMT_IS_DATASET, GMT_IS_NONE, 0, e_dim, NULL, NULL, 0, 0, NULL)) == NULL) { GMT_Report (API, GMT_MSG_ERROR, "Unable to create a data set for saving misfit estimates per eigenvector\n"); @@ -2679,7 +2680,7 @@ EXTERN_MSC int GMT_greenspline (void *V_API, int mode, void *args) { if (Ctrl->E.norm) { /* Compute the history of model misfit */ for (k = 0; k < nm; k++) modnorm += obs[k] * obs[k]; - S->data[4][e] = modnorm; /* Model norm at this point */ + S->data[e_col][e] = sqrt (modnorm); /* Model norm at this point */ } if (Ctrl->Q.active) { /* Derivatives of solution */