-
Notifications
You must be signed in to change notification settings - Fork 375
/
Copy pathgmtmath.c
7080 lines (6379 loc) · 326 KB
/
gmtmath.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
/*--------------------------------------------------------------------
*
* Copyright (c) 1991-2025 by the GMT Team (https://www.generic-mapping-tools.org/team.html)
* See LICENSE.TXT file for copying and redistribution conditions.
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published by
* the Free Software Foundation; version 3 or any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Lesser General Public License for more details.
*
* Contact info: www.generic-mapping-tools.org
*--------------------------------------------------------------------*/
/*
* API functions to support the gmtmath application.
*
* Author: Paul Wessel
* Date: 1-JAN-2010
* Version: 6 API
*
* Brief synopsis: gmtmath.c is a reverse polish calculator that operates
* on table files (and constants) and perform basic mathematical operations
* on them like add, multiply, etc.
* Some operators only work on one operand (e.g., log, exp)
*
* Note on KEYS: AD(= means -A takes an input Dataset as argument which may be followed by optional modifiers.
*/
#include "gmt_dev.h"
#include "longopt/gmtmath_inc.h"
#define THIS_MODULE_CLASSIC_NAME "gmtmath"
#define THIS_MODULE_MODERN_NAME "gmtmath"
#define THIS_MODULE_LIB "core"
#define THIS_MODULE_PURPOSE "Reverse Polish Notation (RPN) calculator for data tables"
#define THIS_MODULE_KEYS "<D(,AD(=,TD(,>D}"
#define THIS_MODULE_NEEDS ""
#define THIS_MODULE_OPTIONS "-:>Vbdefghioqsw" GMT_OPT("HMm")
#define SPECIFIC_OPTIONS "AEILNQST" /* All non-common options except for -C which we will actually process in the loop over args */
#define GMTMATH_ARG_IS_OPERATOR 0
#define GMTMATH_ARG_IS_FILE -1
#define GMTMATH_ARG_IS_NUMBER -2
#define GMTMATH_ARG_IS_PI -3
#define GMTMATH_ARG_IS_E -4
#define GMTMATH_ARG_IS_F_EPS -5
#define GMTMATH_ARG_IS_D_EPS -6
#define GMTMATH_ARG_IS_EULER -7
#define GMTMATH_ARG_IS_PHI -8
#define GMTMATH_ARG_IS_TMIN -9
#define GMTMATH_ARG_IS_TMAX -10
#define GMTMATH_ARG_IS_TRANGE -11
#define GMTMATH_ARG_IS_TINC -12
#define GMTMATH_ARG_IS_N -13
#define GMTMATH_ARG_IS_J_MATRIX -14
#define GMTMATH_ARG_IS_T_MATRIX -15
#define GMTMATH_ARG_IS_t_MATRIX -16
#define GMTMATH_ARG_IS_STORE -50
#define GMTMATH_ARG_IS_RECALL -51
#define GMTMATH_ARG_IS_CLEAR -52
#define GMTMATH_ARG_IS_BAD -99
#define GMTMATH_STACK_SIZE 100
#define GMTMATH_STORE_SIZE 100
#define GMTMATH_STORE_CMD "STO@"
#define GMTMATH_RECALL_CMD "RCL@"
#define GMTMATH_CLEAR_CMD "CLR@"
#define COL_T 0 /* These are the first 3 columns in the Time structure */
#define COL_TN 1
#define COL_TJ 2
#define GMTMATH_COEFFICIENTS 0
#define GMTMATH_EVALUATE 1
#define GMTMATH_WEIGHTS 1
#define GMTMATH_SIGMAS 2
#define DOUBLE_BIT_MASK (~(1023ULL << 54ULL)) /* This will be 00000000 00111111 11111111 .... and sets to 0 anything larger than 2^53 which is max integer in double */
struct GMTMATH_CTRL { /* All control options for this program (except common args) */
/* active is true if the option has been activated */
struct GMTMATH_Out { /* = <filename> */
bool active;
char *file;
} Out;
struct GMTMATH_A { /* -A[-]<t_f(t).d>[+e][+w|s] */
bool active;
bool null;
unsigned int e_mode; /* 0 save coefficients, 1 save predictions and residuals */
unsigned int w_mode; /* 0 no weights, 1 = got weights, 2 = got sigmas */
char *file;
} A;
struct GMTMATH_C { /* -C<cols> */
bool active;
bool *cols;
} C;
struct GMTMATH_E { /* -E<min_eigenvalue> */
bool active;
double eigen;
} E;
struct GMTMATH_I { /* -I */
bool active;
} I;
struct GMTMATH_L { /* -L */
bool active;
} L;
struct GMTMATH_N { /* -N<n_col>/<t_col> */
bool active;
uint64_t ncol, tcol;
} N;
struct GMTMATH_Q { /* -Q[c|i|p|n] */
bool active;
int unit;
} Q;
struct GMTMATH_S { /* -S[f|l] */
bool active;
int mode; /* -1 or +1 */
} S;
struct GMTMATH_T { /* -T[<tmin/tmax/t_inc>[+]] | -T<file> */
bool active;
bool notime;
struct GMT_ARRAY T;
} T;
};
struct GMTMATH_INFO {
bool irregular; /* true if t_inc varies */
bool roots_found; /* true if roots have been solved for */
bool local; /* Per segment operation (true) or global operation (false) */
bool notime; /* No time-array available for operators who depend on that */
bool scalar; /* -Q in effect */
unsigned int n_roots; /* Number of roots found */
unsigned int fit_mode; /* Used for {LSQ|SVD}FIT */
unsigned int w_mode; /* Used for weighted fit */
uint64_t r_col; /* The column used to find roots */
uint64_t n_col; /* Number of columns */
double t_min, t_max, t_inc;
struct GMT_DATATABLE *T; /* Table with all time information */
struct GMT_ORDER **Q; /* For sorting columns */
};
struct GMTMATH_STACK {
struct GMT_DATASET *D; /* The dataset */
bool constant; /* true if a constant (see factor) and S == NULL */
double factor; /* The value if constant is true */
};
struct GMTMATH_STORED {
char *label; /* Name of this stored memory */
struct GMTMATH_STACK stored;
};
static void *New_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */
struct GMTMATH_CTRL *C = gmt_M_memory (GMT, NULL, 1, struct GMTMATH_CTRL);
/* Initialize values whose defaults are not 0/false/NULL */
C->C.cols = gmt_M_memory (GMT, NULL, GMT_MAX_COLUMNS, bool);
C->E.eigen = 1e-7; /* Default cutoff of small eigenvalues */
C->N.ncol = 2;
C->Q.unit = GMT->current.setting.proj_length_unit; /* Default output unit conversion for -Q */
return (C);
}
static void Free_Ctrl (struct GMT_CTRL *GMT, struct GMTMATH_CTRL *C) { /* Deallocate control structure */
if (!C) return;
gmt_M_str_free (C->Out.file);
gmt_M_str_free (C->A.file);
gmt_M_free (GMT, C->C.cols);
gmt_free_array (GMT, &(C->T.T));
gmt_M_free (GMT, C);
}
GMT_LOCAL bool gmtmath_decode_columns (char *txt, bool *skip, uint64_t n_col, uint64_t t_col) {
uint64_t i, start, stop;
unsigned int pos;
char p[GMT_BUFSIZ];
/* gmtmath_decode_columns is used to handle the parsing of -C<cols>. */
if (!txt[0]) { /* Reset to default */
for (i = 0; i < n_col; i++) skip[i] = false;
skip[t_col] = true;
}
else if (txt[0] == 'r' && txt[1] == '\0') { /* Reverse all settings */
for (i = 0; i < n_col; i++) skip[i] = !skip[i];
}
else if (txt[0] == 'a') { /* Select all columns */
for (i = 0; i < n_col; i++) skip[i] = false;
}
else { /* Set the selected columns */
for (i = 0; i < n_col; i++) skip[i] = true;
pos = 0;
while ((gmt_strtok (txt, ",", &pos, p))) {
if (strchr (p, '-'))
sscanf (p, "%" PRIu64 "-%" PRIu64, &start, &stop);
else if (strchr (p, 'x')) /* -Cx is the x/lon column */
start = stop = GMT_X;
else if (strchr (p, 'y')) /* -Cy is the y/lat column */
start = stop = GMT_Y;
else {
sscanf (p, "%" PRIu64, &start);
stop = start;
}
stop = MIN (stop, n_col-1);
for (i = start; i <= stop; i++) skip[i] = false;
}
}
return (!skip[t_col]); /* Returns true if we are changing the time column */
}
GMT_LOCAL int gmtmath_find_stored_item (struct GMTMATH_STORED *recall[], int n_stored, char *label) {
/* Linear search to find the named storage item */
int k = 0;
while (k < n_stored && strcmp (recall[k]->label, label)) k++;
return (k == n_stored ? GMT_NOTSET : k);
}
GMT_LOCAL int gmtmath_load_column (struct GMT_DATASET *to, uint64_t to_col, struct GMT_DATATABLE *from, uint64_t from_col) {
/* Copies data from one column to another */
uint64_t seg;
for (seg = 0; seg < from->n_segments; seg++) {
if (to->table[0]->segment[seg]->n_rows == from->segment[seg]->n_rows)
gmt_M_memcpy (to->table[0]->segment[seg]->data[to_col], from->segment[seg]->data[from_col], from->segment[seg]->n_rows, double);
else
return GMT_NOTSET;
}
return GMT_NOERROR;
}
GMT_LOCAL void gmtmath_load_text (struct GMT_DATASET *to, struct GMT_DATASET *from) {
/* Copies trailing text from all input segments to all output segments */
uint64_t tbl, seg, row;
char **Tf = NULL, **Tt = NULL;
struct GMT_DATASEGMENT_HIDDEN *SH = NULL;
if (from->type != GMT_READ_MIXED) return; /* No text in this one */
for (tbl = 0; tbl < from->n_tables; tbl++) {
for (seg = 0; seg < from->table[tbl]->n_segments; seg++) {
struct GMT_DATASEGMENT *Sf = from->table[tbl]->segment[seg];
struct GMT_DATASEGMENT *St = to->table[tbl]->segment[seg];
if ((Tf = Sf->text) == NULL) continue; /* No text in this segment */
if ((Tt = St->text) == NULL) continue; /* Safety valve for duplicating in to NULL array */
for (row = 0; row < Sf->n_rows; row++)
if (!Tt[row]) Tt[row] = strdup (Tf[row]); /* Only duplicate if not set already */
SH = gmt_get_DS_hidden (St);
SH->alloc_mode_text = GMT_ALLOC_INTERNALLY; /* So we may delete the text we allocated */
}
}
}
/* ---------------------- start convenience functions --------------------- */
GMT_LOCAL int gmtmath_solve_LS_system (struct GMT_CTRL *GMT, struct GMTMATH_INFO *info, struct GMTMATH_STACK *S, uint64_t n_col, bool skip[], char *file, bool svd, double eigen_min, struct GMT_OPTION *options, struct GMT_DATASET *A) {
/* Consider the current table the augmented matrix [A | b], making up the linear system Ax = b.
* We will set up the normal equations, solve for x, and output the solution before quitting.
* This function is special since it operates across columns and returns n_col scalars.
* We try to solve this positive definite & symmetric matrix with Cholesky methods; if that fails
* we do a full SVD decomposition and set small eigenvalues to zero, yielding an approximate solution.
* However, if svd == true then we do the full SVD.
*/
unsigned int i, j, k, k0, i2, j2, n;
int ier = 0;
uint64_t row, seg, rhs, w_col = 0;
double cond, w = 1.0;
double *N = NULL, *r = NULL, *d = NULL, *x = NULL;
struct GMT_DATATABLE *T = S->D->table[0];
struct GMT_DATASET *D = NULL;
gmt_M_unused(A);
for (i = n = 0; i < n_col; i++) if (!skip[i]) n++; /* Need to find how many active columns we have */
if (n < 2) {
char *pre[2] = {"LSQ", "SVD"};
GMT_Report (GMT->parent, GMT_MSG_ERROR, "%sFIT requires at least 2 active columns!\n", pre[svd]);
return (GMT_DIM_TOO_SMALL);
}
rhs = n_col - 1;
while (rhs > 0 && skip[rhs]) rhs--; /* Get last active col number as the rhs vector b */
n--; /* Account for b, the rhs vector, to get row & col dimensions of normal matrix N */
if (info->w_mode) {
w_col = rhs - 1; /* If there are weights, this is the column they are stored in */
while (w_col > 0 && skip[w_col]) w_col--; /* Get next active col number as the weight vector w */
n--; /* Account for w, the rhs vector, to get row & col dimensions of normal matrix N */
}
N = gmt_M_memory (GMT, NULL, n*n, double);
r = gmt_M_memory (GMT, NULL, T->n_records, double);
#if 0
fprintf (stderr, "Printout of A | b matrix\n");
fprintf (stderr, "------------------------------------------------------------------\n");
for (seg = 0; seg < info->T->n_segments; seg++) {
for (row = 0; row < T->segment[seg]->n_rows; row++) {
for (j = 0; j < rhs; j++) {
if (skip[j]) continue;
fprintf (stderr, "%g\t", T->segment[seg]->data[j][row]);
}
fprintf (stderr, "|\t%g\n", T->segment[seg]->data[rhs][row]);
}
}
fprintf (stderr, "------------------------------------------------------------------\n");
#endif
/* Here we build A^T*W*A*x = A^T*W*b ==> N*x = r, where W is the diagonal matrix with squared weights w */
/* Do the row & col dot products, skipping inactive columns as we go along */
for (j = j2 = 0; j < n; j2++) { /* j2 is table column, j is row in N matrix */
if (skip[j2]) continue;
for (i = i2 = 0; i < n; i2++) { /* i2 is table column, i is column in N matrix */
if (skip[i2]) continue;
k0 = j * n + i;
N[k0] = 0.0;
for (seg = k = 0; seg < info->T->n_segments; seg++) {
for (row = 0; row < T->segment[seg]->n_rows; row++, k++) {
if (info->w_mode) {
w = pow (T->segment[seg]->data[w_col][row], 2.0);
if (info->w_mode == GMTMATH_SIGMAS) w = 1.0 / w; /* Got sigma */
}
N[k0] += w * T->segment[seg]->data[j2][row] * T->segment[seg]->data[i2][row];
}
}
i++;
}
r[j] = 0.0;
for (seg = k = 0; seg < info->T->n_segments; seg++) {
for (row = 0; row < T->segment[seg]->n_rows; row++, k++) {
if (info->w_mode) {
w = pow (T->segment[seg]->data[w_col][row], 2.0);
if (info->w_mode == GMTMATH_SIGMAS) w = 1.0 / w; /* Got sigma */
}
r[j] += w * T->segment[seg]->data[j2][row] * T->segment[seg]->data[rhs][row];
}
}
j++;
}
if (svd)
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Solve LS system via SVD decomposition and exclude eigenvalues < %g.\n", eigen_min);
else
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Solve LS system via Cholesky decomposition\n");
#if 0
fprintf (stderr, "Printout of N and r matrix\n");
fprintf (stderr, "------------------------------------------------------------------\n");
for (j = 0; j < n; j++) {
for (k = 0; k < n; k++)
fprintf (stderr, "%g\t", N[j*n+k]);
fprintf (stderr, "\n");
}
for (k = 0; k < n; k++)
fprintf (stderr, "%g\n", r[k]);
fprintf (stderr, "------------------------------------------------------------------\n");
#endif
d = gmt_M_memory (GMT, NULL, n, double);
x = gmt_M_memory (GMT, NULL, n, double);
if (svd || ((ier = gmt_chol_dcmp (GMT, N, d, &cond, n, n) ) != 0)) { /* Cholesky decomposition failed, use SVD method, or use SVD if specified */
unsigned int nrots;
double *b = NULL, *z = NULL, *v = NULL, *lambda = NULL;
if (!svd) {
GMT_Report (GMT->parent, GMT_MSG_WARNING,
"Cholesky decomposition failed, try SVD decomposition instead and exclude eigenvalues < %g.\n", eigen_min);
gmt_chol_recover (GMT, N, d, n, n, ier, true); /* Restore to former matrix N */
}
/* Solve instead using the SVD of a square matrix via gmt_jacobi */
lambda = gmt_M_memory (GMT, NULL, n, double);
b = gmt_M_memory (GMT, NULL, n, double);
z = gmt_M_memory (GMT, NULL, n, double);
v = gmt_M_memory (GMT, NULL, n*n, double);
if (gmt_jacobi (GMT, N, n, n, lambda, v, b, z, &nrots)) {
GMT_Report (GMT->parent, GMT_MSG_ERROR, "Eigenvalue routine failed to converge in 50 sweeps.\n");
GMT_Report (GMT->parent, GMT_MSG_ERROR, "The solution might be inaccurate.\n");
}
/* Solution x = v * lambda^-1 * v' * r */
/* First do d = V' * r, so x = v * lambda^-1 * d */
for (j = 0; j < n; j++) for (k = 0, d[j] = 0.0; k < n; k++) d[j] += v[j*n+k] * r[k];
/* Then do d = lambda^-1 * d by setting small lambda's to zero */
for (j = k = 0; j < n; j++) {
if (lambda[j] < eigen_min) {
d[j] = 0.0;
k++;
}
else
d[j] /= lambda[j];
}
if (k) GMT_Report (GMT->parent, GMT_MSG_WARNING, "%d eigenvalues < %g set to zero to yield a stable solution\n", k, eigen_min);
/* Finally do x = v * d */
for (j = 0; j < n; j++) for (k = 0; k < n; k++) x[j] += v[k*n+j] * d[k];
gmt_M_free (GMT, b);
gmt_M_free (GMT, z);
gmt_M_free (GMT, v);
gmt_M_free (GMT, lambda);
}
else { /* Decomposition worked, now solve system */
gmt_chol_solv (GMT, N, x, r, n, n);
}
if (info->fit_mode == GMTMATH_COEFFICIENTS) { /* Return coefficients only as a single vector */
uint64_t dim[GMT_DIM_SIZE] = {1, 1, 0, 1};
dim[GMT_ROW] = n;
if ((D = GMT_Create_Data (GMT->parent, GMT_IS_DATASET, GMT_IS_NONE, 0, dim, NULL, NULL, 0, 0, NULL)) == NULL)
return (GMT->parent->error);
for (k = 0; k < n; k++) D->table[0]->segment[0]->data[GMT_X][k] = x[k];
D->table[0]->segment[0]->n_rows = n;
GMT_Set_Comment (GMT->parent, GMT_IS_DATASET, GMT_COMMENT_IS_OPTION | GMT_COMMENT_IS_COMMAND, options, D);
if (GMT->common.h.add_colnames) {
char header[GMT_LEN16] = {""};
sprintf (header, "#coefficients");
if (GMT_Set_Comment (GMT->parent, GMT_IS_DATASET, GMT_COMMENT_IS_COLNAMES, header, D)) return (GMT->parent->error);
}
if (GMT_Write_Data (GMT->parent, GMT_IS_DATASET, (file ? GMT_IS_FILE : GMT_IS_STREAM), GMT_IS_NONE, GMT_WRITE_NORMAL, NULL, file, D) != GMT_NOERROR)
return (GMT->parent->error);
}
else { /* Return t, y, p(t), r(t), where p(t) is the predicted solution and r(t) is the residuals */
double value;
struct GMT_DATASET_HIDDEN *DH = gmt_get_DD_hidden (S->D);
k = (unsigned int)DH->dim[GMT_COL];
DH->dim[GMT_COL] = (info->w_mode) ? 5 : 4; /* State we want a different set of columns on output */
D = GMT_Duplicate_Data (GMT->parent, GMT_IS_DATASET, GMT_DUPLICATE_ALLOC, S->D); /* Same table length as S->D, but with up to n_cols columns (lon, lat, dist, g1, g2, ...) */
DH->dim[GMT_COL] = k; /* Reset the original columns */
if (D->table[0]->n_segments > 1) gmt_set_segmentheader (GMT, GMT_OUT, true); /* More than one segment triggers -mo */
gmtmath_load_column (D, 0, info->T, COL_T); /* Place the time-column in first output column */
for (seg = k = 0; seg < info->T->n_segments; seg++) {
for (row = 0; row < T->segment[seg]->n_rows; row++, k++) {
D->table[0]->segment[seg]->data[1][row] = T->segment[seg]->data[rhs][row];
value = 0.0;
for (j2 = j = 0; j2 < n; j2++) { /* j2 is table column, j is entry in x vector */
if (skip[j2]) continue; /* Not included in the fit */
value += T->segment[seg]->data[j2][row] * x[j]; /* Sum up the solution */
j++;
}
D->table[0]->segment[seg]->data[2][row] = value;
D->table[0]->segment[seg]->data[3][row] = T->segment[seg]->data[rhs][row] - value;
if (info->w_mode) D->table[0]->segment[seg]->data[4][row] = T->segment[seg]->data[w_col][row];
}
}
if (GMT->common.h.add_colnames) {
char header[GMT_LEN128] = {""};
sprintf (header, "#t[0]\tobserved(t)[1]\tpredict(t)[2]\tresidual(t)[3]");
if (info->w_mode == GMTMATH_WEIGHTS) strcat (header, "\tweight(t)[4]");
else if (info->w_mode == GMTMATH_SIGMAS) strcat (header, "\tsigma(t)[4]");
if (GMT_Set_Comment (GMT->parent, GMT_IS_DATASET, GMT_COMMENT_IS_COLNAMES, header, D)) return (GMT->parent->error);
}
if (GMT_Write_Data (GMT->parent, GMT_IS_DATASET, (file ? GMT_IS_FILE : GMT_IS_STREAM), GMT_IS_NONE, GMT_WRITE_NORMAL, NULL, file, D) != GMT_NOERROR) {
return (GMT->parent->error);
}
}
gmt_M_free (GMT, x);
gmt_M_free (GMT, d);
gmt_M_free (GMT, N);
gmt_M_free (GMT, r);
return (GMT_NOERROR);
}
GMT_LOCAL int gmtmath_solve_LSQFIT (struct GMT_CTRL *GMT, struct GMTMATH_INFO *info, struct GMTMATH_STACK *S, uint64_t n_col, bool skip[], double eigen, char *file, struct GMT_OPTION *options, struct GMT_DATASET *A) {
return (gmtmath_solve_LS_system (GMT, info, S, n_col, skip, file, false, eigen, options, A));
}
GMT_LOCAL int gmtmath_solve_SVDFIT (struct GMT_CTRL *GMT, struct GMTMATH_INFO *info, struct GMTMATH_STACK *S, uint64_t n_col, bool skip[], double eigen, char *file, struct GMT_OPTION *options, struct GMT_DATASET *A) {
return (gmtmath_solve_LS_system (GMT, info, S, n_col, skip, file, true, eigen, options, A));
}
GMT_LOCAL void gmtmath_load_const_column (struct GMT_DATASET *to, uint64_t to_col, double factor) {
/* Sets all rows in a column to a constant factor */
uint64_t row, seg;
for (seg = 0; seg < to->n_segments; seg++) {
for (row = 0; row < to->table[0]->segment[seg]->n_rows; row++) to->table[0]->segment[seg]->data[to_col][row] = factor;
}
}
GMT_LOCAL bool gmtmath_same_size (struct GMT_DATASET *A, struct GMT_DATASET *B) {
/* Are the two dataset the same size */
uint64_t seg;
if (!(A->table[0]->n_segments == B->table[0]->n_segments && A->table[0]->n_columns == B->table[0]->n_columns)) return (false);
for (seg = 0; seg < A->table[0]->n_segments; seg++) if (A->table[0]->segment[seg]->n_rows != B->table[0]->segment[seg]->n_rows) return (false);
return (true);
}
GMT_LOCAL bool gmtmath_same_domain (struct GMT_DATASET *A, uint64_t t_col, struct GMT_DATATABLE *B) {
/* Are the two dataset the same domain */
uint64_t seg;
for (seg = 0; seg < A->table[0]->n_segments; seg++) {
if (!(doubleAlmostEqualZero (A->table[0]->min[t_col], B->min[COL_T])
&& doubleAlmostEqualZero (A->table[0]->max[t_col], B->max[COL_T])))
return (false);
}
return (true);
}
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 [-A<ftable>[+e][+r][+s|w]] [-C<cols>] [-E<eigen>] [-I] [-L] [-N<n_col>[/<t_col>]] [-Q[%s|n]] [-S[f|l]] "
"[-T[<file>|<list>|<min>/<max>/<inc>[+b|i|l|n]]] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] A B op C op ... = [<outfile>]\n",
name, GMT_DIM_UNITS_DISPLAY, GMT_V_OPT, GMT_b_OPT, GMT_d_OPT, GMT_e_OPT, GMT_f_OPT, GMT_g_OPT, GMT_h_OPT, GMT_i_OPT, GMT_o_OPT, GMT_q_OPT, GMT_s_OPT, GMT_w_OPT, GMT_PAR_OPT);
if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS);
GMT_Message (API, GMT_TIME_NONE, " REQUIRED ARGUMENTS:\n");
GMT_Usage (API, 1, "\n= <outfile>");
GMT_Usage (API, -2, "Writes the final stack result to the named file and pops it off the stack. "
"If no file is named then we write to standard output.");
GMT_Usage (API, 1, "\n<operands>");
GMT_Usage (API, -2, "A, B, etc. are table files, constants, or symbols (see below). "
"To read from standard input, give filename as STDIN (which can appear more than once). "
"The stack can hold up to %d entries (given enough memory).", GMTMATH_STACK_SIZE);
GMT_Usage (API, 1, "\n<operators>");
GMT_Usage (API, -2, "Trigonometric operators expect radians unless noted otherwise. "
"The operator names, the number of input and output arguments, and operation:\n");
/* Should we add an operator that has more characters that these, then the -20 will need to increase */
GMT_Message (API, GMT_TIME_NONE,
" Name #args Returns\n"
" -----------------------\n");
GMT_Message (API, GMT_TIME_NONE, " ABS 1 1 "); GMT_Usage (API, -21, "abs (A)");
GMT_Message (API, GMT_TIME_NONE, " ACOS 1 1 "); GMT_Usage (API, -21, "acos (A)");
GMT_Message (API, GMT_TIME_NONE, " ACOSD 1 1 "); GMT_Usage (API, -21, "acosd (A)");
GMT_Message (API, GMT_TIME_NONE, " ACOSH 1 1 "); GMT_Usage (API, -21, "acosh (A)");
GMT_Message (API, GMT_TIME_NONE, " ACOT 1 1 "); GMT_Usage (API, -21, "acot (A)");
GMT_Message (API, GMT_TIME_NONE, " ACOTD 1 1 "); GMT_Usage (API, -21, "acotd (A)");
GMT_Message (API, GMT_TIME_NONE, " ACOTH 1 1 "); GMT_Usage (API, -21, "acoth (A)");
GMT_Message (API, GMT_TIME_NONE, " ACSC 1 1 "); GMT_Usage (API, -21, "acsc (A)");
GMT_Message (API, GMT_TIME_NONE, " ACSCD 1 1 "); GMT_Usage (API, -21, "acscd (A)");
GMT_Message (API, GMT_TIME_NONE, " ACSCH 1 1 "); GMT_Usage (API, -21, "acsch (A)");
GMT_Message (API, GMT_TIME_NONE, " ADD 2 1 "); GMT_Usage (API, -21, "A + B");
GMT_Message (API, GMT_TIME_NONE, " AND 2 1 "); GMT_Usage (API, -21, "B if A == NaN, else A");
GMT_Message (API, GMT_TIME_NONE, " ASEC 1 1 "); GMT_Usage (API, -21, "asec (A)");
GMT_Message (API, GMT_TIME_NONE, " ASECD 1 1 "); GMT_Usage (API, -21, "asecd (A)");
GMT_Message (API, GMT_TIME_NONE, " ASECH 1 1 "); GMT_Usage (API, -21, "asech (A)");
GMT_Message (API, GMT_TIME_NONE, " ASIN 1 1 "); GMT_Usage (API, -21, "asin (A)");
GMT_Message (API, GMT_TIME_NONE, " ASIND 1 1 "); GMT_Usage (API, -21, "asind (A)");
GMT_Message (API, GMT_TIME_NONE, " ASINH 1 1 "); GMT_Usage (API, -21, "asinh (A)");
GMT_Message (API, GMT_TIME_NONE, " ATAN 1 1 "); GMT_Usage (API, -21, "atan (A)");
GMT_Message (API, GMT_TIME_NONE, " ATAND 1 1 "); GMT_Usage (API, -21, "atand (A)");
GMT_Message (API, GMT_TIME_NONE, " ATAN2 2 1 "); GMT_Usage (API, -21, "atan2 (A, B)");
GMT_Message (API, GMT_TIME_NONE, " ATAN2D 2 1 "); GMT_Usage (API, -21, "atan2d (A, B)");
GMT_Message (API, GMT_TIME_NONE, " ATANH 1 1 "); GMT_Usage (API, -21, "atanh (A)");
GMT_Message (API, GMT_TIME_NONE, " BCDF 3 1 "); GMT_Usage (API, -21, "Binomial cumulative distribution function for p = A, n = B and x = C");
GMT_Message (API, GMT_TIME_NONE, " BEI 1 1 "); GMT_Usage (API, -21, "Kelvin function bei (A)");
GMT_Message (API, GMT_TIME_NONE, " BER 1 1 "); GMT_Usage (API, -21, "Kelvin function ber (A)");
GMT_Message (API, GMT_TIME_NONE, " BPDF 3 1 "); GMT_Usage (API, -21, "Binomial probability density function for p = A, n = B and x = C");
GMT_Message (API, GMT_TIME_NONE, " BITAND 2 1 "); GMT_Usage (API, -21, "A & B (bitwise AND operator)");
GMT_Message (API, GMT_TIME_NONE, " BITLEFT 2 1 "); GMT_Usage (API, -21, "A << B (bitwise left-shift operator)");
GMT_Message (API, GMT_TIME_NONE, " BITNOT 1 1 "); GMT_Usage (API, -21, "~A (bitwise NOT operator, i.e., return two's complement)");
GMT_Message (API, GMT_TIME_NONE, " BITOR 2 1 "); GMT_Usage (API, -21, "A | B (bitwise OR operator)");
GMT_Message (API, GMT_TIME_NONE, " BITRIGHT 2 1 "); GMT_Usage (API, -21, "A >> B (bitwise right-shift operator)");
GMT_Message (API, GMT_TIME_NONE, " BITTEST 2 1 "); GMT_Usage (API, -21, "1 if bit B of A is set, else 0 (bitwise TEST operator)");
GMT_Message (API, GMT_TIME_NONE, " BITXOR 2 1 "); GMT_Usage (API, -21, "A ^ B (bitwise XOR operator)");
GMT_Message (API, GMT_TIME_NONE, " CEIL 1 1 "); GMT_Usage (API, -21, "ceil (A) (smallest integer >= A)");
GMT_Message (API, GMT_TIME_NONE, " CHI2CRIT 2 1 "); GMT_Usage (API, -21, "Chi-squared distribution critical value for alpha = A and nu = B");
GMT_Message (API, GMT_TIME_NONE, " CHI2CDF 2 1 "); GMT_Usage (API, -21, "Chi-squared cumulative distribution function for chi2 = A and nu = B");
GMT_Message (API, GMT_TIME_NONE, " CHI2PDF 2 1 "); GMT_Usage (API, -21, "Chi-squared probability density function for chi = A and nu = B");
GMT_Message (API, GMT_TIME_NONE, " COL 1 1 "); GMT_Usage (API, -21, "Places column A on the stack");
GMT_Message (API, GMT_TIME_NONE, " COMB 2 1 "); GMT_Usage (API, -21, "Combinations n_C_r, with n = A and r = B");
GMT_Message (API, GMT_TIME_NONE, " CORRCOEFF 2 1 "); GMT_Usage (API, -21, "Correlation coefficient r(A, B)");
GMT_Message (API, GMT_TIME_NONE, " COS 1 1 "); GMT_Usage (API, -21, "cos (A) (A in radians)");
GMT_Message (API, GMT_TIME_NONE, " COSD 1 1 "); GMT_Usage (API, -21, "cos (A) (A in degrees)");
GMT_Message (API, GMT_TIME_NONE, " COSH 1 1 "); GMT_Usage (API, -21, "cosh (A)");
GMT_Message (API, GMT_TIME_NONE, " COT 1 1 "); GMT_Usage (API, -21, "cot (A) (A in radians)");
GMT_Message (API, GMT_TIME_NONE, " COTD 1 1 "); GMT_Usage (API, -21, "cot (A) (A in degrees)");
GMT_Message (API, GMT_TIME_NONE, " COTH 1 1 "); GMT_Usage (API, -21, "coth (A)");
GMT_Message (API, GMT_TIME_NONE, " CSC 1 1 "); GMT_Usage (API, -21, "csc (A) (A in radians)");
GMT_Message (API, GMT_TIME_NONE, " CSCD 1 1 "); GMT_Usage (API, -21, "csc (A) (A in degrees)");
GMT_Message (API, GMT_TIME_NONE, " CSCH 1 1 "); GMT_Usage (API, -21, "csch (A)");
GMT_Message (API, GMT_TIME_NONE, " PCDF 2 1 "); GMT_Usage (API, -21, "Poisson cumulative distribution function for x = A and lambda = B");
GMT_Message (API, GMT_TIME_NONE, " DDT 1 1 "); GMT_Usage (API, -21, "d(A)/dt Central 1st derivative");
GMT_Message (API, GMT_TIME_NONE, " D2DT2 1 1 "); GMT_Usage (API, -21, "d^2(A)/dt^2 Central 2nd derivative");
GMT_Message (API, GMT_TIME_NONE, " D2R 1 1 "); GMT_Usage (API, -21, "Converts Degrees to Radians");
GMT_Message (API, GMT_TIME_NONE, " DEG2KM 1 1 "); GMT_Usage (API, -21, "Converts Spherical Degrees to Kilometers");
GMT_Message (API, GMT_TIME_NONE, " DENAN 2 1 "); GMT_Usage (API, -21, "Replace NaNs in A with values from B");
GMT_Message (API, GMT_TIME_NONE, " DILOG 1 1 "); GMT_Usage (API, -21, "dilog (A)");
GMT_Message (API, GMT_TIME_NONE, " DIFF 1 1 "); GMT_Usage (API, -21, "Forward difference between elements of A (A[1]-A[0], A[2]-A[1], ..., NaN)");
GMT_Message (API, GMT_TIME_NONE, " DIV 2 1 "); GMT_Usage (API, -21, "A / B");
GMT_Message (API, GMT_TIME_NONE, " DUP 1 2 "); GMT_Usage (API, -21, "Places duplicate of A on the stack");
GMT_Message (API, GMT_TIME_NONE, " ECDF 2 1 "); GMT_Usage (API, -21, "Exponential cumulative distribution function for x = A and lambda = B");
GMT_Message (API, GMT_TIME_NONE, " ECRIT 2 1 "); GMT_Usage (API, -21, "Exponential distribution critical value for alpha = A and lambda = B");
GMT_Message (API, GMT_TIME_NONE, " EPDF 2 1 "); GMT_Usage (API, -21, "Exponential probability density function for x = A and lambda = B");
GMT_Message (API, GMT_TIME_NONE, " ERF 1 1 "); GMT_Usage (API, -21, "Error function erf (A)");
GMT_Message (API, GMT_TIME_NONE, " ERFC 1 1 "); GMT_Usage (API, -21, "Complementary Error function erfc (A)");
GMT_Message (API, GMT_TIME_NONE, " ERFINV 1 1 "); GMT_Usage (API, -21, "Inverse error function of A");
GMT_Message (API, GMT_TIME_NONE, " EQ 2 1 "); GMT_Usage (API, -21, "1 if A == B, else 0");
GMT_Message (API, GMT_TIME_NONE, " EXCH 2 2 "); GMT_Usage (API, -21, "Exchanges A and B on the stack");
GMT_Message (API, GMT_TIME_NONE, " EXP 1 1 "); GMT_Usage (API, -21, "exp (A)");
GMT_Message (API, GMT_TIME_NONE, " FACT 1 1 "); GMT_Usage (API, -21, "A! (A factorial)");
GMT_Message (API, GMT_TIME_NONE, " FCRIT 3 1 "); GMT_Usage (API, -21, "F distribution critical value for alpha = A, nu1 = B, and nu2 = C");
GMT_Message (API, GMT_TIME_NONE, " FCDF 3 1 "); GMT_Usage (API, -21, "F cumulative distribution function for F = A, nu1 = B, and nu2 = C");
GMT_Message (API, GMT_TIME_NONE, " FLIPUD 1 1 "); GMT_Usage (API, -21, "Reverse order of each column");
GMT_Message (API, GMT_TIME_NONE, " FLOOR 1 1 "); GMT_Usage (API, -21, "floor (A) (greatest integer <= A)");
GMT_Message (API, GMT_TIME_NONE, " FMOD 2 1 "); GMT_Usage (API, -21, "A % B (remainder after truncated division)");
GMT_Message (API, GMT_TIME_NONE, " FPDF 3 1 "); GMT_Usage (API, -21, "F probability density distribution for F = A, nu1 = B and nu2 = C");
GMT_Message (API, GMT_TIME_NONE, " GE 2 1 "); GMT_Usage (API, -21, "1 if A >= B, else 0");
GMT_Message (API, GMT_TIME_NONE, " GT 2 1 "); GMT_Usage (API, -21, "1 if A > B, else 0");
GMT_Message (API, GMT_TIME_NONE, " HSV2LAB 3 3 "); GMT_Usage (API, -21, "Convert hsv to lab, with h = A, s = B and v = C");
GMT_Message (API, GMT_TIME_NONE, " HSV2RGB 3 3 "); GMT_Usage (API, -21, "Convert hsv to rgb, with h = A, s = B and v = C");
GMT_Message (API, GMT_TIME_NONE, " HSV2XYZ 3 3 "); GMT_Usage (API, -21, "Convert hsv to xyz, with h = A, s = B and v = C");
GMT_Message (API, GMT_TIME_NONE, " HYPOT 2 1 "); GMT_Usage (API, -21, "hypot (A, B) = sqrt (A*A + B*B)");
GMT_Message (API, GMT_TIME_NONE, " I0 1 1 "); GMT_Usage (API, -21, "Modified Bessel function of A (1st kind, order 0)");
GMT_Message (API, GMT_TIME_NONE, " I1 1 1 "); GMT_Usage (API, -21, "Modified Bessel function of A (1st kind, order 1)");
GMT_Message (API, GMT_TIME_NONE, " IFELSE 3 1 "); GMT_Usage (API, -21, "B if A != 0, else C");
GMT_Message (API, GMT_TIME_NONE, " IN 2 1 "); GMT_Usage (API, -21, "Modified Bessel function of A (1st kind, order B)");
GMT_Message (API, GMT_TIME_NONE, " INRANGE 3 1 "); GMT_Usage (API, -21, "1 if B <= A <= C, else 0");
GMT_Message (API, GMT_TIME_NONE, " INT 1 1 "); GMT_Usage (API, -21, "Numerically integrate A");
GMT_Message (API, GMT_TIME_NONE, " INV 1 1 "); GMT_Usage (API, -21, "1 / A");
GMT_Message (API, GMT_TIME_NONE, " ISFINITE 1 1 "); GMT_Usage (API, -21, "1 if A is finite, else 0");
GMT_Message (API, GMT_TIME_NONE, " ISNAN 1 1 "); GMT_Usage (API, -21, "1 if A == NaN, else 0");
GMT_Message (API, GMT_TIME_NONE, " J0 1 1 "); GMT_Usage (API, -21, "Bessel function of A (1st kind, order 0)");
GMT_Message (API, GMT_TIME_NONE, " J1 1 1 "); GMT_Usage (API, -21, "Bessel function of A (1st kind, order 1)");
GMT_Message (API, GMT_TIME_NONE, " JN 2 1 "); GMT_Usage (API, -21, "Bessel function of A (1st kind, order B)");
GMT_Message (API, GMT_TIME_NONE, " K0 1 1 "); GMT_Usage (API, -21, "Modified Kelvin function of A (2nd kind, order 0)");
GMT_Message (API, GMT_TIME_NONE, " K1 1 1 "); GMT_Usage (API, -21, "Modified Bessel function of A (2nd kind, order 1)");
GMT_Message (API, GMT_TIME_NONE, " KM2DEG 1 1 "); GMT_Usage (API, -21, "Converts Kilometers to Spherical Degrees");
GMT_Message (API, GMT_TIME_NONE, " KN 2 1 "); GMT_Usage (API, -21, "Modified Bessel function of A (2nd kind, order B)");
GMT_Message (API, GMT_TIME_NONE, " KEI 1 1 "); GMT_Usage (API, -21, "Kelvin function kei (A)");
GMT_Message (API, GMT_TIME_NONE, " KER 1 1 "); GMT_Usage (API, -21, "Kelvin function ker (A)");
GMT_Message (API, GMT_TIME_NONE, " KURT 1 1 "); GMT_Usage (API, -21, "Kurtosis of A");
GMT_Message (API, GMT_TIME_NONE, " LAB2HSV 3 3 "); GMT_Usage (API, -21, "Convert lab to hsv, with l = A, a = B and b = C");
GMT_Message (API, GMT_TIME_NONE, " LAB2RGB 3 3 "); GMT_Usage (API, -21, "Convert lab to rgb, with l = A, a = B and b = C");
GMT_Message (API, GMT_TIME_NONE, " LAB2XYZ 3 3 "); GMT_Usage (API, -21, "Convert lab to xyz, with l = A, a = B and b = C");
GMT_Message (API, GMT_TIME_NONE, " LCDF 1 1 "); GMT_Usage (API, -21, "Laplace cumulative distribution function for z = A");
GMT_Message (API, GMT_TIME_NONE, " LCRIT 1 1 "); GMT_Usage (API, -21, "Laplace distribution critical value for alpha = A");
GMT_Message (API, GMT_TIME_NONE, " LE 2 1 "); GMT_Usage (API, -21, "1 if A <= B, else 0");
GMT_Message (API, GMT_TIME_NONE, " LMSSCL 1 1 "); GMT_Usage (API, -21, "LMS scale estimate (LMS STD) of A");
GMT_Message (API, GMT_TIME_NONE, " LMSSCLW 1 1 "); GMT_Usage (API, -21, "Weighted LMS scale estimate (LMS STD) of A for weights in B");
GMT_Message (API, GMT_TIME_NONE, " LOG 1 1 "); GMT_Usage (API, -21, "log (A) (natural log)");
GMT_Message (API, GMT_TIME_NONE, " LOG10 1 1 "); GMT_Usage (API, -21, "log10 (A) (base 10)");
GMT_Message (API, GMT_TIME_NONE, " LOG1P 1 1 "); GMT_Usage (API, -21, "log (1+A) (accurate for small A)");
GMT_Message (API, GMT_TIME_NONE, " LOG2 1 1 "); GMT_Usage (API, -21, "log2 (A) (base 2)");
GMT_Message (API, GMT_TIME_NONE, " LOWER 1 1 "); GMT_Usage (API, -21, "The lowest (minimum) value of A");
GMT_Message (API, GMT_TIME_NONE, " LPDF 1 1 "); GMT_Usage (API, -21, "Laplace probability density function for z = A");
GMT_Message (API, GMT_TIME_NONE, " LRAND 2 1 "); GMT_Usage (API, -21, "Laplace random noise with mean A and std. deviation B");
GMT_Message (API, GMT_TIME_NONE, " LSQFIT 1 0 "); GMT_Usage (API, -21, "Current stack is [A | b]; Solve A * x = b via Cholesky decomposition");
GMT_Message (API, GMT_TIME_NONE, " LT 2 1 "); GMT_Usage (API, -21, "1 if A < B, else 0");
GMT_Message (API, GMT_TIME_NONE, " MAD 1 1 "); GMT_Usage (API, -21, "Median Absolute Deviation (L1 STD) of A");
GMT_Message (API, GMT_TIME_NONE, " MADW 2 1 "); GMT_Usage (API, -21, "Weighted Median Absolute Deviation (L1 STD) of A for weights in B");
GMT_Message (API, GMT_TIME_NONE, " MAX 2 1 "); GMT_Usage (API, -21, "Maximum of A and B");
GMT_Message (API, GMT_TIME_NONE, " MEAN 1 1 "); GMT_Usage (API, -21, "Mean value of A");
GMT_Message (API, GMT_TIME_NONE, " MEANW 2 1 "); GMT_Usage (API, -21, "Weighted mean value of A for weights in B");
GMT_Message (API, GMT_TIME_NONE, " MEDIAN 1 1 "); GMT_Usage (API, -21, "Median value of A");
GMT_Message (API, GMT_TIME_NONE, " MEDIANW 2 1 "); GMT_Usage (API, -21, "Weighted median value of A for weights in B");
GMT_Message (API, GMT_TIME_NONE, " MIN 2 1 "); GMT_Usage (API, -21, "Minimum of A and B");
GMT_Message (API, GMT_TIME_NONE, " MOD 2 1 "); GMT_Usage (API, -21, "A mod B (remainder after floored division)");
GMT_Message (API, GMT_TIME_NONE, " MODE 1 1 "); GMT_Usage (API, -21, "Mode value (Least Median of Squares) of A");
GMT_Message (API, GMT_TIME_NONE, " MODEW 2 1 "); GMT_Usage (API, -21, "Weighted mode value of A for weights in B");
GMT_Message (API, GMT_TIME_NONE, " MUL 2 1 "); GMT_Usage (API, -21, "A * B");
GMT_Message (API, GMT_TIME_NONE, " NAN 2 1 "); GMT_Usage (API, -21, "NaN if A == B, else A");
GMT_Message (API, GMT_TIME_NONE, " NEG 1 1 "); GMT_Usage (API, -21, "-A");
GMT_Message (API, GMT_TIME_NONE, " NEQ 2 1 "); GMT_Usage (API, -21, "1 if A != B, else 0");
GMT_Message (API, GMT_TIME_NONE, " NORM 1 1 "); GMT_Usage (API, -21, "Normalize (A) so min(A) = 0 and max(A) = 1");
GMT_Message (API, GMT_TIME_NONE, " NOT 1 1 "); GMT_Usage (API, -21, "NaN if A == NaN, 1 if A == 0, else 0");
GMT_Message (API, GMT_TIME_NONE, " NRAND 2 1 "); GMT_Usage (API, -21, "Normal, random values with mean A and std. deviation B");
GMT_Message (API, GMT_TIME_NONE, " OR 2 1 "); GMT_Usage (API, -21, "NaN if B == NaN, else A");
GMT_Message (API, GMT_TIME_NONE, " PERM 2 1 "); GMT_Usage (API, -21, "Permutations n_P_r, with n = A and r = B");
GMT_Message (API, GMT_TIME_NONE, " PLM 3 1 "); GMT_Usage (API, -21, "Associated Legendre polynomial P(A) degree B order C");
GMT_Message (API, GMT_TIME_NONE, " PLMg 3 1 "); GMT_Usage (API, -21, "Normalized associated Legendre polynomial P(A) degree B order C (geophysical convention)");
GMT_Message (API, GMT_TIME_NONE, " POP 1 0 "); GMT_Usage (API, -21, "Delete top element from the stack");
GMT_Message (API, GMT_TIME_NONE, " POW 2 1 "); GMT_Usage (API, -21, "A ^ B");
GMT_Message (API, GMT_TIME_NONE, " PPDF 2 1 "); GMT_Usage (API, -21, "Poisson probability density function for x = A and lambda = B");
GMT_Message (API, GMT_TIME_NONE, " PQUANT 2 1 "); GMT_Usage (API, -21, "The B'th Quantile (0-100%) of A");
GMT_Message (API, GMT_TIME_NONE, " PQUANTW 3 1 "); GMT_Usage (API, -21, "The C'th Quantile (0-100%) of A for weights in B");
GMT_Message (API, GMT_TIME_NONE, " PSI 1 1 "); GMT_Usage (API, -21, "Psi (or Digamma) of A");
GMT_Message (API, GMT_TIME_NONE, " PV 3 1 "); GMT_Usage (API, -21, "Legendre function Pv(A) of degree v = real(B) + imag(C)");
GMT_Message (API, GMT_TIME_NONE, " QV 3 1 "); GMT_Usage (API, -21, "Legendre function Qv(A) of degree v = real(B) + imag(C)");
GMT_Message (API, GMT_TIME_NONE, " R2 2 1 "); GMT_Usage (API, -21, "R2 = A^2 + B^2");
GMT_Message (API, GMT_TIME_NONE, " R2D 1 1 "); GMT_Usage (API, -21, "Convert Radians to Degrees");
GMT_Message (API, GMT_TIME_NONE, " RAND 2 1 "); GMT_Usage (API, -21, "Uniform random values between A and B");
GMT_Message (API, GMT_TIME_NONE, " RCDF 1 1 "); GMT_Usage (API, -21, "Rayleigh cumulative distribution function for z = A");
GMT_Message (API, GMT_TIME_NONE, " RCRIT 1 1 "); GMT_Usage (API, -21, "Rayleigh distribution critical value for alpha = A");
GMT_Message (API, GMT_TIME_NONE, " RGB2HSV 3 3 "); GMT_Usage (API, -21, "Convert rgb to hsv, with r = A, g = B and b = C");
GMT_Message (API, GMT_TIME_NONE, " RGB2LAB 3 3 "); GMT_Usage (API, -21, "Convert rgb to lab, with r = A, g = B and b = C");
GMT_Message (API, GMT_TIME_NONE, " RGB2XYZ 3 3 "); GMT_Usage (API, -21, "Convert rgb to xyz, with r = A, g = B and b = C");
GMT_Message (API, GMT_TIME_NONE, " RPDF 1 1 "); GMT_Usage (API, -21, "Rayleigh probability density function for z = A");
GMT_Message (API, GMT_TIME_NONE, " RINT 1 1 "); GMT_Usage (API, -21, "rint (A) (round to integral value nearest to A)");
GMT_Message (API, GMT_TIME_NONE, " RMS 1 1 "); GMT_Usage (API, -21, "Root-mean-square of A");
GMT_Message (API, GMT_TIME_NONE, " RMSW 2 1 "); GMT_Usage (API, -21, "Weighted Root-mean-square of A for weights in B");
GMT_Message (API, GMT_TIME_NONE, " ROLL 2 0 "); GMT_Usage (API, -21, "Cyclically shifts the top A stack items by an amount B");
GMT_Message (API, GMT_TIME_NONE, " ROTT 2 1 "); GMT_Usage (API, -21, "Rotate A by the (constant) shift B in the t-direction");
GMT_Message (API, GMT_TIME_NONE, " SEC 1 1 "); GMT_Usage (API, -21, "sec (A) (A in radians)");
GMT_Message (API, GMT_TIME_NONE, " SECD 1 1 "); GMT_Usage (API, -21, "sec (A) (A in degrees)");
GMT_Message (API, GMT_TIME_NONE, " SECH 1 1 "); GMT_Usage (API, -21, "sech (A)");
GMT_Message (API, GMT_TIME_NONE, " SIGN 1 1 "); GMT_Usage (API, -21, "sign (+1 or -1) of A");
GMT_Message (API, GMT_TIME_NONE, " SIN 1 1 "); GMT_Usage (API, -21, "sin (A) (A in radians)");
GMT_Message (API, GMT_TIME_NONE, " SINC 1 1 "); GMT_Usage (API, -21, "sinc (A) (sin (pi*A)/(pi*A))");
GMT_Message (API, GMT_TIME_NONE, " SIND 1 1 "); GMT_Usage (API, -21, "sin (A) (A in degrees)");
GMT_Message (API, GMT_TIME_NONE, " SINH 1 1 "); GMT_Usage (API, -21, "sinh (A)");
GMT_Message (API, GMT_TIME_NONE, " SKEW 1 1 "); GMT_Usage (API, -21, "Skewness of A");
GMT_Message (API, GMT_TIME_NONE, " SORT 3 1 "); GMT_Usage (API, -21, "Sort all columns based on column A in direction of B (-1 descending |+1 ascending)");
GMT_Message (API, GMT_TIME_NONE, " SQR 1 1 "); GMT_Usage (API, -21, "A^2");
GMT_Message (API, GMT_TIME_NONE, " SQRT 1 1 "); GMT_Usage (API, -21, "sqrt (A)");
GMT_Message (API, GMT_TIME_NONE, " STD 1 1 "); GMT_Usage (API, -21, "Standard deviation of A");
GMT_Message (API, GMT_TIME_NONE, " STDW 2 1 "); GMT_Usage (API, -21, "Weighted standard deviation of A for weights in B");
GMT_Message (API, GMT_TIME_NONE, " STEP 1 1 "); GMT_Usage (API, -21, "Heaviside step function H(A)");
GMT_Message (API, GMT_TIME_NONE, " STEPT 1 1 "); GMT_Usage (API, -21, "Heaviside step function H(t-A)");
GMT_Message (API, GMT_TIME_NONE, " SUB 2 1 "); GMT_Usage (API, -21, "A - B");
GMT_Message (API, GMT_TIME_NONE, " SUM 1 1 "); GMT_Usage (API, -21, "Cumulative sum of A");
GMT_Message (API, GMT_TIME_NONE, " SVDFIT 1 0 "); GMT_Usage (API, -21, "Current stack is [A | b]; solve A * x = B via SVD decomposition (see -E)");
GMT_Message (API, GMT_TIME_NONE, " TAN 1 1 "); GMT_Usage (API, -21, "tan (A) (A in radians)");
GMT_Message (API, GMT_TIME_NONE, " TAND 1 1 "); GMT_Usage (API, -21, "tan (A) (A in degrees)");
GMT_Message (API, GMT_TIME_NONE, " TANH 1 1 "); GMT_Usage (API, -21, "tanh (A)");
GMT_Message (API, GMT_TIME_NONE, " TAPER 1 1 "); GMT_Usage (API, -21, "Unit weights cosine-tapered to zero within A of end margins");
GMT_Message (API, GMT_TIME_NONE, " TCDF 2 1 "); GMT_Usage (API, -21, "Student's t cumulative distribution function for t = A and nu = B");
GMT_Message (API, GMT_TIME_NONE, " TN 2 1 "); GMT_Usage (API, -21, "Chebyshev polynomial Tn(-1<A<+1) of degree B");
GMT_Message (API, GMT_TIME_NONE, " TPDF 2 1 "); GMT_Usage (API, -21, "Student's t probability density function for t = A and nu = B");
GMT_Message (API, GMT_TIME_NONE, " TCRIT 2 1 "); GMT_Usage (API, -21, "Student's t distribution critical value for alpha = A and nu = B");
GMT_Message (API, GMT_TIME_NONE, " UPPER 1 1 "); GMT_Usage (API, -21, "The highest (maximum) value of A");
GMT_Message (API, GMT_TIME_NONE, " VAR 1 1 "); GMT_Usage (API, -21, "Variance of A");
GMT_Message (API, GMT_TIME_NONE, " VARW 2 1 "); GMT_Usage (API, -21, "Weighted variance of A for weights in B");
GMT_Message (API, GMT_TIME_NONE, " VPDF 3 1 "); GMT_Usage (API, -21, "Von Mises probability density function for angles = A, mu = B, and kappa = C");
GMT_Message (API, GMT_TIME_NONE, " WCDF 3 1 "); GMT_Usage (API, -21, "Weibull cumulative distribution function for x = A, scale = B, and shape = C");
GMT_Message (API, GMT_TIME_NONE, " WCRIT 3 1 "); GMT_Usage (API, -21, "Weibull distribution critical value for alpha = A, scale = B, and shape = C");
GMT_Message (API, GMT_TIME_NONE, " WPDF 3 1 "); GMT_Usage (API, -21, "Weibull probability density function for x = A, scale = B and shape = C");
GMT_Message (API, GMT_TIME_NONE, " XOR 2 1 "); GMT_Usage (API, -21, "B if A == NaN, else A");
GMT_Message (API, GMT_TIME_NONE, " XYZ2HSV 3 3 "); GMT_Usage (API, -21, "Convert xyz to hsv, with x = A, y = B and z = C");
GMT_Message (API, GMT_TIME_NONE, " XYZ2LAB 3 3 "); GMT_Usage (API, -21, "Convert xyz to lab, with x = A, y = B and z = C");
GMT_Message (API, GMT_TIME_NONE, " XYZ2RGB 3 3 "); GMT_Usage (API, -21, "Convert xyz to rgb, with x = A, y = B and z = C");
GMT_Message (API, GMT_TIME_NONE, " Y0 1 1 "); GMT_Usage (API, -21, "Bessel function of A (2nd kind, order 0)");
GMT_Message (API, GMT_TIME_NONE, " Y1 1 1 "); GMT_Usage (API, -21, "Bessel function of A (2nd kind, order 1)");
GMT_Message (API, GMT_TIME_NONE, " YN 2 1 "); GMT_Usage (API, -21, "Bessel function of A (2nd kind, order B)");
GMT_Message (API, GMT_TIME_NONE, " ZCRIT 1 1 "); GMT_Usage (API, -21, "Normal distribution critical value for alpha = A");
GMT_Message (API, GMT_TIME_NONE, " ZCDF 1 1 "); GMT_Usage (API, -21, "Normal cumulative distribution function for z = A");
GMT_Message (API, GMT_TIME_NONE, " ZPDF 1 1 "); GMT_Usage (API, -21, "Normal probability density function for z = A");
GMT_Message (API, GMT_TIME_NONE, " ROOTS 2 1 "); GMT_Usage (API, -21, "Treat col A as f(t) = 0 and returns its roots");
GMT_Usage (API, -2, "\nThe special symbols are:\n");
GMT_Message (API, GMT_TIME_NONE, " PI = "); GMT_Usage (API, -26, "3.1415926...");
GMT_Message (API, GMT_TIME_NONE, " E = "); GMT_Usage (API, -26, "2.7182818...");
GMT_Message (API, GMT_TIME_NONE, " EULER = "); GMT_Usage (API, -26, "0.5772156...");
GMT_Message (API, GMT_TIME_NONE, " PHI (golden ratio) = "); GMT_Usage (API, -26, "1.6180339...");
GMT_Message (API, GMT_TIME_NONE, " F_EPS (single eps) = "); GMT_Usage (API, -26, "1.192092896e-07");
GMT_Message (API, GMT_TIME_NONE, " D_EPS (double eps) = "); GMT_Usage (API, -26, "2.2204460492503131e-16");
GMT_Message (API, GMT_TIME_NONE, " TMIN = "); GMT_Usage (API, -26, "min \"time\" value set via -T");
GMT_Message (API, GMT_TIME_NONE, " TMAX = "); GMT_Usage (API, -26, "max \"time\" value set via -T");
GMT_Message (API, GMT_TIME_NONE, " TRANGE = "); GMT_Usage (API, -26, "range of \"time\" values");
GMT_Message (API, GMT_TIME_NONE, " TINC = "); GMT_Usage (API, -26, "increment of \"time\" values");
GMT_Message (API, GMT_TIME_NONE, " N = "); GMT_Usage (API, -26, "number of records");
GMT_Message (API, GMT_TIME_NONE, " T = "); GMT_Usage (API, -26, "table with t-coordinates");
GMT_Message (API, GMT_TIME_NONE, " TNORM = "); GMT_Usage (API, -26, "table with normalized [-1 to +1] t-coordinates");
GMT_Message (API, GMT_TIME_NONE, " TROW = "); GMT_Usage (API, -26, "table with row numbers 0, 1, ..., N-1");
GMT_Usage (API, -2, "\nUse macros for frequently used long expressions; see the gmtmath man page. "
"Store stack to named variable via STO@<label>, recall via [RCL]@<label>, clear via CLR@<label>.");
GMT_Message (API, GMT_TIME_NONE, "\n OPTIONAL ARGUMENTS:\n");
GMT_Usage (API, 1, "\n-A<ftable>[+e][+r][+s|w]");
GMT_Usage (API, -2, "Set up and solve a linear system A x = b, and return vector x. "
"Requires -N and initializes extended matrix [A | b] from <ftable> holding t and f(t) only. "
"Vector t goes into column <t_col> while f(t) goes into column <n_col> - 1 (i.e., r.h.s. vector b). "
"No additional data files are read. Output will be a single column with coefficients.");
GMT_Usage (API, 3, "+e Evaluate solution and write t, f(t), the solution, residuals[, weight|sigma].");
GMT_Usage (API, 3, "+r Only place f(t) in b and leave A initialized to zeros.");
GMT_Usage (API, 3, "+s Third column contains 1-sigmas.");
GMT_Usage (API, 3, "+w Third column contains weights");
GMT_Usage (API, -2, "Use either LSQFIT or SVDFIT to solve the [weighted] linear system.");
GMT_Usage (API, 1, "\n-C<cols>");
GMT_Usage (API, -2, "Change which columns to operate on [Default is all except time]. "
"Plain -C reverts to the default, -Cr toggles current settings, and -Ca selects all columns.");
GMT_Usage (API, 1, "\n-E<eigen>");
GMT_Usage (API, -2, "Set minimum eigenvalue used by LSQFIT and SVDFIT [1e-7].");
GMT_Usage (API, 1, "\n-I Reverse the output sequence into descending order [ascending].");
GMT_Usage (API, 1, "\n-L Apply operators on a per-segment basis [accumulates operations across file].");
GMT_Usage (API, 1, "\n-N<n_col>[/<t_col>]");
GMT_Usage (API, -2, "Set the number of columns and optionally the id of the time column (0 is first) [2/0]. "
"If input files are given, -N will add extra columns initialized to zero, if needed.");
GMT_Usage (API, 1, "\n-Q[%s|n]", GMT_DIM_UNITS_DISPLAY);
GMT_Usage (API, -2, "Quick scalar calculator (Shorthand for -Ca -N1/0 -T0/0/1). "
"Allows constants to have dimensional units (i.e., %s); if so the answer is converted using PROJ_LENGTH_UNIT. "
"Optionally, append another unit or n for no unit conversion on output.", GMT_DIM_UNITS_DISPLAY);
GMT_Usage (API, 1, "\n-S[f|l]");
GMT_Usage (API, -2, "Only write first row upon completion of calculations [write all rows]. "
"Optionally, append l for last row or f for first row [Default].");
GMT_Usage (API, 1, "\n-T[<file>|<list>|<min>/<max>/<inc>[+b|i|l|n]]");
GMT_Usage (API, -2, "Set domain from <min> to <max> in steps of <inc>. Control setup via modifiers:");
GMT_Usage (API, 3, "+b Select log2 spacing in <inc>");
GMT_Usage (API, 3, "+i Indicate <inc> is the reciprocal of desired <inc> (e.g., 3 for 0.3333.....).");
GMT_Usage (API, 3, "+l Select log10 spacing via <inc> = 1,2,3.");
GMT_Usage (API, 3, "+n Let <inc> mean the number of points instead. of increment");
GMT_Usage (API, -2, "Alternatively, give a <file> with output times in the first column, or a comma-separated <list>. "
"If no domain is appended then we assume no time, i.e., only data columns are present. "
"This choice implicitly sets -Ca.");
GMT_Option (API, "V,bi,bo,d,e,f,g,h,i,o,q,s,w,.");
return (GMT_MODULE_USAGE);
}
static int parse (struct GMT_CTRL *GMT, struct GMTMATH_CTRL *Ctrl, struct GMT_OPTION *options) {
/* This parses the options provided to gmtmath and sets parameters in CTRL.
* Any GMT common options will override values set previously by other commands.
* It also replaces any file names specified as input or output with the data ID
* returned when registering these sources/destinations with the API.
*/
unsigned int n_errors = 0, k = 0;
bool missing_equal = true;
char *c = NULL, *t_arg = NULL;
struct GMT_OPTION *opt = NULL;
struct GMTAPI_CTRL *API = GMT->parent;
int gmt_parse_o_option (struct GMT_CTRL *GMT, char *arg);
for (opt = options; opt; opt = opt->next) {
switch (opt->option) {
case '<': /* Input files */
if (opt->arg[0] == '=' && opt->arg[1] == 0) { /* No it was an = [outfile] sequence */
missing_equal = false;
opt->option = GMT_OPT_OUTFILE; /* Prevents further use later */
if (opt->next && (opt->next->option == GMT_OPT_INFILE || opt->next->option == GMT_OPT_OUTFILE)) {
opt = opt->next; /* Skip ahead */
n_errors += gmt_M_repeated_module_option (API, Ctrl->Out.active);
if (opt->arg[0]) Ctrl->Out.file = strdup (opt->arg); /* Optional output file since stdout is default */
opt->option = GMT_OPT_OUTFILE; /* Prevents further use later */
}
}
break;
case '>': /* Output file specified via an API; set required output file here */
n_errors += gmt_M_repeated_module_option (API, Ctrl->Out.active);
n_errors += gmt_get_required_file (GMT, opt->arg, opt->option, 0, GMT_IS_DATASET, GMT_OUT, GMT_FILE_LOCAL, &(Ctrl->Out.file));
missing_equal = false;
break;
case '#': /* Skip numbers */
break;
/* Processes program-specific parameters */
case 'A': /* y(x) table for LSQFIT/SVDFIT operations */
n_errors += gmt_M_repeated_module_option (API, Ctrl->A.active);
if (opt->arg[0] == '-') { /* Old-style leading hyphen to the filename has been replaced by modifier +r */
if (gmt_M_compat_check (GMT, 5)) {
GMT_Report (API, GMT_MSG_COMPAT, "The leading hyphen in -A is deprecated. Append modifier +r instead.\n");
Ctrl->A.null = true;
k = 1;
}
else {
GMT_Report (API, GMT_MSG_ERROR, "Option -A: Unable to decode arguments\n");
n_errors++;
}
}
if ((c = strchr(opt->arg, '+')) != NULL && strchr("ersw", c[1]) != NULL) { /* Got a valid modifier */
unsigned int pos = 0;
char p[GMT_LEN256] = {""};
c[0] = '\0'; /* Temporarily chop off modifiers */
Ctrl->A.file = strdup (&opt->arg[k]);
c[0] = '+'; /* Restore the modifier */
while (gmt_strtok (c, "+", &pos, p)) {
switch (p[0]) {
case 'e': Ctrl->A.e_mode = GMTMATH_EVALUATE; break; /* Evaluate solution */
case 'r': Ctrl->A.null = true; break; /* Only set rhs of equation */
case 's': Ctrl->A.w_mode = GMTMATH_SIGMAS; break; /* Got t,y,s */
case 'w': Ctrl->A.w_mode = GMTMATH_WEIGHTS; break; /* Got t,y,w */
default: n_errors++; break;
}
}
}
else /* No modifiers, selected default output of coefficient column */
Ctrl->A.file = strdup (&opt->arg[k]);
break;
case 'C': /* Processed in the main loop but not here; just skip */
break;
case 'E': /* Set minimum eigenvalue cutoff */
n_errors += gmt_M_repeated_module_option (API, Ctrl->E.active);
n_errors += gmt_get_required_double (GMT, opt->arg, opt->option, 0, &Ctrl->E.eigen);
break;
case 'F': /* Now obsolete, using -o instead */
if (gmt_M_compat_check (GMT, 4)) {
GMT_Report (API, GMT_MSG_COMPAT, "Option -F is deprecated; use -o instead\n");
gmt_parse_o_option (GMT, opt->arg);
}
else
n_errors += gmt_default_option_error (GMT, opt);
break;
case 'I': /* Reverse output order */
n_errors += gmt_M_repeated_module_option (API, Ctrl->I.active);
n_errors += gmt_get_no_argument (GMT, opt->arg, opt->option, 0);
break;
case 'L': /* Apply operator per segment basis */
n_errors += gmt_M_repeated_module_option (API, Ctrl->L.active);
n_errors += gmt_get_no_argument (GMT, opt->arg, opt->option, 0);
break;
case 'N': /* Sets no of columns and optionally the time column [0] */
n_errors += gmt_M_repeated_module_option (API, Ctrl->N.active);
if (sscanf (opt->arg, "%" PRIu64 "/%" PRIu64, &Ctrl->N.ncol, &Ctrl->N.tcol) == 1) Ctrl->N.tcol = 0;
break;
case 'Q': /* Quick for -Ca -N1/0 -T0/0/1 */
n_errors += gmt_M_repeated_module_option (API, Ctrl->Q.active);
if (opt->arg[0] == 'n') /* Want no unit conversion on output */
Ctrl->Q.unit = GMT_INCH; /* We do this which will convert from inch to inch, i.e., no change */
else if (opt->arg[0] && strchr (GMT_DIM_UNITS, opt->arg[0])) /* Want a specific unit conversion from inch to this unit on output */
Ctrl->Q.unit = gmt_get_dim_unit (GMT, opt->arg[0]);
/* else: Default GMT unit on output */
break;
case 'S': /* Only want one row (first or last) */
n_errors += gmt_M_repeated_module_option (API, Ctrl->S.active);
switch (opt->arg[0]) {
case 'f': case 'F': case '\0':
Ctrl->S.mode = -1; break;
case 'l': case 'L':
Ctrl->S.mode = +1; break;
default:
GMT_Report (API, GMT_MSG_ERROR, "Option -S: Syntax is -S[f|l]\n");
n_errors++;
break;
}
break;
case 'T': /* Either get a file with time coordinate or a min/max/dt setting */
n_errors += gmt_M_repeated_module_option (API, Ctrl->T.active);
t_arg = opt->arg;
break;
default: /* Report bad options */
n_errors += gmt_default_option_error (GMT, opt);
}
}
if (Ctrl->T.active) { /* Do this one here since we need Ctrl->N.col to be set first, if selected */
if (t_arg && t_arg[0]) { /* Gave an argument, so we can parse and create the array */
n_errors += gmt_parse_array (GMT, 'T', t_arg, &(Ctrl->T.T), GMT_ARRAY_TIME | GMT_ARRAY_SCALAR |
GMT_ARRAY_RANGE, (unsigned int)Ctrl->N.tcol);
n_errors += gmt_create_array (GMT, 'T', &(Ctrl->T.T), NULL, NULL);
}
else
Ctrl->T.notime = true;
}
n_errors += gmt_M_check_condition (GMT, Ctrl->A.active && gmt_access (GMT, Ctrl->A.file, R_OK),
"Option -A: Cannot read file %s!\n", Ctrl->A.file);
n_errors += gmt_M_check_condition (GMT, missing_equal, "Usage is <operations> = [outfile]\n");
n_errors += gmt_M_check_condition (GMT, Ctrl->Q.active && (Ctrl->T.active || Ctrl->N.active || Ctrl->C.active),
"Cannot use -T, -N, or -C when -Q has been set\n");
n_errors += gmt_M_check_condition (GMT, Ctrl->N.active && (Ctrl->N.ncol == 0 || Ctrl->N.tcol >= Ctrl->N.ncol),
"Option -N must have positive n_cols and 0 <= t_col < n_col\n");
return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR);
}
GMT_LOCAL unsigned int gmtmath_assign_ptrs (struct GMT_CTRL *GMT, unsigned int last, struct GMTMATH_STACK *S[], struct GMT_DATATABLE **T, struct GMT_DATATABLE **T_prev) { /* Centralize the assignment of previous stack ID and the current and previous stack tables */
unsigned int prev;
if (last == 0) { /* User error in requesting more items that presently on the stack */
GMT_Report (GMT->parent, GMT_MSG_ERROR, "Not enough items on the stack\n");
return UINT_MAX; /* Error flag */
}
prev = last - 1;
*T = (S[last]->constant && !S[last]->D) ? NULL : S[last]->D->table[0];
*T_prev = S[prev]->D->table[0];
return prev;
}
/* -----------------------------------------------------------------
* Definitions of all operator functions
* -----------------------------------------------------------------*/
/* Note: The OPERATOR: **** lines were used to extract syntax for documentation in the old days.
* the first number is input args and the second is the output args. */
GMT_LOCAL int gmtmath_ABS (struct GMT_CTRL *GMT, struct GMTMATH_INFO *info, struct GMTMATH_STACK *S[], unsigned int last, unsigned int col) {
/*OPERATOR: ABS 1 1 abs (A). */
uint64_t s, row;
double a = 0.0;
struct GMT_DATATABLE *T = S[last]->D->table[0];
gmt_M_unused (GMT);
if (S[last]->constant) a = fabs (S[last]->factor);
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] = (S[last]->constant) ? a : fabs (T->segment[s]->data[col][row]);
return 0;
}
GMT_LOCAL int gmtmath_ACOS (struct GMT_CTRL *GMT, struct GMTMATH_INFO *info, struct GMTMATH_STACK *S[], unsigned int last, unsigned int col) {
/*OPERATOR: ACOS 1 1 acos (A). */
uint64_t s, row;
double a = 0.0;
struct GMT_DATATABLE *T = S[last]->D->table[0];
if (S[last]->constant && fabs (S[last]->factor) > 1.0) GMT_Report (GMT->parent, GMT_MSG_WARNING, "|Operand| > 1 for ACOS!\n");
if (S[last]->constant) a = d_acos (S[last]->factor);
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] = (S[last]->constant) ? a : d_acos (T->segment[s]->data[col][row]);
return 0;
}
GMT_LOCAL int gmtmath_ACOSD (struct GMT_CTRL *GMT, struct GMTMATH_INFO *info, struct GMTMATH_STACK *S[], unsigned int last, unsigned int col) {
/*OPERATOR: ACOSD 1 1 acosd (A). */
uint64_t s, row;
double a = 0.0;
struct GMT_DATATABLE *T = S[last]->D->table[0];
if (S[last]->constant && fabs (S[last]->factor) > 1.0) GMT_Report (GMT->parent, GMT_MSG_WARNING, "|Operand| > 1 for ACOSD!\n");
if (S[last]->constant) a = R2D * d_acos (S[last]->factor);
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] = (S[last]->constant) ? a : R2D * d_acos (T->segment[s]->data[col][row]);
return 0;
}
GMT_LOCAL int gmtmath_ACOSH (struct GMT_CTRL *GMT, struct GMTMATH_INFO *info, struct GMTMATH_STACK *S[], unsigned int last, unsigned int col)
/*OPERATOR: ACOSH 1 1 acosh (A). */