Changeset 37208 in project


Ignore:
Timestamp:
02/05/19 08:11:01 (2 weeks ago)
Author:
iraikov
Message:

statistics: replacing gsl dependency with toms/asa implementations

Location:
release/5/statistics/trunk
Files:
8 added
4 edited

Legend:

Unmodified
Added
Removed
  • release/5/statistics/trunk/build.sh

    r36068 r37208  
    11# -*- sh -*-
    22
    3 "$CHICKEN_CSI" -q -s build.scm "$CHICKEN_CSC $@"
     3$CHICKEN_CSC $@ toms743.c toms179.c asa239.c asa245.c
    44
     5
  • release/5/statistics/trunk/statistics.egg

    r36443 r37208  
    66 (category math)
    77 (build-dependencies compile-file srfi-13 )
    8  (dependencies srfi-1 srfi-25 srfi-69 vector-lib yasos)
     8 (dependencies srfi-1 srfi-25 srfi-69 vector-lib random-mtzig yasos)
    99 (test-dependencies test srfi-1)
    1010 (components (extension statistics (custom-build "build.sh"))))
  • release/5/statistics/trunk/statistics.scm

    r37172 r37208  
    8585    square
    8686    cumsum
     87    set-random-seed!
     88    random-state
    8789   
    8890    ; descriptive statistics
     
    180182          (chicken keyword) (prefix (chicken sort) list.)
    181183          (prefix (only srfi-1 fold iota filter find delete-duplicates every first last) list.)
    182           (only srfi-13 string<) srfi-25 srfi-69 vector-lib
     184          srfi-4 (only srfi-13 string<) srfi-25 srfi-69 vector-lib
     185          (prefix random-mtzig random:)
    183186          yasos yasos-collections)
    184187
     
    187190
    188191  #>
    189   #include <assert.h>
    190   #include <gsl/gsl_rng.h>
    191   #include <gsl/gsl_sf_erf.h>
    192   #include <gsl/gsl_sf_gamma.h>
    193   #include <gsl/gsl_sf_lambert.h>
     192  #include "asa239.h"
     193  #include "asa245.h"
     194  #include "toms179.h"
     195  #include "toms743.h"
     196  <#
     197
    194198 
    195   const gsl_rng_type *statistics_rng_T;
    196   gsl_rng *statistics_rng_st;
    197 
    198   void statistics_rng_init()
    199   {
    200     gsl_rng_env_setup();
    201                    
    202     statistics_rng_T = gsl_rng_default;
    203     assert(statistics_rng_T != NULL);
    204     statistics_rng_st = gsl_rng_alloc (statistics_rng_T);
    205     assert(statistics_rng_st != NULL);
    206     }
    207 
    208   gsl_rng *statistics_rng_get_state ()
    209   {
    210      return statistics_rng_st;
    211   }
    212   <#
    213 
    214   (define rng-init (foreign-lambda void "statistics_rng_init"))
    215   (define rng-state (foreign-lambda nonnull-c-pointer "statistics_rng_get_state"))
    216   (define gsl-random-uniform (foreign-lambda double "gsl_rng_uniform" nonnull-c-pointer))
    217   (define gsl-random-uniform-int (foreign-lambda unsigned-int "gsl_rng_uniform_int" nonnull-c-pointer unsigned-int))
    218 
    219   ;; The principal branch of the Lambert W function
    220   (define lambert-W0
    221     (foreign-lambda double "gsl_sf_lambert_W0" double))
    222 
     199  (define random-state (make-parameter (random:init)))
     200  (define (set-random-seed! x) (random-state (random:init x)))
     201
     202  (define (random-uniform)
     203    (random:randu! (random-state)))
     204
     205  (define (random-uniform-int n)
     206    (modulo (random:random! (random-state)) n))
     207 
     208  ;; Approximation to the Lambert W function
     209  (define wapr
     210    (foreign-lambda double "wapr" double int s32vector int))
     211
     212  (define (lambert-W0 x)
     213    (let ((err (s32vector 0)))
     214      (wapr x 0 err 0)))
     215 
    223216  ;;The secondary real-valued branch of the Lambert W function
    224   (define lambert-Wm1
    225     (foreign-lambda double "gsl_sf_lambert_Wm1" double))
    226 
    227   (define beta-incomplete
    228     (foreign-lambda double "gsl_sf_beta_inc" double double double))
    229 
    230   (define (random-uniform)
    231     (gsl-random-uniform (rng-state)))
    232 
    233   (define (random-uniform-int n)
    234     (gsl-random-uniform-int (rng-state) n))
    235 
     217  (define (lambert-Wm1 x)
     218    (let ((err (s32vector 0)))
     219      (wapr x 1 err 0)))
     220
     221
     222  (define mdbeta
     223    (foreign-lambda double "mdbeta" double double double s32vector))
     224
     225  (define (beta-incomplete p q x)
     226    (let ((ierr (s32vector 0)))
     227      (mdbeta x p q ierr)))
     228 
     229  (define gammad (foreign-lambda double "gammad" double double s32vector))
     230   
    236231  ;; returns complementary normalised incomplete Gamma Function
    237232  (define (gamma-incomplete a x)
    238     (define gamma-inc (foreign-lambda double "gsl_sf_gamma_inc_P" double double))
    239     (values (if (= x 0.0) 0.0 (gamma-inc a x))
    240             (gamma-ln a)))
     233    (let ((ierr (s32vector 0)))
     234      (values (if (= x 0.0) 0.0 (gammad x a ierr))
     235              (gamma-ln a))))
     236
     237  (define alngam (foreign-lambda double "alngam" double s32vector))
    241238
    242239  (define (gamma-ln x)
    243     (define lngamma (foreign-lambda double "gsl_sf_lngamma" double))
    244     (cond ((<= x 0)
    245            (error "Argument to gamma-ln must be positive"))
    246           ((> x 1.0e302)
    247            (error "Argument to gamma-ln too large"))
    248           (else
    249             (lngamma x))))
     240    (let ((ierr (s32vector 0)))
     241      (cond ((<= x 0)
     242             (error "Argument to gamma-ln must be positive"))
     243            ((> x 1.0e302)
     244             (error "Argument to gamma-ln too large"))
     245            (else
     246             (alngam x ierr)))))
    250247
    251248  ;; ---------------------------------------------------------------------------
     
    16001597                   tail-area))))))
    16011598
    1602   (rng-init)
    16031599 
    16041600  ) ; end of library
  • release/5/statistics/trunk/tests/run.scm

    r37172 r37208  
    6565            (test "gamma-incomplete"
    6666                  (list 0.44217 0.0)
    67                   (let-values (((a b) (gamma-incomplete 2.0 1.5))) (list (to-5-dp a) b)))
     67                  (let-values (((a b) (gamma-incomplete 2.0 1.5)))
     68                    (list (to-5-dp a) b)))
    6869           
    6970            (test-assert "gamma-ln"
     
    329330            (test-assert "lambert-W0" (=5 (lambert-W0 1.0) 0.567143290410))
    330331           
    331             (test-assert "lambert-Wm1" (=5 (lambert-Wm1 1.0) 0.567143290410))
    332             )
     332            (test-assert "lambert-Wm1" (=5 (lambert-Wm1 -0.1) -3.57715206))
     333            )
     334
    333335
    334336(test-exit)
Note: See TracChangeset for help on using the changeset viewer.