1/* $Id$Revision: */
2/* vim:set shiftwidth=4 ts=8: */
3
4/*************************************************************************
5 * Copyright (c) 2011 AT&T Intellectual Property
6 * All rights reserved. This program and the accompanying materials
7 * are made available under the terms of the Eclipse Public License v1.0
8 * which accompanies this distribution, and is available at
9 * http://www.eclipse.org/legal/epl-v10.html
10 *
11 * Contributors: See CVS logs. Details at http://www.graphviz.org/
12 *************************************************************************/
13
14#include "general.h"
15#include <errno.h>
16
17#ifdef DEBUG
18double _statistics[10];
19#endif
20
21real vector_median(int n, real *x){
22 /* find the median value in a list of real */
23 int *p = NULL;
24 real res;
25 vector_ordering(n, x, &p, TRUE);
26
27 if ((n/2)*2 == n){
28 res = 0.5*(x[p[n/2-1]] + x[p[n/2]]);
29 } else {
30 res = x[p[n/2]];
31 }
32 FREE(p);
33 return res;
34}
35real vector_percentile(int n, real *x, real y){
36 /* find the value such that y% of element of vector x is <= that value.
37 y: a value between 0 and 1.
38 */
39 int *p = NULL, i;
40 real res;
41 vector_ordering(n, x, &p, TRUE);
42
43
44 y = MIN(y, 1);
45 y = MAX(0, y);
46
47 i = n*y;
48 res = x[p[i]];
49 FREE(p); return res;
50}
51
52real drand(){
53 return rand()/(real) RAND_MAX;
54}
55
56int irand(int n){
57 /* 0, 1, ..., n-1 */
58 assert(n > 1);
59 /*return (int) MIN(floor(drand()*n),n-1);*/
60 return rand()%n;
61}
62
63int *random_permutation(int n){
64 int *p;
65 int i, j, pp, len;
66 if (n <= 0) return NULL;
67 p = MALLOC(sizeof(int)*n);
68 for (i = 0; i < n; i++) p[i] = i;
69
70 len = n;
71 while (len > 1){
72 j = irand(len);
73 pp = p[len-1];
74 p[len-1] = p[j];
75 p[j] = pp;
76 len--;
77 }
78 return p;
79}
80
81
82real* vector_subtract_from(int n, real *x, real *y){
83 /* y = x-y */
84 int i;
85 for (i = 0; i < n; i++) y[i] = y[i] - x[i];
86 return y;
87}
88real* vector_subtract_to(int n, real *x, real *y){
89 /* y = x-y */
90 int i;
91 for (i = 0; i < n; i++) y[i] = x[i] - y[i];
92 return y;
93}
94real* vector_add_to(int n, real *x, real *y){
95 /* y = x-y */
96 int i;
97 for (i = 0; i < n; i++) y[i] = x[i] + y[i];
98 return y;
99}
100
101real vector_product(int n, real *x, real *y){
102 real res = 0;
103 int i;
104 for (i = 0; i < n; i++) res += x[i]*y[i];
105 return res;
106}
107
108real* vector_saxpy(int n, real *x, real *y, real beta){
109 /* y = x+beta*y */
110 int i;
111 for (i = 0; i < n; i++) y[i] = x[i] + beta*y[i];
112 return y;
113}
114
115real* vector_saxpy2(int n, real *x, real *y, real beta){
116 /* x = x+beta*y */
117 int i;
118 for (i = 0; i < n; i++) x[i] = x[i] + beta*y[i];
119 return x;
120}
121
122void vector_print(char *s, int n, real *x){
123 int i;
124 printf("%s{",s);
125 for (i = 0; i < n; i++) {
126 if (i > 0) printf(",");
127 printf("%f",x[i]);
128 }
129 printf("}\n");
130}
131
132void vector_take(int n, real *v, int m, int *p, real **u){
133 /* take m elements v[p[i]]],i=1,...,m and oput in u */
134 int i;
135
136 if (!*u) *u = MALLOC(sizeof(real)*m);
137
138 for (i = 0; i < m; i++) {
139 assert(p[i] < n && p[i] >= 0);
140 (*u)[i] = v[p[i]];
141 }
142
143}
144
145void vector_float_take(int n, float *v, int m, int *p, float **u){
146 /* take m elements v[p[i]]],i=1,...,m and oput in u */
147 int i;
148
149 if (!*u) *u = MALLOC(sizeof(float)*m);
150
151 for (i = 0; i < m; i++) {
152 assert(p[i] < n && p[i] >= 0);
153 (*u)[i] = v[p[i]];
154 }
155
156}
157
158int comp_ascend(const void *s1, const void *s2){
159 real *ss1, *ss2;
160 ss1 = (real*) s1;
161 ss2 = (real*) s2;
162
163 if ((ss1)[0] > (ss2)[0]){
164 return 1;
165 } else if ((ss1)[0] < (ss2)[0]){
166 return -1;
167 }
168 return 0;
169}
170
171int comp_descend(const void *s1, const void *s2){
172 real *ss1, *ss2;
173 ss1 = (real*) s1;
174 ss2 = (real*) s2;
175
176 if ((ss1)[0] > (ss2)[0]){
177 return -1;
178 } else if ((ss1)[0] < (ss2)[0]){
179 return 1;
180 }
181 return 0;
182}
183int comp_descend_int(const void *s1, const void *s2){
184 int *ss1, *ss2;
185 ss1 = (int*) s1;
186 ss2 = (int*) s2;
187
188 if ((ss1)[0] > (ss2)[0]){
189 return -1;
190 } else if ((ss1)[0] < (ss2)[0]){
191 return 1;
192 }
193 return 0;
194}
195
196int comp_ascend_int(const void *s1, const void *s2){
197 int *ss1, *ss2;
198 ss1 = (int*) s1;
199 ss2 = (int*) s2;
200
201 if ((ss1)[0] > (ss2)[0]){
202 return 1;
203 } else if ((ss1)[0] < (ss2)[0]){
204 return -1;
205 }
206 return 0;
207}
208
209
210void vector_ordering(int n, real *v, int **p, int ascending){
211 /* give the position of the lagest, second largest etc in vector v if ascending = FALSE
212
213 or
214
215 give the position of the smallest, second smallest etc in vector v if ascending = TRUE.
216 results in p. If *p == NULL, p is assigned.
217
218 ascending: TRUE if v[p] is from small to large.
219 */
220
221 real *u;
222 int i;
223
224 if (!*p) *p = MALLOC(sizeof(int)*n);
225 u = MALLOC(sizeof(real)*2*n);
226
227 for (i = 0; i < n; i++) {
228 u[2*i+1] = i;
229 u[2*i] = v[i];
230 }
231
232 if (ascending){
233 qsort(u, n, sizeof(real)*2, comp_ascend);
234 } else {
235 qsort(u, n, sizeof(real)*2, comp_descend);
236 }
237
238 for (i = 0; i < n; i++) (*p)[i] = (int) u[2*i+1];
239 FREE(u);
240
241}
242
243void vector_sort_real(int n, real *v, int ascending){
244 if (ascending){
245 qsort(v, n, sizeof(real), comp_ascend);
246 } else {
247 qsort(v, n, sizeof(real), comp_descend);
248 }
249}
250void vector_sort_int(int n, int *v, int ascending){
251 if (ascending){
252 qsort(v, n, sizeof(int), comp_ascend_int);
253 } else {
254 qsort(v, n, sizeof(int), comp_descend_int);
255 }
256}
257
258int excute_system_command3(char *s1, char *s2, char *s3){
259 char c[1000];
260
261 strcpy(c, s1);
262 strcat(c, s2);
263 strcat(c, s3);
264 return system(c);
265}
266
267int excute_system_command(char *s1, char *s2){
268 char c[1000];
269
270 strcpy(c, s1);
271 strcat(c, s2);
272 return system(c);
273}
274
275real distance_cropped(real *x, int dim, int i, int j){
276 int k;
277 real dist = 0.;
278 for (k = 0; k < dim; k++) dist += (x[i*dim+k] - x[j*dim + k])*(x[i*dim+k] - x[j*dim + k]);
279 dist = sqrt(dist);
280 return MAX(dist, MINDIST);
281}
282
283real distance(real *x, int dim, int i, int j){
284 int k;
285 real dist = 0.;
286 for (k = 0; k < dim; k++) dist += (x[i*dim+k] - x[j*dim + k])*(x[i*dim+k] - x[j*dim + k]);
287 dist = sqrt(dist);
288 return dist;
289}
290
291real point_distance(real *p1, real *p2, int dim){
292 int i;
293 real dist;
294 dist = 0;
295 for (i = 0; i < dim; i++) dist += (p1[i] - p2[i])*(p1[i] - p2[i]);
296 return sqrt(dist);
297}
298
299char *strip_dir(char *s){
300 int i, first = TRUE;
301 if (!s) return s;
302 for (i = strlen(s); i >= 0; i--) {
303 if (first && s[i] == '.') {/* get rid of .mtx */
304 s[i] = '\0';
305 first = FALSE;
306 }
307 if (s[i] == '/') return (char*) &(s[i+1]);
308 }
309 return s;
310}
311
312void scale_to_box(real xmin, real ymin, real xmax, real ymax, int n, int dim, real *x){
313 real min[3], max[3], min0[3], ratio = 1;
314 int i, k;
315
316 for (i = 0; i < dim; i++) {
317 min[i] = x[i];
318 max[i] = x[i];
319 }
320
321 for (i = 0; i < n; i++){
322 for (k = 0; k < dim; k++) {
323 min[k] = MIN(x[i*dim+k], min[k]);
324 max[k] = MAX(x[i*dim+k], max[k]);
325 }
326 }
327
328 if (max[0] - min[0] != 0) {
329 ratio = (xmax-xmin)/(max[0] - min[0]);
330 }
331 if (max[1] - min[1] != 0) {
332 ratio = MIN(ratio, (ymax-ymin)/(max[1] - min[1]));
333 }
334
335 min0[0] = xmin;
336 min0[1] = ymin;
337 min0[2] = 0;
338 for (i = 0; i < n; i++){
339 for (k = 0; k < dim; k++) {
340 x[i*dim+k] = min0[k] + (x[i*dim+k] - min[k])*ratio;
341 }
342 }
343
344
345}
346
347int digitsQ(char *s){
348 while (*s && *s - '0' >= 0 && *s - '0' <= 9) {
349 s++;
350 }
351 if (*s) return 0;
352 return 1;
353}
354int validQ_int_string(char *to_convert, int *v){
355 /* check to see if this is a string is integer */
356 char *p = to_convert;
357 uint64_t val;
358 errno = 0;
359 val = strtoul(to_convert, &p, 10);
360 if (errno != 0 ||// conversion failed (EINVAL, ERANGE)
361 to_convert == p || // conversion failed (no characters consumed)
362 *p != 0
363 ) return 0;
364 if (val > INT_MAX || val < INT_MIN) return 0;
365 *v = (int) val;
366 return 1;
367}
368