Function gmtmath_MAD

src/gmtmath.c:3455–3498

}

GMT_LOCAL int gmtmath_MAD (struct GMT_CTRL *GMT, struct GMTMATH_INFO *info, struct GMTMATH_STACK *S[], unsigned int last, unsigned int col)
/*OPERATOR: MAD 1 1 Median Absolute Deviation (L1 STD) of A.  */
{
	uint64_t s, row, k;
	double mad, med, *z = NULL;
	struct GMT_DATATABLE *T = S[last]->D->table[0];

	if (S[last]->constant) {	/* Trivial case */
		for (s = 0; s < info->T->n_segments; s++) gmt_M_memset (T->segment[s]->data[col], info->T->segment[s]->n_rows, double);
		return 0;
	}

	if (!info->local) z = gmt_M_memory (GMT, NULL, info->T->n_records, double);

	for (s = k = 0; s < info->T->n_segments; s++) {
		if (info->local) {
			gmt_sort_array (GMT, T->segment[s]->data[col], info->T->segment[s]->n_rows, GMT_DOUBLE);
			for (row = info->T->segment[s]->n_rows; row > 1 && gmt_M_is_dnan (T->segment[s]->data[col][row-1]); row--);
			if (row) {
				med = (row%2) ? T->segment[s]->data[col][row/2] : 0.5 * (T->segment[s]->data[col][(row-1)/2] + T->segment[s]->data[col][row/2]);
				gmt_getmad (GMT, T->segment[s]->data[col], row, med, &mad);
			}
			else
				mad = GMT->session.d_NaN;
			for (row = 0; row < info->T->segment[s]->n_rows; row++) T->segment[s]->data[col][row] = mad;
		}
		else {	/* Just accumulate the total table */
			gmt_M_memcpy (&z[k], T->segment[s]->data[col], info->T->segment[s]->n_rows, double);
			k += info->T->segment[s]->n_rows;
		}
	}
	if (info->local) return 0;	/* Done with local */
	gmt_sort_array (GMT, z, info->T->n_records, GMT_DOUBLE);
	for (row = info->T->n_records; row > 1 && gmt_M_is_dnan (z[row-1]); row--);
	if (row) {
		med = (row%2) ? z[row/2] : 0.5 * (z[(row-1)/2] + z[row/2]);
		gmt_getmad (GMT, z, row, med, &mad);
	}
	else
		mad = GMT->session.d_NaN;
	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] = mad;
	gmt_M_free (GMT, z);
	return 0;
}

GMT_LOCAL int gmtmath_MADW (struct GMT_CTRL *GMT, struct GMTMATH_INFO *info, struct GMTMATH_STACK *S[], unsigned int last, unsigned int col) {
/*OPERATOR: MADW 2 1 Weighted Median Absolute Deviation (L1 STD) of A for weights in B.  */