Implement RNG using FFI instead of C implementation - gsl-shell.git - gsl-shell

diff --git a/Makefile b/Makefile
index 8a7a4f07..772fef03 100644
--- a/Makefile
+++ b/Makefile

@@ -45,8 +45,7 @@ endif

SUBDIRS = $(LUADIR)

-C_SRC_FILES = str.c gs-types.c lua-utils.c lua-rng.c randist.c sf.c \

- lua-graph.c lua-gsl.c

+C_SRC_FILES = str.c gs-types.c lua-utils.c sf.c lua-graph.c lua-gsl.c

ifeq ($(strip $(USE_READLINE)),yes)

C_SRC_FILES += completion.c

diff --git a/gsl.lua b/gsl.lua
index e0a3864c..d4e67617 100644
--- a/gsl.lua
+++ b/gsl.lua

@@ -1522,8 +1522,55 @@ ffi.cdef[[

double * chisq,

gsl_multifit_linear_workspace * work);

-struct _gsl_rng;

-typedef struct _gsl_rng gsl_rng;

+typedef struct

+ {

+ const char *name;

+ unsigned long int max;

+ unsigned long int min;

+ size_t size;

+ void (*set) (void *state, unsigned long int seed);

+ unsigned long int (*get) (void *state);

+ double (*get_double) (void *state);

+ }

+gsl_rng_type;

+

+typedef struct

+ {

+ const gsl_rng_type * type;

+ void *state;

+ }

+gsl_rng;

+

+const gsl_rng_type ** gsl_rng_types_setup(void);

+

+const gsl_rng_type *gsl_rng_default;

+unsigned long int gsl_rng_default_seed;

+

+gsl_rng *gsl_rng_alloc (const gsl_rng_type * T);

+int gsl_rng_memcpy (gsl_rng * dest, const gsl_rng * src);

+gsl_rng *gsl_rng_clone (const gsl_rng * r);

+

+void gsl_rng_free (gsl_rng * r);

+

+void gsl_rng_set (const gsl_rng * r, unsigned long int seed);

+unsigned long int gsl_rng_max (const gsl_rng * r);

+unsigned long int gsl_rng_min (const gsl_rng * r);

+const char *gsl_rng_name (const gsl_rng * r);

+

+int gsl_rng_fread (FILE * stream, gsl_rng * r);

+int gsl_rng_fwrite (FILE * stream, const gsl_rng * r);

+

+size_t gsl_rng_size (const gsl_rng * r);

+void * gsl_rng_state (const gsl_rng * r);

+

+void gsl_rng_print_state (const gsl_rng * r);

+

+const gsl_rng_type * gsl_rng_env_setup (void);

+

+unsigned long int gsl_rng_get (const gsl_rng * r);

+double gsl_rng_uniform (const gsl_rng * r);

+double gsl_rng_uniform_pos (const gsl_rng * r);

+unsigned long int gsl_rng_uniform_int (const gsl_rng * r, unsigned long int n);

unsigned int gsl_ran_bernoulli (const gsl_rng * r, double p);

double gsl_ran_bernoulli_pdf (const unsigned int k, double p);

diff --git a/gslext.lua b/gslext.lua
index db16ee9f..79adceee 100644
--- a/gslext.lua
+++ b/gslext.lua

@@ -3,6 +3,8 @@

require('iter')

require('matrix')

require('num')

+require('rng')

+require('rnd')

require('integ-init')

require('fft-init')

require('graph-init')

diff --git a/lua-gsl.c b/lua-gsl.c
index 82687f11..cd850170 100644
--- a/lua-gsl.c
+++ b/lua-gsl.c

@@ -21,11 +21,12 @@

#include <lua.h>

#include <lauxlib.h>

+#include <gsl/gsl_types.h>

+#include <gsl/gsl_errno.h>

+

#include "lua-gsl.h"

#include "gs-types.h"

#include "lua-utils.h"

-#include "lua-rng.h"

-#include "randist.h"

#include "sf.h"

#include "lua-graph.h"

@@ -53,8 +54,6 @@ luaopen_gsl (lua_State *L)

lua_pop (L, 1);

register_graph (L);

- rng_register (L);

- randist_register (L);

sf_register (L);

return 1;

diff --git a/lua-rng.c b/lua-rng.c
deleted file mode 100644
index 248b4788..00000000
--- a/lua-rng.c
+++ /dev/null

@@ -1,169 +0,0 @@

-

-/* random.c

- *

- * Copyright (C) 2009 Francesco Abbate

- *

- * This program is free software; you can redistribute it and/or modify

- * it under the terms of the GNU General Public License as published by

- * the Free Software Foundation; either version 3 of the License, or (at

- * your option) any later version.

- *

- * This program is distributed in the hope that it will be useful, but

- * WITHOUT ANY WARRANTY; without even the implied warranty of

- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU

- * General Public License for more details.

- *

- * You should have received a copy of the GNU General Public License

- * along with this program; if not, write to the Free Software

- * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.

- */

-

-#include <lua.h>

-#include <lauxlib.h>

-#include <string.h>

-#include <math.h>

-#include <gsl/gsl_rng.h>

-

-#include "gs-types.h"

-#include "lua-rng.h"

-

-static int rng_free (lua_State *L);

-static int rng_new (lua_State *L);

-static int rng_type_list (lua_State *L);

-

-static int rng_get (lua_State *L);

-static int rng_getint (lua_State *L);

-static int rng_set (lua_State *L);

-

-static gsl_rng *push_rng (lua_State *L, const gsl_rng_type *T);

-

-static const struct luaL_Reg rng_functions[] = {

- {"new", rng_new},

- {"list", rng_type_list},

- {NULL, NULL}

-};

-

-static const struct luaL_Reg rng_methods[] = {

- {"get", rng_get},

- {"getint", rng_getint},

- {"set", rng_set},

- {"__gc", rng_free},

- {NULL, NULL}

-};

-

-int

-rng_free (lua_State *L)

-{

- struct lua_rng *rng_udata = gs_check_userdata (L, 1, GS_RNG);

- gsl_rng_free (rng_udata->rng);

- return 0;

-}

-

-int

-rng_type_list (lua_State *L)

-{

- const gsl_rng_type **t, **t0;

- size_t k;

-

- t0 = gsl_rng_types_setup ();

-

- lua_newtable (L);

- for (t = t0, k = 0; *t != NULL; t++, k++)

- {

- lua_pushstring (L, (*t)->name);

- lua_rawseti (L, -2, k+1);

- }

-

- return 1;

-}

-

-gsl_rng *

-push_rng (lua_State *L, const gsl_rng_type *T)

-{

- struct lua_rng *r;

-

- r = lua_newuserdata (L, sizeof(struct lua_rng));

- r->rng = gsl_rng_alloc (T);

-

- r->min = gsl_rng_min (r->rng);

- r->max = gsl_rng_max (r->rng);

-

- gs_set_metatable (L, GS_RNG);

-

- return r->rng;

-}

-

-int

-rng_new (lua_State *L)

-{

- const char *reqname = luaL_optstring (L, 1, "mt19937");

- const gsl_rng_type **t, **t0;

-

- t0 = gsl_rng_types_setup ();

-

- for (t = t0; *t != NULL; t++)

- {

- if (strcmp ((*t)->name, reqname) == 0)

- {

- push_rng (L, *t);

- return 1;

- }

- }

-

- luaL_error (L, "the requested generator does not exist");

-

- return 0;

-}

-

-int

-rng_get (lua_State *L)

-{

- struct lua_rng *udata = gs_check_userdata (L, 1, GS_RNG);

- double v = gsl_rng_uniform (udata->rng);

- lua_pushnumber (L, v);

- return 1;

-}

-

-int

-rng_set (lua_State *L)

-{

- struct lua_rng *udata = gs_check_userdata (L, 1, GS_RNG);

- unsigned long int seed = luaL_checkinteger (L, 2);

- gsl_rng_set (udata->rng, seed);

- return 0;

-}

-

-int

-rng_getint (lua_State *L)

-{

- struct lua_rng *udata = gs_check_userdata (L, 1, GS_RNG);

- unsigned long int lmt = luaL_checkinteger (L, 2);

- unsigned long int genlmt = udata->max - udata->min;

- unsigned long int j;

- double v = 0.0, vmod;

-

- for (j = lmt; j > 0; j = j / genlmt)

- {

- v *= genlmt;

- v += gsl_rng_uniform_int (udata->rng, genlmt);

- }

-

- vmod = fmod (v, (double) lmt);

-

- lua_pushnumber (L, vmod);

- return 1;

-}

-

-void

-rng_register (lua_State *L)

-{

- luaL_newmetatable (L, GS_METATABLE(GS_RNG));

- lua_pushvalue (L, -1);

- lua_setfield (L, -2, "__index");

- luaL_register (L, NULL, rng_methods);

- lua_pop (L, 1);

-

- /* gsl module registration */

- luaL_register (L, "rng", rng_functions);

- lua_pop(L, 1);

-}

diff --git a/lua-rng.h b/lua-rng.h
deleted file mode 100644
index 7ed53f0d..00000000
--- a/lua-rng.h
+++ /dev/null

@@ -1,14 +0,0 @@

-#ifndef LUA_RNG_H

-#define LUA_RNG_H

-

-#include <lua.h>

-#include <gsl/gsl_rng.h>

-

-extern void rng_register (lua_State *L);

-

-struct lua_rng {

- gsl_rng *rng;

- unsigned long int min, max;

-};

-

-#endif

diff --git a/randist-source.c b/randist-source.c
deleted file mode 100644
index c9ed51a3..00000000
--- a/randist-source.c
+++ /dev/null

@@ -1,102 +0,0 @@

-

-typedef double (*my_gsl_func_t)(double, double);

-typedef double (*my_gsl_func_2p_t)(double, double, double);

-

-#define MY_FUNC(name) CONCAT2 (my, DECLINE (name))

-#define MY_DECLARE(name) static int MY_FUNC (name) (lua_State *L);

-#define REG_MY_DECLARE(name) {#name, MY_FUNC (name)},

-

-#define MY_FUNC_IMPLEMENT(name) \

-int MY_FUNC (name) (lua_State *L) \

- { \

- return DECLINE (gener_raw) (L, MY_GSL_FUNC (name)); \

- }

-

-#define MY_FUNC_IMPLEMENT_2P(name) \

-int MY_FUNC (name) (lua_State *L) \

- { \

- return DECLINE (gener_raw_2p) (L, MY_GSL_FUNC (name)); \

- }

-

-#define EXPAND(n) MY_DECLARE(n)

-#define EXPAND_2P(n) MY_DECLARE(n)

-#define EXPAND_OTHER(n) MY_DECLARE(n)

-#include "distributions-list.c"

-#undef EXPAND

-#undef EXPAND_2P

-#undef EXPAND_OTHER

-

-

-static const struct luaL_Reg DECLINE (functions)[] = {

-#define EXPAND(n) REG_MY_DECLARE(n)

-#define EXPAND_2P(n) REG_MY_DECLARE(n)

-#define EXPAND_OTHER(n) REG_MY_DECLARE(n)

-#include "distributions-list.c"

-#undef EXPAND

-#undef EXPAND_2P

-#undef EXPAND_OTHER

- {NULL, NULL}

-};

-

-static int

-DECLINE (gener_raw) (lua_State *L, my_gsl_func_t func)

-{

- double x = luaL_checknumber (L, 1);

- double param = luaL_optnumber (L, 2, 1.0);

- double v = func (x, param);

- lua_pushnumber (L, v);

- return 1;

-}

-

-static int

-DECLINE (gener_raw_2p) (lua_State *L, my_gsl_func_2p_t func)

-{

- double x = luaL_checknumber (L, 1);

- double p1 = luaL_checknumber (L, 2);

- double p2 = luaL_checknumber (L, 3);

- double v = func (x, p1, p2);

- lua_pushnumber (L, v);

- return 1;

-}

-

-#define EXPAND(n) MY_FUNC_IMPLEMENT(n)

-#define EXPAND_2P(n) MY_FUNC_IMPLEMENT_2P(n)

-#define EXPAND_OTHER(n)

-#include "distributions-list.c"

-#undef EXPAND

-#undef EXPAND_2P

-#undef EXPAND_OTHER

-

-int

-MY_FUNC (binomial) (lua_State *L)

-{

- int n = luaL_checkinteger (L, 1);

- double p1 = luaL_checknumber (L, 2);

- int p2 = luaL_checkinteger (L, 3);

- double v;

-

- if (p2 < 0)

- luaL_error (L, "parameter n cannot be negative for binomial distribution");

-

- v = MY_GSL_FUNC (binomial) (n, p1, (unsigned int) p2);

- lua_pushnumber (L, v);

- return 1;

-}

-

-static int

-MY_FUNC (poisson) (lua_State *L)

-{

- int n = luaL_checkinteger (L, 1);

- double param = luaL_optnumber (L, 2, 1.0);

- double v = MY_GSL_FUNC (poisson) (n, param);

- lua_pushnumber (L, v);

- return 1;

-}

-

-void

-INCLINE (register) (lua_State *L)

-{

- lua_newtable (L);

- luaL_register (L, NULL, DECLINE(functions));

- lua_setfield (L, -2, MODULE_NAME);

-}

diff --git a/randist.c b/randist.c
deleted file mode 100644
index 1b51bb88..00000000
--- a/randist.c
+++ /dev/null

@@ -1,155 +0,0 @@

-

-/* randist.c

- *

- * Copyright (C) 2009 Francesco Abbate

- *

- * This program is free software; you can redistribute it and/or modify

- * it under the terms of the GNU General Public License as published by

- * the Free Software Foundation; either version 3 of the License, or (at

- * your option) any later version.

- *

- * This program is distributed in the hope that it will be useful, but

- * WITHOUT ANY WARRANTY; without even the implied warranty of

- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU

- * General Public License for more details.

- *

- * You should have received a copy of the GNU General Public License

- * along with this program; if not, write to the Free Software

- * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.

- */

-

-#include <lua.h>

-#include <lauxlib.h>

-#include <string.h>

-#include <gsl/gsl_rng.h>

-#include <gsl/gsl_randist.h>

-#include <gsl/gsl_cdf.h>

-

-#include "gs-types.h"

-#include "lua-rng.h"

-#include "randist.h"

-

-typedef double (*gsl_ran_func_t)(const gsl_rng *, double);

-typedef double (*gsl_ran_func_2p_t)(const gsl_rng *, double, double);

-

-#define RAN_IMPLEMENT(name) \

-int randist_ ## name (lua_State *L) \

- { \

- return randist_gener_raw (L, gsl_ran_ ## name); \

- }

-

-#define RAN_IMPLEMENT_2P(name) \

-int randist_ ## name (lua_State *L) \

- { \

- return randist_gener_raw_2p (L, gsl_ran_ ## name); \

- }

-

-#define RAN_DECLARE(name) static int randist_ ## name (lua_State *L)

-#define LUAREG_DECLARE(name) {#name, randist_ ## name}

-

-RAN_DECLARE (gaussian);

-RAN_DECLARE (exponential);

-RAN_DECLARE (chisq);

-RAN_DECLARE (laplace);

-RAN_DECLARE (tdist);

-RAN_DECLARE (cauchy);

-RAN_DECLARE (poisson);

-RAN_DECLARE (fdist);

-RAN_DECLARE (gamma);

-RAN_DECLARE (beta);

-RAN_DECLARE (binomial);

-RAN_DECLARE (gaussian_tail);

-RAN_DECLARE (exppow);

-RAN_DECLARE (rayleigh);

-RAN_DECLARE (lognormal);

-RAN_DECLARE (flat);

-

-static const struct luaL_Reg randist_functions[] = {

- LUAREG_DECLARE (gaussian),

- LUAREG_DECLARE (exponential),

- LUAREG_DECLARE (chisq),

- LUAREG_DECLARE (laplace),

- LUAREG_DECLARE (tdist),

- LUAREG_DECLARE (poisson),

- LUAREG_DECLARE (fdist),

- LUAREG_DECLARE (gamma),

- LUAREG_DECLARE (beta),

- LUAREG_DECLARE (binomial),

- LUAREG_DECLARE (gaussian_tail),

- LUAREG_DECLARE (exppow),

- LUAREG_DECLARE (cauchy),

- LUAREG_DECLARE (rayleigh),

- LUAREG_DECLARE (lognormal),

- LUAREG_DECLARE (flat),

- {NULL, NULL}

-};

-

-static int

-randist_gener_raw (lua_State *L, gsl_ran_func_t func)

-{

- struct lua_rng *urng = gs_check_userdata (L, 1, GS_RNG);

- double param = luaL_optnumber (L, 2, 1.0);

- double v = func (urng->rng, param);

- lua_pushnumber (L, v);

- return 1;

-}

-

-static int

-randist_gener_raw_2p (lua_State *L, gsl_ran_func_2p_t func)

-{

- struct lua_rng *urng = gs_check_userdata (L, 1, GS_RNG);

- double p1 = luaL_checknumber (L, 2);

- double p2 = luaL_checknumber (L, 3);

- double v = func (urng->rng, p1, p2);

- lua_pushnumber (L, v);

- return 1;

-}

-

-RAN_IMPLEMENT(gaussian)

-RAN_IMPLEMENT(exponential)

-RAN_IMPLEMENT(chisq)

-RAN_IMPLEMENT(laplace)

-RAN_IMPLEMENT(tdist)

-RAN_IMPLEMENT(cauchy)

-RAN_IMPLEMENT(rayleigh)

-

-RAN_IMPLEMENT_2P(fdist)

-RAN_IMPLEMENT_2P(gamma)

-RAN_IMPLEMENT_2P(beta)

-RAN_IMPLEMENT_2P(gaussian_tail)

-RAN_IMPLEMENT_2P(exppow)

-RAN_IMPLEMENT_2P(lognormal)

-RAN_IMPLEMENT_2P(flat)

-

-int

-randist_binomial (lua_State *L)

-{

- struct lua_rng *urng = gs_check_userdata (L, 1, GS_RNG);

- double p1 = luaL_checknumber (L, 2);

- int p2 = luaL_checkinteger (L, 3);

- double v;

-

- if (p2 < 0)

- luaL_error (L, "parameter n cannot be negative for binomial distribution");

-

- v = gsl_ran_binomial (urng->rng, p1, (unsigned int) p2);

- lua_pushnumber (L, v);

- return 1;

-}

-

-static int

-randist_poisson (lua_State *L)

-{

- struct lua_rng *urng = gs_check_userdata (L, 1, GS_RNG);

- double param = luaL_optnumber (L, 2, 1.0);

- unsigned int v = gsl_ran_poisson (urng->rng, param);

- lua_pushnumber (L, (double) v);

- return 1;

-}

-

-void

-randist_register (lua_State *L)

-{

- luaL_register (L, "rnd", randist_functions);

- lua_pop(L, 1);

-}

diff --git a/randist.h b/randist.h
deleted file mode 100644
index 0f428c45..00000000
--- a/randist.h
+++ /dev/null

@@ -1,9 +0,0 @@

-#ifndef RANDIST_H

-#define RANDIST_H

-

-#include <lua.h>

-#include "defs.h"

-

-extern void randist_register (lua_State *L);

-

-#endif

diff --git a/randist-init.lua b/randist.lua
index 0a1c71d2..0a1c71d2 100644
--- a/randist-init.lua
+++ b/randist.lua

diff --git a/rnd.lua b/rnd.lua
new file mode 100644
index 00000000..8e017ceb
--- /dev/null
+++ b/rnd.lua

@@ -0,0 +1,51 @@

+

+local gsl = require 'gsl'

+local ran = {}

+

+ran.bernoulli = gsl.gsl_ran_bernoulli

+ran.beta = gsl.gsl_ran_beta

+ran.binomial = gsl.gsl_ran_binomial

+ran.binomial_knuth = gsl.gsl_ran_binomial_knuth

+ran.binomial_tpe = gsl.gsl_ran_binomial_tpe

+ran.exponential = gsl.gsl_ran_exponential

+ran.exppow = gsl.gsl_ran_exppow

+ran.cauchy = gsl.gsl_ran_cauchy

+ran.chisq = gsl.gsl_ran_chisq

+ran.erlang = gsl.gsl_ran_erlang

+ran.fdist = gsl.gsl_ran_fdist

+ran.flat = gsl.gsl_ran_flat

+ran.gamma = gsl.gsl_ran_gamma

+ran.gamma_int = gsl.gsl_ran_gamma_int

+ran.gamma_mt = gsl.gsl_ran_gamma_mt

+ran.gamma_knuth = gsl.gsl_ran_gamma_knuth

+ran.gaussian = gsl.gsl_ran_gaussian

+ran.gaussian_ratio_method = gsl.gsl_ran_gaussian_ratio_method

+ran.gaussian_ziggurat = gsl.gsl_ran_gaussian_ziggurat

+ran.ugaussian = gsl.gsl_ran_ugaussian

+ran.ugaussian_ratio_method = gsl.gsl_ran_ugaussian_ratio_method

+ran.gaussian_tail = gsl.gsl_ran_gaussian_tail

+ran.ugaussian_tail = gsl.gsl_ran_ugaussian_tail

+ran.bivariate_gaussian = gsl.gsl_ran_bivariate_gaussian

+ran.landau = gsl.gsl_ran_landau

+ran.geometric = gsl.gsl_ran_geometric

+ran.hypergeometric = gsl.gsl_ran_hypergeometric

+ran.gumbel1 = gsl.gsl_ran_gumbel1

+ran.gumbel2 = gsl.gsl_ran_gumbel2

+ran.logistic = gsl.gsl_ran_logistic

+ran.lognormal = gsl.gsl_ran_lognormal

+ran.logarithmic = gsl.gsl_ran_logarithmic

+ran.pascal = gsl.gsl_ran_pascal

+ran.pareto = gsl.gsl_ran_pareto

+ran.poisson = gsl.gsl_ran_poisson

+ran.rayleigh = gsl.gsl_ran_rayleigh

+ran.rayleigh_tail = gsl.gsl_ran_rayleigh_tail

+ran.tdist = gsl.gsl_ran_tdist

+ran.laplace = gsl.gsl_ran_laplace

+ran.levy = gsl.gsl_ran_levy

+ran.levy_skew = gsl.gsl_ran_levy_skew

+ran.weibull = gsl.gsl_ran_weibull

+

+-- set global variable

+rnd = ran

+

+return ran

diff --git a/rng.lua b/rng.lua
new file mode 100644
index 00000000..7587a1ef
--- /dev/null
+++ b/rng.lua

@@ -0,0 +1,58 @@

+

+local gsl = require 'gsl'

+local ffi = require 'ffi'

+

+local format = string.format

+

+local M = {}

+

+local rng_type = ffi.typeof('gsl_rng')

+

+local rng_mt = {

+ __tostring = function(s)

+ return format("<Random number generator: %p>", s)

+ end,

+

+ __index = {

+ getint = gsl.gsl_rng_uniform_int,

+ get = gsl.gsl_rng_uniform,

+ set = gsl.gsl_rng_set,

+ },

+}

+

+ffi.metatype(rng_type, rng_mt)

+

+local function rng_type_lookup(s)

+ if s then

+ local ts = gsl.gsl_rng_types_setup()

+ while ts[0] ~= nil do

+ local t = ts[0]

+ if ffi.string(t.name) == s then

+ return t

+ end

+ ts = ts+1

+ end

+ error('unknown generator type: ' .. s)

+ else

+ return gsl.gsl_rng_default

+ end

+end

+

+function M.new(s)

+ local T = rng_type_lookup(s)

+ return gsl.gsl_rng_alloc(T)

+end

+

+function M.list()

+ local t = {}

+ local ts = gsl.gsl_rng_types_setup()

+ while ts[0] ~= nil do

+ t[#t+1] = ffi.string(ts[0].name)

+ ts = ts+1

+ end

+ return t

+end

+

+rng = M

+

+return M