MCPcopy Create free account
hub / github.com/GenericMappingTools/gmt / gmtmath_LMSSCL

Function gmtmath_LMSSCL

src/gmtmath.c:3202–3250  ·  view source on GitHub ↗

Source from the content-addressed store, hash-verified

3200}
3201
3202GMT_LOCAL int gmtmath_LMSSCL (struct GMT_CTRL *GMT, struct GMTMATH_INFO *info, struct GMTMATH_STACK *S[], unsigned int last, unsigned int col)
3203/*OPERATOR: LMSSCL 1 1 LMS scale estimate (LMS STD) of A. */
3204{
3205 uint64_t s, row, k;
3206 unsigned int gmt_mode_selection = 0, GMT_n_multiples = 0;
3207 double lmsscl, mode, *z = NULL;
3208 struct GMT_DATATABLE *T = S[last]->D->table[0];
3209
3210 if (S[last]->constant) { /* Trivial case */
3211 for (s = 0; s < info->T->n_segments; s++) gmt_M_memset (T->segment[s]->data[col], info->T->segment[s]->n_rows, double);
3212 return 0;
3213 }
3214
3215 if (!info->local) z = gmt_M_memory (GMT, NULL, info->T->n_records, double);
3216
3217 for (s = k = 0; s < info->T->n_segments; s++) {
3218 if (info->local) {
3219 gmt_sort_array (GMT, T->segment[s]->data[col], info->T->segment[s]->n_rows, GMT_DOUBLE);
3220 for (row = info->T->segment[s]->n_rows; row > 1 && gmt_M_is_dnan (T->segment[s]->data[col][row-1]); row--);
3221 if (row) {
3222 gmt_mode (GMT, T->segment[s]->data[col], row, row/2, 0, gmt_mode_selection, &GMT_n_multiples, &mode);
3223 gmt_getmad (GMT, T->segment[s]->data[col], row, mode, &lmsscl);
3224 }
3225 else
3226 lmsscl = GMT->session.d_NaN;
3227
3228 for (row = 0; row < info->T->segment[s]->n_rows; row++) T->segment[s]->data[col][row] = lmsscl;
3229 if (GMT_n_multiples > 0) GMT_Report (GMT->parent, GMT_MSG_WARNING, "%d Multiple modes found for segment %" PRIu64 "\n", GMT_n_multiples, s);
3230 }
3231 else { /* Just accumulate the total table */
3232 gmt_M_memcpy (&z[k], T->segment[s]->data[col], info->T->segment[s]->n_rows, double);
3233 k += info->T->segment[s]->n_rows;
3234 }
3235 }
3236 if (info->local) return 0; /* Done with local */
3237 gmt_sort_array (GMT, z, info->T->n_records, GMT_DOUBLE);
3238 for (row = info->T->n_records; row > 1 && gmt_M_is_dnan (z[row-1]); row--);
3239 if (row) {
3240 gmt_mode (GMT, z, row, row/2, 0, gmt_mode_selection, &GMT_n_multiples, &mode);
3241 gmt_getmad (GMT, z, row, mode, &lmsscl);
3242 }
3243 else
3244 lmsscl = GMT->session.d_NaN;
3245
3246 for (s = 0; s < info->T->n_segments; s++) for (row = 0; row < info->T->segment[s]->n_rows; row++) T->segment[s]->data[col][row] = lmsscl;
3247 if (GMT_n_multiples > 0) GMT_Report (GMT->parent, GMT_MSG_WARNING, "%d Multiple modes found\n", GMT_n_multiples);
3248 gmt_M_free (GMT, z);
3249 return 0;
3250}
3251
3252GMT_LOCAL int gmtmath_LMSSCLW (struct GMT_CTRL *GMT, struct GMTMATH_INFO *info, struct GMTMATH_STACK *S[], unsigned int last, unsigned int col) {
3253/*OPERATOR: LMSSCLW 1 1 Weighted LMS scale estimate (LMS STD) of A for weights in B. */

Callers

nothing calls this directly

Calls 4

gmt_sort_arrayFunction · 0.85
gmt_modeFunction · 0.85
gmt_getmadFunction · 0.85
GMT_ReportFunction · 0.85

Tested by

no test coverage detected