list.c

00001 /* $Id: list.c,v 1.14 2008/02/07 14:39:57 cizzo Exp $
00002  *
00003  * This program is free software; you can redistribute it and/or modify
00004  * it under the terms of the GNU General Public License as published by
00005  * the Free Software Foundation; either version 2 of the License, or
00006  * (at your option) any later version.
00007  *
00008  * This program is distributed in the hope that it will be useful,
00009  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00010  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00011  * GNU General Public License for more details.
00012  *
00013  * You should have received a copy of the GNU General Public License
00014  * along with this program; if not, write to the Free Software
00015  * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301 USA
00016  */
00017 
00018 /*
00019  * $Author: cizzo $
00020  * $Date: 2008/02/07 14:39:57 $
00021  * $Revision: 1.14 $
00022  * $Name:  $
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 /* Same data structure as C++ STL's vector */
00073 struct list
00074 {
00075     void **elements;
00076     int size;
00077     int back;    /* Extra allocated space */
00078 
00079     int current; /* 1 element iteration */
00080 
00081     int current_p1, current_p2; /* pair iteration */
00082 };
00083 
00087 //static void *(*const list_malloc)(size_t) = malloc;
00088 //static void *(*const list_calloc)(size_t, size_t) = calloc;
00089 //static void *(*const list_realloc)(void *, size_t) = realloc;
00090 //static void  (*const list_free)(const void *) = (void (*)(const void *))free;
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  * @brief Destructor
00157  * @param l         list to delete
00158  * @param delete    element copy constructor. If NULL,
00159  *                  elements are not deleted.
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  * @brief Get list size
00188  * @param l         list
00189  * @return  number of elements in list 
00190  *
00191  * Time: O(1)
00192  */
00193 int
00194 list_size(const list *l)
00195 {
00196     assure( l != NULL );
00197 
00198     return l->size;
00199 }
00200 
00201 /*
00202  * @brief Insert element
00203  * @param l         list 
00204  * @param e         element to insert. Must be non-NULL
00205  * @return  number of elements in list 
00206  *
00207  * Time: Amortized O(1)
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  * @brief Remove element
00228  * @param l         list
00229  * @param e         element to remove. Note: pointer comparison
00230  * @return e
00231  *
00232  * The provided element must exist in the list.
00233  * Only one occurrence of e is removed.
00234  * The element is removed from the list, but not deallocated.
00235  *
00236  * For convenience, the function returns the provided element pointer.
00237  * This is to allow code like
00238  * @code
00239  *     element *e = element_list_remove(l, element_list_first(l));
00240  * @endcode
00241  *
00242  * Time: Worst case O(n), O(1) for removing the element returned by list_first()
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         /* Note: amortized constant time */
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  * @brief Iterate
00287  * @param l         list
00288  * @return first element, or NULL if list empty
00289  *
00290  * @code
00291  *     for(element *e = list_first(l);
00292  *         e != NULL; 
00293  *         e = list_next(l)) {...}
00294  * @endcode
00295  *
00296  * The list must not be modified between calls to list_first() or list_next()
00297  *
00298  * Note: It is not possible to have more simultaneous iterations over the same
00299  * list. This functionality can be achived by duplicating the list.
00300  *
00301  * Time: O(1)
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     /* Loop backwards, faster if user
00311        erases the first element */
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  * @brief Iterate
00325  * @param l         list
00326  * @return next element, or NULL if no more elements
00327  *
00328  * See list_first().
00329  *
00330  * Time: O(1)
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  * @brief Iterate through pairs
00353  * @param l         list
00354  * @param e1        (output) first pair 1st element, or NULL if list
00355                     has less than two elements
00356  * @param e2        (output) first pair 2nd element
00357  *
00358  * The iteration is through the K(n,2) different pairs, i.e. the
00359  * pair (e1,e2) is considered equal to the pair (e2,e1) and only
00360  * visited once.
00361  *
00362  * @code
00363  *     for(list_first_pair(l, &e1, &e2);
00364  *         e1 != NULL; 
00365  *         list_next_pair(l, &e1, &e2))
00366  * @endcode
00367  *
00368  * The list must not be modified between calls to
00369  * list_first_pair() or list_next_pair()
00370  *
00371  * The current position is cached in the list object. Therefore simultaneous
00372  * pair iterations over the same list is not allowed.
00373  *
00374  * However, iterating pairs simultaneously with a 1 element iteration (using
00375  * list_first() and list_next()) is allowed.
00376  *
00377  * Time: O(1)
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  * @brief Iterate through pairs
00417  * @param l         list
00418  * @param e1        (output) next pair 1st element, or NULL if no more pairs
00419  * @param e2        (output) next pair 2nd element, or NULL if no more pairs
00420  *
00421  * See list_first_pair().
00422  *
00423  * Time: O(1)
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  * @brief Find minimum element
00506  * @param l          non-empty list
00507  * @param less_than  comparison function which must return true, iff
00508  *                   its first argument is considered less than the second argument.
00509  *                   The 3rd argument is auxillary data used for the comparison.
00510  *                   The provided function should be a strict ordering (i.e. behave
00511  *                   like <)
00512  * @param data       Auxillary data sent to the comparison function. May be NULL.
00513  * @return minimum element
00514  *
00515  * Time: O(n)
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  * @brief Find minimum element
00538  * @param l          non-empty list
00539  * @param eval       evaluation function
00540  * @param data       Auxillary data sent to the evaluation function. May be NULL.
00541  * @return minimum (as defined by the evaluation function) element
00542  *
00543  * Time: O(n)
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  * @brief Find maximum element
00571  * @param l          non-empty list
00572  * @param eval       evaluation function
00573  * @param data       Auxillary data sent to the evaluation function. May be NULL.
00574  * @return maximum (as defined by the evaluation function) element
00575  *
00576  * Time: O(n)
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  * @brief Find maximum element
00604  * @param l          see list_min()
00605  * @param less_than  see list_min()
00606  * @param data       see list_min()
00607  * @return maximum element
00608  *
00609  * Time: O(n)
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  * @brief Find k'th element
00639  * @param l          non-empty list
00640  * @param k          between 1 and list size, inclusive.
00641  * @param less_than  see list_min()
00642  * @param data       see list_min()
00643  * @return k'th element
00644  *
00645  * Be careful to provide an irreflexive comparison function (i.e.
00646  * x < x must always be false), or this function may not return.
00647  *
00648  * Time: Worst case O(n*n). Average over all input: O(n)
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];//fixme select randomly, swap with 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  * @brief Determine order of elements given an evaluation function
00700  * @param e1         first element
00701  * @param e2         second element
00702  * @param data       containing the evaluation function, and
00703  *                   additional data which is passed to the evaluation function
00704  * @return true iff e1 evaluates to a number less than e2
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     /* Cast is safe, see caller of this function */
00714 
00715     
00716     /* Unfortunately, as the following commented code demonstrated,
00717        with GCC-4.2.0 calling the evaluation function two times with the same
00718        input can give two numbers which satisfy both equality (==) *and*
00719        less than (<) (!) but not greater than (>). This causes the kth function
00720        to loop infinitely.
00721 
00722        Avoid that by handling explicitly this special case.
00723     */
00724 
00725     if (e1 == e2) return false;
00726 
00727     /*
00728     double d1 = d->f(e1, d->aux_data);
00729     fprintf(stderr, "%s got %f %g \n", __func__, d1, d1);
00730     fprintf(stderr, "%.100f\n", d1);
00731 
00732     double d2 = d->f(e2, d->aux_data);
00733     fprintf(stderr, "%s %d %d %d\n", __func__, i1, i2, i3);
00734 
00735     fprintf(stderr, "%s got %f %g %d  %d\n", __func__, d2, d2, d1 < d2, e1 == e2);
00736     fprintf(stderr, "%.100f  %d %d  %d  %d %d\n", d2, 
00737             d1 < d2, d2 > d1, d1 == d2, d1 > d2, d2 < d1);
00738 
00739     fprintf(stderr, "l1 = %ld\n", d1);
00740     fprintf(stderr, "l2 = %ld\n", d2);
00741     */    
00742 
00743     return (d->f(e1, d->aux_data) < d->f(e2, d->aux_data));
00744 }
00745 
00746 /*
00747  * @brief k'th element
00748  * @param    l        list
00749  * @param    k        counting from 1
00750  * @param    eval     returns the value of an element
00751  * @param    data     sent to eval
00752  * @return   a k'th element according to the evaluation function
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  * @brief Compute median
00778  * @param l          list
00779  * @param eval       returns the value of an element
00780  * @param data       additional data passed to eval
00781  * @return median
00782  *
00783  * Time: Average O(n).
00784  *
00785  * For an even number of elements, the median is the
00786  * average of the two central values.
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  * @brief Compute mean
00812  * @param l          list
00813  * @param eval       returns the value of an element
00814  * @param data       additional data passed to eval
00815  * @return arithmetic average
00816  *
00817  * Time: O(n).
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  * @brief Compute optimally weighted mean
00840  * @param l          list with two or more elements
00841  * @param eval       returns the value of an element
00842  * @param data       additional data passed to eval
00843  * @param eval_err   returns the error (1 sigma) of an element,
00844  *                   must be positive
00845  * @param data_err   additional data passed to eval_err
00846  * @param err        (output) square root variance
00847  * @param red_chisq  (output) reduced chi square, may be NULL
00848  * @return optimal mean
00849  *
00850  *  average  =  [ sum  x_i / sigma_i^2 ] / [sum 1 / sigma_i^2]
00851  *  variance =                         1 / [sum 1 / sigma_i^2]
00852  *
00853  *  chi^2/(n-1) = sum ((x_i - average) / sigma_i)^2  / (n-1)
00854  *
00855  * Time: O(n).
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  * @brief Compute absolute deviation
00906  * @param e1        element
00907  * @param data      reference value, and evaluation function
00908  *                  and additional data passed to the evaluation function
00909  * @return absolute difference between e1 and reference point
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     /* Cast is safe, see caller */
00920 
00921     return fabs(d->f(e1, d->aux_data) - d->ref);
00922 }
00923 
00924 /*
00925  * @brief Compute median absolute deviation wrt median
00926  * @param l          list
00927  * @param eval       returns the value of an element
00928  * @param data       additional data passed to eval
00929  * @return mad
00930  *
00931  * Time: Average O(n).
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 

Generated on Wed Sep 10 07:31:52 2008 for FORS Pipeline Reference Manual by  doxygen 1.4.6