diff --git a/doc/rst/source/supplements/x2sys/x2sys_solve.rst b/doc/rst/source/supplements/x2sys/x2sys_solve.rst index 698ada7c803..a52aa1dcf7a 100644 --- a/doc/rst/source/supplements/x2sys/x2sys_solve.rst +++ b/doc/rst/source/supplements/x2sys/x2sys_solve.rst @@ -12,7 +12,7 @@ Synopsis .. include:: ../../common_SYN_OPTs.rst_ -**gmt x2sys_solve** |-C|\ *column* |-T|\ *TAG* |-E|\ **c**\|\ **d**\|\ **g**\|\ **h**\|\ **s**\|\ **y**\|\ **z** +**gmt x2sys_solve** |-C|\ *column* |-T|\ *TAG* |-E|\ **c**\|\ **d**\|\ **g**\|\ **h**\|\ **s**\|\ **t**\|\ **z**\ [**+r**\ [*K*]] [ *COE_list.txt* ] [ |SYN_OPT-V| ] [ |-W|\ [**+u**] ] @@ -53,7 +53,7 @@ Required Arguments .. _-E: -**-E**\ **c**\|\ **d**\|\ **g**\|\ **h**\|\ **s**\|\ **y**\|\ **z** +**-E**\ **c**\|\ **d**\|\ **g**\|\ **h**\|\ **s**\|\ **t**\|\ **z**\ [**+r**\ [*K*]] The correction type you wish to model. Choose among the following functions f(**p**) , where **p** are the *m* parameters per track that we will fit simultaneously using a least @@ -83,6 +83,14 @@ Required Arguments **z** will fit f(**p**) = *a* + *b* \* z (an offset plus a unit scale correction); *z* is the data value at the crossover; records must contain z1, z2, ID1, ID2. + Append **+r**\ [*K*] to **d** or **t** to ridge-regularize the solve + (see the Warning further below): this pulls a poorly-constrained + track's offset back towards a plain, well-behaved correction instead + of letting it blow up, while leaving well-constrained tracks nearly + unaffected. *K* is a trust multiplier: the prior belief is that + |offset| should not need to exceed roughly *K* times the + pre-correction crossover-error standard deviation [10]. + Optional Arguments ------------------ @@ -126,6 +134,26 @@ constraint equations. If you need a particular reference track to have a particular offset (e.g., 0) then you can subtract the offset you found from every track correction and add in the desired offset. +**Warning**: Models that solve for more than one parameter per track +(**d**, **t**, **h**, **g**, **z**) need each track to have not just +at least as many crossings as parameters, but *enough of them, well +distributed along the track*, to separate the parameters reliably. A +track whose crossings are few and/or clustered together leaves its +part of the normal-equation system poorly conditioned. Because that +system is not exactly singular no error is raised, and the +least-squares solution can come back with an offset many times larger +than the crossover errors it was supposed to remove -- silently +corrupting that track instead of correcting it, even though the fit +still looks fine (crossover residuals stay small). **x2sys_solve** +flags such tracks with a warning naming the track, the ratio of its +solved offset to the pre-correction crossover-error scale, and how +many crossings back it. If you see this warning, inspect the flagged +track's connectivity (:doc:`x2sys_report` reports crossings per +track); try appending **+r** to **-Ed** or **-Et** to regularize the +solve (see |-E| above), or fall back to **-Ec** if your survey does +not have enough well-spread crossings per track to support a reliable +drift (or other multi-parameter) estimate. + Input Format ------------ diff --git a/src/x2sys/x2sys_solve.c b/src/x2sys/x2sys_solve.c index 20d8f4ed72c..65458bd1a12 100644 --- a/src/x2sys/x2sys_solve.c +++ b/src/x2sys/x2sys_solve.c @@ -67,6 +67,10 @@ #define F_IS_DRIFT_T 6 /* Subtract a trend with time from each track */ #define F_IS_SCALE_OFF 7 /* Apply a scale and offset to the observations for each track */ +/* Few or clustered crossings leave a track's offset poorly constrained: the fit still closes but the + * offset can come back orders of magnitude too large, and the matrix is not singular so nothing errors */ +#define X2SYS_SOLVE_UNSTABLE_RATIO 20.0 /* Offsets exceeding this many old_stdev's are suspect */ + struct X2SYS_SOLVE_CTRL { struct X2SYS_SOLVE_In { bool active; @@ -76,9 +80,11 @@ struct X2SYS_SOLVE_CTRL { bool active; char *col; } C; - struct X2SYS_SOLVE_E { /* -E */ + struct X2SYS_SOLVE_E { /* -E[+r[]] */ bool active; int mode; + bool regularize; /* Ridge-regularize the d|t solve */ + double K; /* Trust multiplier: |offset| not expected to exceed K * COE st.dev. */ } E; struct X2SYS_SOLVE_T { /* -T */ bool active; @@ -162,6 +168,8 @@ static void *New_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new /* Initialize values whose defaults are not 0/false/NULL */ + C->E.K = 10.0; /* Default +r trust multiplier */ + return (C); } @@ -176,7 +184,7 @@ static void Free_Ctrl (struct GMT_CTRL *GMT, struct X2SYS_SOLVE_CTRL *C) { /* De 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 [] -C -Ec|d|g|h|s|t|z -T [%s] [-W[+u]] [%s] [%s]%s[%s]\n", + GMT_Usage (API, 0, "usage: %s [] -C -Ec|d|g|h|s|t|z[+r[]] -T [%s] [-W[+u]] [%s] [%s]%s[%s]\n", name, GMT_V_OPT, GMT_bi_OPT, GMT_di_OPT, GMT_x_OPT, GMT_PAR_OPT); if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS); @@ -185,7 +193,7 @@ static int usage (struct GMTAPI_CTRL *API, int level) { GMT_Usage (API, 1, "\n File with crossover error data base [standard input]."); GMT_Usage (API, 1, "\n-C"); GMT_Usage (API, -2, "Specify the column name to process (e.g., faa, mag)."); - GMT_Usage (API, 1, "\n-Ec|d|g|h|s|t|z"); + GMT_Usage (API, 1, "\n-Ec|d|g|h|s|t|z[+r[]]"); GMT_Usage (API, -2, "Equation to fit: specify to indicate model to fit per track:"); GMT_Usage (API, 3, "c: Constant offset... Determine offset [Default]."); GMT_Usage (API, 3, "d: Drift by distance. Determine offset and drift-vs-distance rate."); @@ -194,6 +202,13 @@ static int usage (struct GMTAPI_CTRL *API, int level) { GMT_Usage (API, 3, "s: Data scale........ Determine scaling factor."); GMT_Usage (API, 3, "t: Drift over time... Determine offset and drift-vs-time rate."); GMT_Usage (API, 3, "z: Data scale/offset. Determine offset and scaling factor."); + GMT_Usage (API, -2, "Note: d|t|g|h|z solve for more than one parameter per track and need each " + "track to have enough well-distributed crossings to separate them reliably; too few or " + "clustered crossings can give an unstable, oversized offset (a warning is issued when this happens). " + "For d|t, append +r to ridge-regularize the solve against this instability: this pulls poorly-" + "constrained tracks back towards a plain offset (small effect on well-constrained tracks) by " + "assuming |offset| is not expected to exceed about times the pre-correction crossover-error " + "scale [10]."); GMT_Usage (API, 1, "\n-T"); GMT_Usage (API, -2, "Set the system tag for this compilation."); GMT_Message (API, GMT_TIME_NONE, "\n OPTIONAL ARGUMENTS:\n"); @@ -264,6 +279,14 @@ static int parse (struct GMT_CTRL *GMT, struct X2SYS_SOLVE_CTRL *Ctrl, struct GM n_errors++; break; } + if (gmt_validate_modifiers (GMT, opt->arg, 'E', "r", GMT_MSG_ERROR)) n_errors++; + { /* Optional +r[] modifier */ + char txt[GMT_LEN64] = {""}; + if (gmt_get_modifier (opt->arg, 'r', txt)) { + Ctrl->E.regularize = true; + if (txt[0]) Ctrl->E.K = atof (txt); + } + } break; case 'T': n_errors += gmt_M_repeated_module_option (API, Ctrl->T.active); @@ -283,6 +306,11 @@ static int parse (struct GMT_CTRL *GMT, struct X2SYS_SOLVE_CTRL *Ctrl, struct GM n_errors += gmt_M_check_condition (GMT, !Ctrl->T.active || !Ctrl->T.TAG, "Option -T must be used to set the TAG\n"); n_errors += gmt_M_check_condition (GMT, Ctrl->E.mode < 0, "Option -E: Choose among c, d, g, h, s, t and z\n"); + n_errors += gmt_M_check_condition (GMT, Ctrl->E.regularize && Ctrl->E.K <= 0.0, "Option -E: The +r trust multiplier must be positive\n"); + if (Ctrl->E.regularize && Ctrl->E.mode != F_IS_DRIFT_D && Ctrl->E.mode != F_IS_DRIFT_T) { + GMT_Report (API, GMT_MSG_WARNING, "Option -E: +r is only implemented for the d and t models; ignored for this model\n"); + Ctrl->E.regularize = false; + } return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR); } @@ -306,7 +334,7 @@ EXTERN_MSC int GMT_x2sys_solve (void *V_API, int mode, void *args) { int min_ID, max_ID; unsigned int rec_mode; uint64_t n_par = 0, n_in = 0, n, m, n_tracks = 0, n_active, n_constraints = 0; - uint64_t i, p, j, k, r, s, row_off, row, n_COE = 0, w_col, id_col, bin_expect, *R = NULL, *col_off = NULL, *cluster = NULL; + uint64_t i, p, j, k, r, s, row_off, row, n_COE = 0, w_col, id_col, bin_expect, *R = NULL, *col_off = NULL, *cluster = NULL, *n_cross = NULL; size_t n_alloc = GMT_INITIAL_MEM_ROW_ALLOC, n_alloc_t = GMT_CHUNK; double *N = NULL, *a = NULL, *b = NULL, *in = NULL, *data[N_COE_PARS], sgn, old_mean, new_mean, sw2, C_i, C_j; @@ -645,6 +673,8 @@ EXTERN_MSC int GMT_x2sys_solve (void *V_API, int mode, void *args) { R[i]++; /* Increase COE count for track i */ R[j]++; /* Increase COE count for track j */ } + n_cross = gmt_M_memory (GMT, NULL, n_tracks, uint64_t); /* Raw crossing count, before R[] is capped */ + for (p = 0; p < n_tracks; p++) n_cross[p] = R[p]; for (p = n = 0; p < n_tracks; p++) { /* For each track, determine R[track], total number of parameters, and the column offsets */ (GMT->common.b.active[GMT_IN]) ? sprintf (trk[0], "%" PRIu64, p) : sprintf (trk[0], "%s", trk_list[p]); if (R[p] < n_par) /* Came up short */ @@ -778,6 +808,17 @@ EXTERN_MSC int GMT_x2sys_solve (void *V_API, int mode, void *args) { } } + if (Ctrl->E.regularize && n > 0) { + /* Only the real unknowns get lambda; the Lagrange rows must stay exact equality constraints. + * Scaling lambda to the mean diagonal keeps it invariant to the COE and weight units in use */ + double mean_diag = 0.0, lambda; + for (i = 0; i < n; i++) mean_diag += N[i*m+i]; + mean_diag /= n; + lambda = mean_diag / (Ctrl->E.K * Ctrl->E.K); + for (i = 0; i < n; i++) N[i*m+i] += lambda; + GMT_Report (API, GMT_MSG_INFORMATION, "Ridge-regularizing the solve: mean diagonal = %g, K = %g, lambda = %g\n", mean_diag, Ctrl->E.K, lambda); + } + /* Get LS solution */ if ((error = gmt_gaussjordan (GMT, N, (unsigned int)m, b)) != 0) { @@ -838,12 +879,24 @@ EXTERN_MSC int GMT_x2sys_solve (void *V_API, int mode, void *args) { sprintf (frmt_name, "%%-%ds", max_len+2); for (p = 0; p < n_tracks; p++) { - if (normalize) a[col_off[p]+1] /= range; /* Unnormalize slopes */ + if (normalize && R[p] > 1) a[col_off[p]+1] /= range; /* Unnormalize slopes; a track short of crossings has no slope, and +1 would be the next track's offset */ (GMT->common.b.active[GMT_IN]) ? sprintf (line, "%" PRIu64, p) : sprintf (line, frmt_name, trk_list[p]); strcat (line, "\t"); strcat (line, Ctrl->C.col); gmt_M_memset (var, N_BASIS, double); /* Reset all parameters to zero */ for (r = 0; r < R[p]; r++) var[r] = a[col_off[p]+r]; /* Just get the first R(p) items; the rest are set to 0 */ + if (Ctrl->E.mode != F_IS_SCALE && old_stdev > 0.0 && fabs (var[0]) > X2SYS_SOLVE_UNSTABLE_RATIO * old_stdev) { + (GMT->common.b.active[GMT_IN]) ? sprintf (trk[0], "%" PRIu64, p) : sprintf (trk[0], "%s", trk_list[p]); + GMT_Report (API, GMT_MSG_WARNING, + "Track %s: solved offset %.4g is %.0fx the pre-correction COE st.dev. (%.4g), backed by only %" PRIu64 " crossing(s) -- " + "likely unstable/unreliable, not a real systematic offset. %s\n", + trk[0], var[0], fabs (var[0]) / old_stdev, old_stdev, n_cross[p], + Ctrl->E.regularize ? + "Still unstable even with +r; try a smaller K, or fall back to -Ec." : + ((Ctrl->E.mode == F_IS_DRIFT_D || Ctrl->E.mode == F_IS_DRIFT_T) ? + "Try appending +r to this -E option, or fall back to -Ec." : + "Consider -Ec instead, or improve this track's crossover coverage.")); + } switch (Ctrl->E.mode) { /* Set up pointers to basis functions and assign constants */ case F_IS_CONSTANT: sprintf (text, "\t%g", var[0]); @@ -883,6 +936,7 @@ EXTERN_MSC int GMT_x2sys_solve (void *V_API, int mode, void *args) { gmt_M_free (GMT, b); gmt_M_free (GMT, R); gmt_M_free (GMT, col_off); + gmt_M_free (GMT, n_cross); if (!GMT->common.b.active[GMT_IN]) x2sys_free_list (GMT, trk_list, n_tracks); gmt_M_free (GMT, Out);