Mercurial > hg > medcouple

annotate medcouple.c++ @ 74:305b7361a5bd default tip @

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
showalgo: save a snapshot instead of waiting for keyboard input
author Jordi Gutiérrez Hermoso <jordigh@octave.org>
date Sun, 29 May 2016 19:05:01 -0400
parents 9b6502f3f3bc
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
1 /*
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
2 Algorithm for the skewness estimator medcouple (MC)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
3 --------------------------------------------------
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
4 ( originally matlabmc.c and also mc/mcrsoft/spmc.c )
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
5 */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
6
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
7 #include <stdio.h>
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
8 #include <stdlib.h>
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
9 #include <math.h>
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
10 #include <algorithm>
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
11
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
12 #include <inttypes.h>
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
13 // -> int64_t
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
14
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
15 template <typename T> int sign(T val) {
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
16 return (T(0) < val) - (val < T(0));
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
17 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
18
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
19 /* Interface routines to be called via .C() and those from API : */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
20 #include "wgt_himed.c"
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
21 /*
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
22 including
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
23
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
24 whimed_i(a,iw,n): the weighted high median of an array a of length n,
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
25 using the positive integer weights iw[].
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
26
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
27 * which is in ./wgt_himed.c_templ
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
28 * ~~~~~~~~~~~~~~~~~
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
29 */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
30
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
31 /* Includes the auxiliary function
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
32
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
33 h_kern(a,b, ai,bi,ab, eps): the values h(a,b) needed to compute the mc
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
34 */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
35 static
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
36 double h_kern(double a, double b, int ai, int bi, int ab, double eps);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
37
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
38 /* MM: The tolerance 'eps1' and 'eps2' can now be passed from R;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
39 * the original code had only one 'eps' for both and hardcoded
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
40 * eps = 0.0000000000001; /.* == 1e-13 )) *./
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
41 */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
42 /* MK: eps1: for (relative) "equality" checks
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
43 * eps2: used to check for over- and underflow, respectively
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
44 * therefore I suggest eps1 = DBL_EPS and eps2 = DBL_MIN
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
45 */
2
5fa249d59d58 Allow to set the trace level
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 0
diff changeset
46 double mc_C(double *z, int n, double *eps, int trace_lev)
0
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
47 {
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
48 /* NOTE:
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
49 eps = c(eps1, eps2)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
50 iter := c(maxit, trace.lev) as input
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
51 = c(it, converged) as output
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
52 */
2
5fa249d59d58 Allow to set the trace level
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 0
diff changeset
53 int it = 0;
0
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
54 bool converged = true;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
55 double medc; // "the" result
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
56
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
57 if (n < 3) {
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
58 return 0.;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
59 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
60 /* copy data before sort()ing in place, also reflecting it */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
61 /* NOTE: x[0] is empty --- so 1-indexing is used (MM: danger!) */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
62 double *x = (double *) malloc((n+1)*sizeof(double));
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
63 x[0] = 0;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
64 for (int i = 0; i < n; i++)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
65 x[i+1] = -z[i];
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
66
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
67 std::sort(&x[1], &x[1] + n); /* full sort */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
68
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
69 double xmed; // := median( x[1:n] ) = - median( z[0:(n-1)] ):
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
70 if (n%2) { /* n even */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
71 xmed = x[(n/2)+1];
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
72 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
73 else { /* n odd */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
74 int ind = (n/2);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
75 xmed = (x[ind] + x[ind+1])/2;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
76 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
77
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
78 if (fabs(x[1] - xmed) < eps[0] * (eps[0] + fabs(xmed))) {
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
79 return -1;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
80 } else if (fabs(x[n] - xmed) < eps[0] * (eps[0] + fabs(xmed))) {
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
81 return 1;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
82 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
83 /* else : median is not at the border ------------------- */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
84
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
85 if(trace_lev)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
86 printf("mc_C(z[1:%d], trace_lev=%d): Median = %g (not at the border)\n",
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
87 n, trace_lev, -xmed);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
88
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
89 int i,j;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
90 /* center x[] wrt median --> such that then median( x[1:n] ) == 0 */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
91 for (i = 1; i <= n; i++)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
92 x[i] -= xmed;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
93
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
94 /* Now scale to inside [-0.5, 0.5] and flip sign such that afterwards
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
95 * x[1] >= x[2] >= ... >= x[n] */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
96 double xden = -2 * std::max(-x[1], x[n]);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
97 for (i = 1; i <= n; i++)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
98 x[i] /= xden;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
99 xmed /= xden;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
100 if(trace_lev >= 2)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
101 printf(" x[] has been rescaled (* 1/s) with s = %g\n", -xden);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
102
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
103 j = 1;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
104 double x_eps = eps[0] * (eps[0] + fabs(xmed));
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
105 while (x[j] > x_eps && j <= n) { /* test relative to xmed */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
106 /* x1[j] = x[j]; */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
107 j++;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
108 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
109 if(trace_lev >= 2)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
110 printf(" x1[] := {x | x_j > x_eps = %g} has %d (='j-1') entries\n",
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
111 x_eps, j-1);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
112 i = 1;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
113 double *x2 = x+j-1; /* pointer -- corresponding to x2[i] = x[j]; */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
114 while (x[j] > -x_eps && j <= n) { /* test relative to xmed */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
115 /* x1[j] = x[j]; */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
116 /* x2[i] = x[j]; */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
117 j++;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
118 i++;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
119 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
120 /* now x1[] := {x | x_j > -eps} also includes the median (0) */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
121 if(trace_lev >= 2)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
122 printf("'median-x' {x | -eps < x_i <= eps} has %d (= 'k') entries\n",
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
123 i-1);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
124 int h1 = j-1, /* == size of x1[] == the sum of those two sizes above */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
125 /* conceptually, x2[] := {x | x_j <= eps} (which includes the median 0) */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
126 h2 = i + (n-j);// == size of x2[] == maximal size of whimed() arrays
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
127
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
128 if(trace_lev)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
129 printf(" now allocating 2+5 work arrays of size (1+) h2=%d each:\n", h2);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
130 /* work arrays for whimed_i() : allocate *once* only !! */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
131 double *acand = (double *) malloc(h2*sizeof(double)),
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
132 *a_srt = (double *) malloc(h2*sizeof(double));
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
133
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
134 int *iw_cand= (int *) malloc(h2*sizeof(int)),
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
135 /* work arrays for the fast-median-of-table algorithm:
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
136 * currently still with 1-indexing */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
137 *left = (int *) malloc((h2+1)*sizeof(int)),
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
138 *right = (int *) malloc((h2+1)*sizeof(int)),
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
139 *p = (int *) malloc((h2+1)*sizeof(int)),
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
140 *q = (int *) malloc((h2+1)*sizeof(int));
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
141
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
142 for (i = 1; i <= h2; i++) {
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
143 left [i] = 1;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
144 right[i] = h1;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
145 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
146 int64_t nr = ((int64_t) h1) * ((int64_t) h2), // <-- careful to *NOT* overflow
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
147 knew = nr/2 +1;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
148 if(trace_lev >= 2)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
149 printf(" (h1,h2, nr, knew) = (%d,%d, %.0f, %.0f)\n",
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
150 h1,h2, (double)nr, (double)knew);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
151
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
152 double trial = -2./* -Wall */;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
153 double *work= (double *) malloc(n*sizeof(double));
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
154 int *iwt = (int *) malloc(n*sizeof(int));
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
155 bool IsFound = false;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
156 int nl = 0,
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
157 neq = 0;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
158 /* MK: 'neq' counts the number of observations in the
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
159 * inside the tolerance range, i.e., where left > right + 1,
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
160 * since we would miss those when just using 'nl-nr'.
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
161 * This is to prevent index overflow in work[] later on.
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
162 * left might be larger than right + 1 since we are only
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
163 * testing with accuracy eps_trial and therefore there might
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
164 * be more than one observation in the `tolerance range`
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
165 * between < and <=.
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
166 */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
167 while (!IsFound && (nr-nl+neq > n) && it < 1e5)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
168 {
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
169 int64_t sum_p, sum_q;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
170 it++;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
171 j = 0;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
172 for (i = 1; i <= h2; i++)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
173 if (left[i] <= right[i]) {
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
174 iwt[j] = right[i] - left[i]+1;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
175 int k = left[i] + (iwt[j]/2);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
176 work[j] = h_kern(x[k], x2[i], k, i, h1+1, eps[1]);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
177 j++;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
178 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
179 if(trace_lev >= 4) {
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
180 printf(" before whimed(): work and iwt, each [0:(%d-1)]:\n", j);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
181 if(j >= 100) {
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
182 for(i=0; i < 90; i++) printf(" %8g", work[i]); printf("\n ... ");
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
183 for(i=j-4; i < j; i++)printf(" %8g", work[i]); printf("\n");
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
184 for(i=0; i < 90; i++) printf(" %8d", iwt [i]); printf("\n ... ");
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
185 for(i=j-4; i < j; i++)printf(" %8d", iwt [i]); printf("\n");
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
186 } else { // j <= 99
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
187 for(i=0; i < j; i++) printf(" %8g", work[i]); printf("\n");
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
188 for(i=0; i < j; i++) printf(" %8d", iwt [i]); printf("\n");
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
189 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
190 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
191 trial = whimed_i(work, iwt, j, acand, a_srt, iw_cand);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
192 double eps_trial = eps[0] * (eps[0] + fabs(trial));
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
193 if(trace_lev >= 3)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
194 printf("%2s it=%2d, whimed(*, n=%6d)= %8g ", " ", it, j, trial);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
195
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
196 j = 1;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
197 for (i = h2; i >= 1; i--) {
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
198 while (j <= h1 && h_kern(x[j],x2[i],j,i,h1+1,eps[1]) - trial > eps_trial) {
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
199 // while (j <= h1 && h_kern(x[j],x2[i],j,i,h1+1,eps[1]) > trial) {
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
200 if (trace_lev >= 5)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
201 printf("\nj=%3d, i=%3d, x[j]=%g, x2[i]=%g, h=%g",
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
202 j, i, x[j], x2[i],
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
203 h_kern(x[j],x2[i],j,i,h1+1,eps[1]));
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
204 j++;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
205 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
206 /* for(; j <= h1; j++) { */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
207 /* register double h = h_kern(x[j],x2[i],j,i,h1+1,eps[1]); */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
208 /* if(h > trial) break; */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
209 /* } */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
210 p[i] = j-1;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
211 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
212 j = h1;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
213 for (i = 1, sum_p=0, sum_q=0; i <= h2; i++) {
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
214 while (j >= 1 && trial - h_kern(x[j],x2[i],j,i,h1+1,eps[1]) > eps_trial)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
215 // while (j >= 1 && h_kern(x[j],x2[i],j,i,h1+1,eps[1]) < trial)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
216 j--;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
217 q[i] = j+1;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
218
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
219 sum_p += p[i];
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
220 sum_q += j;/* = q[i]-1 */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
221 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
222
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
223 if(trace_lev >= 3) {
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
224 if (trace_lev == 3)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
225 printf("sum_(p,q)= (%.0f,%.0f)", (double)sum_p, (double)sum_q);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
226 else { /* trace_lev >= 4 */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
227 printf("\n%3s p[1:%d]:", "", h2);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
228 bool lrg = h2 >= 100;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
229 int i_m = lrg ? 95 : h2;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
230 for(i = 1; i <= i_m; i++) printf(" %2d", p[i]); if(lrg) printf(" ...");
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
231 printf(" sum=%4.0f\n%3s q[1:%d]:", (double)sum_p, "", h2);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
232 for(i = 1; i <= i_m; i++) printf(" %2d", q[i]); if(lrg) printf(" ...");
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
233 printf(" sum=%4.0f\n", (double)sum_q);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
234 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
235 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
236
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
237 if (knew <= sum_p) {
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
238 if(trace_lev >= 3)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
239 printf("; sum_p >= kn\n");
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
240 for (i = 1, neq = 0; i <= h2; i++) {
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
241 right[i] = p[i];
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
242 if (left[i] > right[i]+1) neq += left[i]-right[i]-1;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
243 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
244 nr = sum_p;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
245 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
246 else { /* knew > sum_p */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
247 IsFound = (knew <= sum_q); /* i.e. sum_p < knew <= sum_q */;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
248
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
249 if(trace_lev >= 3)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
250 printf("; s_p < kn ?<=? s_q: %s\n", IsFound ? "TRUE": "no");
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
251 if(IsFound) {
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
252 medc = trial;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
253 } else { /* knew > sum_q */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
254 for (i = 1; i <= h2; i++) {
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
255 left[i] = q[i];
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
256 if (left[i] > right[i]+1) neq += left[i]-right[i]-1;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
257 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
258 nl = sum_q;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
259 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
260 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
261
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
262 } /* end while loop */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
263
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
264 converged = IsFound || (nr-nl+neq <= n);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
265 if(!converged) {
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
266 /* still: */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
267 medc = trial;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
268 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
269
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
270 if (converged && !IsFound) { /* e.g., for mc(1:4) : */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
271 j = 0;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
272 for (i = 1; i <= h2; i++) {
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
273 if (left[i] <= right[i]) {
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
274 for (int k = left[i]; k <= right[i]; k++) {
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
275 work[j] = -h_kern(x[k],x2[i],k,i,h1+1,eps[1]);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
276 j++;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
277 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
278 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
279 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
280 if(trace_lev)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
281 printf(" not found [it=%d, (nr,nl) = (%d,%d)],"
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
282 " -> (knew-nl, j) = (%d,%d)\n",
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
283 it, nr, nl, knew-nl, j);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
284 /* using rPsort(work, n,k), since we don't need work[] anymore:*/
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
285 std::nth_element(work, work + knew-nl-1, work + j);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
286 medc = - work[knew-nl-1];
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
287 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
288
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
289 if(converged && trace_lev >= 2)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
290 printf("converged in %d iterations\n", it);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
291
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
292 return medc;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
293
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
294 } /* end{ mc_C } */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
295
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
296
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
297 static
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
298 double h_kern(double a, double b, int ai, int bi, int ab, double eps)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
299 {
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
300 /* if (fabs(a-b) <= DBL_MIN) */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
301 /* check for zero division and positive b */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
302 if (fabs(a-b) < 2.0*eps || b > 0)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
303 return sign((double)(ab - (ai+bi)));
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
304
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
305 /* else */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
306 return (a+b)/(a-b);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
307 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
308
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
309 #include <vector>
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
310 #include <iostream>
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
311 #include <string>
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
312 #include <fstream>
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
313 #include <iomanip>
2
5fa249d59d58 Allow to set the trace level
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 0
diff changeset
314 #include <sstream>
0
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
315
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
316 int main(int argc, char** argv) {
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
317 using namespace std;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
318
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
319 vector<double> eps = {numeric_limits<double>::epsilon(),
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
320 numeric_limits<double>::min()};
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
321
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
322 string fname;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
323 if (argc > 1) {
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
324 fname = argv[1];
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
325 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
326 else {
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
327 cerr << "Please provide input file." << endl;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
328 return 1;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
329 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
330
2
5fa249d59d58 Allow to set the trace level
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 0
diff changeset
331 double trace_lev = 0;
5fa249d59d58 Allow to set the trace level
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 0
diff changeset
332 if(argc > 2) {
5fa249d59d58 Allow to set the trace level
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 0
diff changeset
333 stringstream ss;
5fa249d59d58 Allow to set the trace level
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 0
diff changeset
334 ss << argv[2];
5fa249d59d58 Allow to set the trace level
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 0
diff changeset
335 ss >> trace_lev;
5fa249d59d58 Allow to set the trace level
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 0
diff changeset
336 }
5fa249d59d58 Allow to set the trace level
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 0
diff changeset
337
0
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
338 ifstream ifs(fname);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
339 if (not ifs) {
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
340 cerr << "File " << fname << " not found." << endl;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
341 return 1;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
342 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
343
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
344 vector<double> data;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
345 double datum;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
346 while (ifs >> datum) {
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
347 data.push_back(datum);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
348 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
349
2
5fa249d59d58 Allow to set the trace level
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 0
diff changeset
350 cout << setprecision(16) << mc_C(&data[0], data.size(), &eps[0], trace_lev) << endl;
0
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
351
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
352 return 0;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
353 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
354