Return to ko-perturbation.c CVS log

Up to [local] / OpenXM / src / hgm / fisher-bingham / src

Annotation of OpenXM/src/hgm/fisher-bingham/src/ko-perturbation.c, Revision 1.1

1.1     ! takayama    1: #include <stdio.h>
        !             2: #include <stdlib.h>
        !             3: #include <gsl/gsl_rng.h>
        !             4: #include <gsl/gsl_randist.h>
        !             5:
        !             6: static int init=1;
        !             7: static const gsl_rng_type *T;
        !             8: static gsl_rng * r;
        !             9: static int length;
        !            10: static double eps;
        !            11:
        !            12: double
        !            13: *perturbation(const double v[], double err[])
        !            14: {
        !            15:   if(init){
        !            16:     fprintf(stderr, "perturbation: erorr! Forget perturbation_init?\n");
        !            17:     exit(1);
        !            18:   }
        !            19:
        !            20:   int i;
        !            21:   double *vp = (double *) malloc(sizeof(double) * length);
        !            22:   for(i=0; i<length; i++){
        !            23:     err[i]= gsl_ran_gaussian (r, eps);
        !            24:     vp[i] = v[i]+err[i];
        !            25:   }
        !            26:   return vp;
        !            27: }
        !            28:
        !            29: void
        !            30: perturbation_init(int l,  double e)
        !            31: {
        !            32:   if(init){
        !            33:     gsl_rng_env_setup();
        !            34:     T = gsl_rng_default;
        !            35:     r = gsl_rng_alloc (T);
        !            36:     init = 0;
        !            37:   }
        !            38:   length = l;
        !            39:   eps = e;
        !            40: }
        !            41:
        !            42: void
        !            43: perturbation_free(void)
        !            44: {
        !            45:   gsl_rng_free (r);
        !            46: }