00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #ifdef HAVE_CONFIG_H
00026 #include <config.h>
00027 #endif
00028
00029 #include <list_void.h>
00030
00031 #include <math.h>
00032 #include <stdbool.h>
00033 #include <stdlib.h>
00034 #include <stdio.h>
00035
00072
00073 struct list
00074 {
00075 void **elements;
00076 int size;
00077 int back;
00078
00079 int current;
00080
00081 int current_p1, current_p2;
00082 };
00083
00087
00088
00089
00090
00091 static void *(*list_malloc)(size_t) = malloc;
00092 static void *(*list_calloc)(size_t, size_t) = calloc;
00093 static void *(*list_realloc)(void *, size_t) = realloc;
00094 static void (*list_free)(const void *) = (void (*)(const void *))free;
00095
00099 #include <assert.h>
00100 #define assure(EXPR) assert(EXPR)
00101
00106 list *
00107 list_new(void)
00108 {
00109 list *l = list_malloc(sizeof(*l));
00110
00111 l->elements = NULL;
00112 l->size = 0;
00113 l->back = 0;
00114
00115 return l;
00116 }
00117
00124 list *
00125 list_duplicate(const list *l, void * (*duplicate)(const void *))
00126 {
00127 assure( l != NULL );
00128
00129 {
00130 list *dupl = list_malloc(sizeof(*dupl));
00131
00132 dupl->elements = list_malloc((l->size+l->back) * sizeof(*dupl->elements));
00133 dupl->size = l->size;
00134 dupl->back = l->back;
00135 dupl->current = l->current;
00136 dupl->current_p1 = l->current_p1;
00137 dupl->current_p2 = l->current_p2;
00138
00139 {
00140 int i;
00141 for (i = 0; i < l->size; i++) {
00142 if (duplicate != NULL) {
00143 dupl->elements[i] = duplicate(l->elements[i]);
00144 }
00145 else {
00146 dupl->elements[i] = l->elements[i];
00147 }
00148 }
00149 }
00150
00151 return dupl;
00152 }
00153 }
00154
00155
00156
00157
00158
00159
00160
00161 void
00162 list_delete_const(const list **l, void (*delete)(void **))
00163 {
00164 if (l != NULL && *l != NULL) {
00165
00166 if (delete != NULL) {
00167
00168 int i;
00169 for (i = 0; i < (*l)->size; i++) {
00170 delete(&((*l)->elements[i]));
00171 }
00172 }
00173 list_free((*l)->elements);
00174 list_free((*l)); *l = NULL;
00175 }
00176 return;
00177 }
00178
00179 void
00180 list_delete(list **l, void (*delete)(void **))
00181 {
00182 list_delete_const((const list **)l, delete);
00183 return;
00184 }
00185
00186
00187
00188
00189
00190
00191
00192
00193 int
00194 list_size(const list *l)
00195 {
00196 assure( l != NULL );
00197
00198 return l->size;
00199 }
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209 void
00210 list_insert(list *l, void *e)
00211 {
00212 assure( e != NULL );
00213
00214 if (l->back == 0) {
00215 l->back = l->size + 1;
00216 l->elements = list_realloc(l->elements, (l->size + l->back) * sizeof(*l->elements));
00217 }
00218
00219 l->size++;
00220 l->back--;
00221 l->elements[l->size - 1] = e;
00222
00223 return;
00224 }
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244 const void *
00245 list_remove_const(list *l, const void *e)
00246 {
00247 assure( l != NULL );
00248 assure( e != NULL );
00249
00250 {
00251 int i;
00252 int indx = -1;
00253 for (i = l->size - 1; i >= 0 && indx < 0; i--) {
00254 if (l->elements[i] == e) {
00255 indx = i;
00256 }
00257 }
00258
00259 assure( indx >= 0 );
00260
00261 for (i = indx; i < l->size-1; i++) {
00262 l->elements[i] = l->elements[i+1];
00263 }
00264 }
00265
00266 l->size--;
00267 l->back++;
00268
00269 if (l->back > 4 * l->size) {
00270
00271 l->back = l->size;
00272 l->elements = list_realloc(l->elements,
00273 (l->size + l->back) * sizeof(*l->elements));
00274 }
00275
00276 return e;
00277 }
00278
00279 void *
00280 list_remove(list *l, void *e)
00281 {
00282 return (void *)list_remove_const(l, e);
00283 }
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303 const void *
00304 list_first_const(const list *l)
00305 {
00306 assure( l != NULL );
00307
00308 if (l->size == 0) return NULL;
00309
00310
00311
00312
00313 ((list *)l)->current = l->size - 1;
00314 return l->elements[l->current];
00315 }
00316
00317 void *
00318 list_first(list *l)
00319 {
00320 return (void *)list_first_const(l);
00321 }
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332 const void *
00333 list_next_const(const list *l)
00334 {
00335 assure( l != NULL );
00336
00337 if (l->size == 0) return NULL;
00338
00339 ((list *)l)->current -= 1;
00340
00341 if (l->current < 0) return NULL;
00342 else return l->elements[l->current];
00343 }
00344
00345 void *
00346 list_next(list *l)
00347 {
00348 return (void *)list_next_const(l);
00349 }
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379 void
00380 list_first_pair_const(const list *l,
00381 const void **e1,
00382 const void **e2)
00383 {
00384 assure( l != NULL );
00385 assure( e1 != NULL );
00386 assure( e2 != NULL );
00387
00388 if (l->size <= 1) {
00389 *e1 = NULL;
00390 *e2 = NULL;
00391 return;
00392 }
00393
00394 ((list *)l)->current_p1 = l->size - 1;
00395 ((list *)l)->current_p2 = l->size - 2;
00396
00397 *e1 = l->elements[l->current_p1];
00398 *e2 = l->elements[l->current_p2];
00399
00400 return;
00401 }
00402
00403 void
00404 list_first_pair(list *l,
00405 void **e1,
00406 void **e2)
00407 {
00408 list_first_pair_const(l,
00409 (const void **)e1,
00410 (const void **)e2);
00411
00412 return;
00413 }
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425 void
00426 list_next_pair_const(const list *l,
00427 const void **e1,
00428 const void **e2)
00429 {
00430 assure( l != NULL );
00431 assure( e1 != NULL );
00432 assure( e2 != NULL );
00433
00434 if (l->size <= 1) {
00435 *e1 = NULL;
00436 *e2 = NULL;
00437 return;
00438 }
00439
00440 ((list *)l)->current_p2 -= 1;
00441
00442 if (l->current_p2 < 0) {
00443 ((list *)l)->current_p1 -= 1;
00444 ((list *)l)->current_p2 = l->current_p1 - 1;
00445 if (l->current_p2 < 0) {
00446 *e1 = NULL;
00447 *e2 = NULL;
00448 return;
00449 }
00450 *e1 = l->elements[l->current_p1];
00451 *e2 = l->elements[l->current_p2];
00452 return;
00453 }
00454
00455 *e2 = l->elements[l->current_p2];
00456 return;
00457 }
00458
00459 void
00460 list_next_pair(list *l,
00461 void **e1,
00462 void **e2)
00463 {
00464 list_next_pair_const(l, (const void **)e1, (const void **)e2);
00465 return;
00466 }
00467
00480 list *
00481 list_extract(const list *l,
00482 void *(*duplicate)(const void *),
00483 bool (*predicate)(const void *, void *),
00484 void *data)
00485 {
00486 assure( l != NULL );
00487 assure( duplicate != NULL);
00488 assure( predicate != NULL);
00489
00490 {
00491 list *ex = list_new();
00492 int i;
00493
00494 for (i = 0; i < l->size; i++) {
00495 if (predicate(l->elements[i], data)) {
00496 list_insert(ex, duplicate(l->elements[i]));
00497 }
00498 }
00499
00500 return ex;
00501 }
00502 }
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517 void *
00518 list_min(list *l, list_func_lt less_than, void *data)
00519 {
00520 assure( l != NULL );
00521 assure( less_than != NULL);
00522 assure( list_size(l) > 0);
00523
00524 {
00525 int minindex = 0;
00526 int i;
00527 for (i = 1; i < l->size; i++) {
00528 if (less_than(l->elements[i], l->elements[minindex], data))
00529 minindex = i;
00530 }
00531
00532 return l->elements[minindex];
00533 }
00534 }
00535
00536
00537
00538
00539
00540
00541
00542
00543
00544
00545 void *
00546 list_min_val(list *l, list_func_eval eval, void *data)
00547 {
00548 assure( l != NULL );
00549 assure( eval != NULL);
00550 assure( list_size(l) > 0);
00551
00552 {
00553 int minindex = 0;
00554 double minval = eval(l->elements[0], data);
00555 int i;
00556
00557 for (i = 1; i < l->size; i++) {
00558 double val = eval(l->elements[i], data);
00559 if (val < minval) {
00560 minval = val;
00561 minindex = i;
00562 }
00563 }
00564
00565 return l->elements[minindex];
00566 }
00567 }
00568
00569
00570
00571
00572
00573
00574
00575
00576
00577
00578 void *
00579 list_max_val(list *l, list_func_eval eval, void *data)
00580 {
00581 assure( l != NULL );
00582 assure( eval != NULL);
00583 assure( list_size(l) > 0);
00584
00585 {
00586 int maxindex = 0;
00587 double maxval = eval(l->elements[0], data);
00588 int i;
00589
00590 for (i = 1; i < l->size; i++) {
00591 double val = eval(l->elements[i], data);
00592 if (val > maxval) {
00593 maxval = val;
00594 maxindex = i;
00595 }
00596 }
00597
00598 return l->elements[maxindex];
00599 }
00600 }
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611 const void *
00612 list_max_const(const list *l, list_func_lt less_than, void *data)
00613 {
00614 assure( l != NULL );
00615 assure( less_than != NULL);
00616 assure( list_size(l) > 0);
00617
00618 {
00619 int maxindex = 0;
00620 int i;
00621 for (i = 1; i < l->size; i++) {
00622 if (!less_than(l->elements[i], l->elements[maxindex], data))
00623 maxindex = i;
00624 }
00625
00626 return l->elements[maxindex];
00627 }
00628 }
00629
00630 void *
00631 list_max(list *l, list_func_lt less_than, void *data)
00632 {
00633 return (void *)list_max_const(l, less_than, data);
00634 }
00635
00636
00637
00638
00639
00640
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650 static const void *
00651 kth(const void *a[], int k, int n, list_func_lt less_than, void *data)
00652 {
00653 int i, j, lo, hi;
00654
00655 k -= 1;
00656
00657 lo = 0;
00658 hi = n - 1;
00659 while (lo < hi) {
00660 const void *pivot = a[k];
00661 i = lo;
00662 j = hi;
00663 do {
00664 while (less_than(a[i], pivot, data)) i++;
00665 while (less_than(pivot, a[j], data)) {
00666 j--;
00667 }
00668 if (i <= j) {
00669 const void *tmp = a[i];
00670 a[i] = a[j];
00671 a[j] = tmp;
00672 i++; j--;
00673 }
00674 } while (i <= j);
00675 if (j < k) lo = i;
00676 if (k < i) hi = j;
00677 }
00678 return a[k];
00679 }
00680 const void *
00681 list_kth_const(const list *l, int k,
00682 list_func_lt less_than, void *data)
00683 {
00684 assure( l != NULL );
00685 assure( 1 <= k && k <= l->size );
00686
00687 return kth((const void **)l->elements, k, l->size, less_than, data);
00688 }
00689
00690 void *
00691 list_kth(list *l, int k,
00692 list_func_lt less_than,
00693 void *data)
00694 {
00695 return (void *)list_kth_const(l, k, less_than, data);
00696 }
00697
00698
00699
00700
00701
00702
00703
00704
00705
00706 static
00707 bool val_less_than(const void *e1, const void *e2, void *data)
00708 {
00709 struct {
00710 list_func_eval f;
00711 void *aux_data;
00712 } *d = data;
00713
00714
00715
00716
00717
00718
00719
00720
00721
00722
00723
00724
00725 if (e1 == e2) return false;
00726
00727
00728
00729
00730
00731
00732
00733
00734
00735
00736
00737
00738
00739
00740
00741
00742
00743 return (d->f(e1, d->aux_data) < d->f(e2, d->aux_data));
00744 }
00745
00746
00747
00748
00749
00750
00751
00752
00753
00754 const void *
00755 list_kth_val_const(const list *l, int k, list_func_eval eval, void *data)
00756 {
00757 assure( l != NULL );
00758 assure( 1 <= k && k <= l->size );
00759 assure( eval != NULL );
00760
00761 struct {
00762 list_func_eval f;
00763 void *aux_data;
00764 } d;
00765 d.f = eval;
00766 d.aux_data = data;
00767
00768 return list_kth_const(l, k, val_less_than, &d);
00769 }
00770 void *
00771 list_kth_val(list *l, int k, list_func_eval eval, void *data)
00772 {
00773 return (void *) list_kth_val_const(l, k, eval, data);
00774 }
00775
00776
00777
00778
00779
00780
00781
00782
00783
00784
00785
00786
00787
00788 double
00789 list_median(const list *l, list_func_eval eval, void *data)
00790 {
00791 assure( l != NULL );
00792 assure( eval != NULL );
00793 assure( l->size > 0 );
00794
00795 const void *median_element = list_kth_val_const(l, (l->size+1)/2, eval, data);
00796
00797 double median_val = eval(median_element, data);
00798
00799 if (list_size(l) && 2 == 0)
00800 {
00801 const void *other_median_element =
00802 list_kth_val_const(l, (l->size+2)/2, eval, data);
00803
00804 median_val = (median_val + eval(other_median_element, data) ) / 2.0;
00805 }
00806
00807 return median_val;
00808 }
00809
00810
00811
00812
00813
00814
00815
00816
00817
00818
00819 double
00820 list_mean(const list *l, list_func_eval eval, void *data)
00821 {
00822 assure( l != NULL );
00823 assure( eval != NULL );
00824 assure( l->size > 0 );
00825
00826 double result = eval(l->elements[0], data);
00827 int i;
00828
00829 for (i = 1; i < l->size; i++) {
00830 result += eval(l->elements[i], data);
00831 }
00832
00833 result /= l->size;
00834
00835 return result;
00836 }
00837
00838
00839
00840
00841
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857 double
00858 list_mean_optimal(const list *l,
00859 list_func_eval eval, void *data_eval,
00860 list_func_eval eval_err, void *data_err,
00861 double *err,
00862 double *red_chisq)
00863 {
00864 assure( l != NULL );
00865 assure( l->size >= 1 );
00866 assure( red_chisq == NULL || l->size >= 2 );
00867 assure( eval != NULL );
00868 assure( eval_err != NULL );
00869 assure( err != NULL );
00870
00871 double sum_weights = 0;
00872 double opt_average = 0;
00873 int i;
00874
00875 for (i = 0; i < l->size; i++) {
00876 void *e = l->elements[i];
00877 double sigma = eval_err(e, data_err);
00878
00879 assure( sigma > 0 );
00880
00881 double weight = 1/(sigma*sigma);
00882
00883 opt_average += eval(e, data_eval) * weight;
00884 sum_weights += weight;
00885 }
00886 opt_average /= sum_weights;
00887 *err = 1 / sqrt(sum_weights);
00888
00889 if (red_chisq != NULL) {
00890 *red_chisq = 0;
00891 for (i = 0; i < l->size; i++) {
00892 void *e = l->elements[i];
00893 double residual = ((eval(e, data_eval) - opt_average)) / eval_err(e, data_err);
00894 *red_chisq += residual * residual;
00895 }
00896 *red_chisq /= (l->size - 1);
00897 }
00898
00899 return opt_average;
00900 }
00901
00902
00903
00904
00905
00906
00907
00908
00909
00910
00911 static double
00912 abs_dev(const void *e1, void *data)
00913 {
00914 struct {
00915 double ref;
00916 list_func_eval f;
00917 void *aux_data;
00918 } *d = data;
00919
00920
00921 return fabs(d->f(e1, d->aux_data) - d->ref);
00922 }
00923
00924
00925
00926
00927
00928
00929
00930
00931
00932
00933 double
00934 list_mad(list *l, list_func_eval eval, void *data)
00935 {
00936 assure( l != NULL );
00937 assure( eval != NULL );
00938
00939 double median = list_median(l, eval, data);
00940
00941 struct {
00942 double ref;
00943 list_func_eval f;
00944 void *aux_data;
00945 } d;
00946
00947 d.ref = median;
00948 d.f = eval;
00949 d.aux_data = data;
00950
00951 return list_median(l, abs_dev, &d);
00952 }
00953