library(probdist) for sampling from probability distributions

[dcnorris]

In advance of completing this work and submitting a PR, I'd like to seek feedback on the approach I've taken to sampling from probability distributions in this code from the DEDUCTION project.

I know other Scyerites have thought long and hard about random number generation and other related issues, and would really like to hear their input.

[bakaq]

I haven't done a lot of probability distribution sampling in my explorations of property based testing in Prolog, but I'll share here a snippet of the PRNG with explicit state that I made together with some basic generators for a lot of basic Prolog things.

Notice that the RNG state is always threaded explicitly, so you can specify a seed just by giving a RNG state that you created.

:- use_module(library(lists)).
:- use_module(library(dcgs)).
:- use_module(library(random)).
:- use_module(library(format)).
:- use_module(library(lambda)).
:- use_module(library(charsio)).
:- use_module(library(debug)).

seed_rng(Seed, Rng) :-
    Rng = xorshift64(Seed).

% Gets the next number of the generator
rng_next(xorshift64(Rng0), xorshift64(Rng), Rng) :-
    Rng1 is (Rng0 xor (Rng0 << 37)) mod 2^64,
    Rng2 is (Rng1 xor (Rng1 >> 11)) mod 2^64,
    Rng is (Rng2 xor (Rng2 << 17)) mod 2^64.

% Generates an integer between [Lower, Upper)
rng_randint(Rng0, Rng, Lower, Upper, RandInt) :-
    rng_next(Rng0, Rng, N),
    RandInt is Lower + (N mod max(1, Upper - Lower)).

% Generates a float in the interval [0, 1)
rng_uniform(Rng0, Rng, Uniform) :-
    Max is 2^64 - 1,
    rng_randint(Rng0, Rng, 0, Max, N),
    Uniform is N / Max.

% Generates a float in an normal distribution
rng_normal(Rng0, Rng, Normal) :-
    rng_uniform(Rng0, Rng1, U1),
    rng_uniform(Rng1, Rng, U2),
    Normal is sqrt(-2 * log(U1)) * cos(2 * pi * U2).

% Gets a random element from a list
rng_choice(Rng0, Rng, Alternatives, Choice) :-
    length(Alternatives, N),
    rng_randint(Rng0, Rng, 0, N, Idx),
    nth0(Idx, Alternatives, Choice).

% Generates an arbitrary integer
arbitrary_integer(Rng0, Rng, Int) :-
    rng_normal(Rng0, Rng, X),
    Int is round(X * 1000).

% Generates an arbitrary char
arbitrary_char(Rng0, Rng, Char) :-
    SurrogatesStart is 0xD800,
    SurrogatesEnd is 0xDFFF,
    MaxUnicode is 0x11000,
    SurrogatesSize is SurrogatesEnd - SurrogatesStart,
    Upper is MaxUnicode - SurrogatesSize,
    rng_randint(Rng0, Rng, 0, Upper, Code0),
    (   Code0 < SurrogatesStart ->
        Code = Code0
    ;   Code is Code0 + SurrogatesSize + 1
    ),
    char_code(Char, Code).
    
% Generates an arbitrary list of chars
arbitrary_chars(Rng0, Rng, Chars) :-
    rng_normal(Rng0, Rng1, Len0),
    Len is abs(round(Len0 * 10)),
    length(Chars, Len),
    foldl(\C^R0^R^arbitrary_char(R0, R, C), Chars, Rng1, Rng).

% Generates an arbitrary atom
arbitrary_atom(Rng0, Rng, Atom) :-
    arbitrary_chars(Rng0, Rng, Chars),
    atom_chars(Atom, Chars).

% Generates an arbitrary float
arbitrary_float(Rng0, Rng, Float) :-
    rng_uniform(Rng0, Rng, X),
    Float is tan((X + 0.5) * pi).

% Generates and arbitrary term
arbitrary_term(Rng0, Rng, Term) :-
    rng_randint(Rng0, Rng1, 10, 20, Budget),
    rng_randint(Rng1, Rng2, 1, 5, VarBudget),
    arbitrary_term(Rng2, Rng3, Term, Budget),
    connect_vars(Rng3, Rng, Term, VarBudget).

connect_vars(Rng0, Rng, Term, VarBudget) :-
    term_variables(Term, Variables),
    length(Variables, NVariables),
    (   NVariables =< VarBudget ->
        Rng = Rng0
    ;   rng_choice(Rng0, Rng1, Variables, Var1),
        rng_choice(Rng1, Rng2, Variables, Var2),
        Var1 = Var2,
        connect_vars(Rng2, Rng, Term, VarBudget)
    ).

reduce_budget(Rng0, Rng, Budget0, Budget) :-
    rng_randint(Rng0, Rng, 0, Budget0, Budget1),
    Budget is max(0, Budget1 - 1).

arbitrary_term(Rng0, Rng, Term, Budget) :-
    NoBudgetVariants = [variable, integer, float, atom, chars, char],
    BudgetVariants = [list, compound],
    append(NoBudgetVariants, BudgetVariants, AllVariants),
    (   Budget =:= 0 ->
        rng_choice(Rng0, Rng1, NoBudgetVariants, Variant)
    ;   rng_choice(Rng0, Rng1, AllVariants, Variant)
    ),
    reduce_budget(Rng1, Rng2, Budget, Budget1),
    arbitrary_term(Variant, Rng2, Rng, Term, Budget1).

arbitrary_term(variable, Rng, Rng, _, _).
arbitrary_term(integer, Rng0, Rng, Integer, _) :-
    arbitrary_integer(Rng0, Rng, Integer).
arbitrary_term(float, Rng0, Rng, Float, _) :-
    arbitrary_float(Rng0, Rng, Float).
arbitrary_term(atom, Rng0, Rng, Atom, _) :-
    arbitrary_atom(Rng0, Rng, Atom).
arbitrary_term(chars, Rng0, Rng, Chars, _) :-
    arbitrary_chars(Rng0, Rng, Chars).
arbitrary_term(char, Rng0, Rng, Char, _) :-
    arbitrary_char(Rng0, Rng, Char).

arbitrary_term(list, Rng0, Rng, List, Budget) :-
    arbitrary_list(Rng0, Rng, List, Budget).
arbitrary_term(compound, Rng0, Rng, Compound, Budget) :-
    arbitrary_compound(Rng0, Rng, Compound, Budget).

arbitrary_list(Rng0, Rng, List, Budget) :-
    rng_normal(Rng0, Rng1, Len0),
    Len is abs(round(Len0 * 5)),
    length(List, Len),
    reduce_budget(Rng1, Rng2, Budget, Budget1),
    foldl(Budget1+\L^R0^R^arbitrary_term(R0, R, L, Budget1), List, Rng2, Rng).

arbitrary_compound(Rng0, Rng, Compound, Budget) :-
    arbitrary_atom(Rng0, Rng1, Functor),
    arbitrary_list(Rng1, Rng, Args, Budget),
    Compound =.. [Functor|Args].

?- seed_rng(15, Rng0), arbitrary_term(Rng0, _, Term).
   Rng0 = xorshift64(15), Term = '㋽⳨'(童㐈ﳃ죀ᬬ햴딿㗘崶("⹞瀇἗䲳縇@흷䃪ҭ獗",'༎鸇ߕꗒᯣ谦졪⏃ب鈴鈁',-593.7140112390217,902,뙸촒첩醑闐괃,"璖ၑ爼圵園ṍ鳮푕"),'䪓礒ꮌ칧𐾬쿃𐞪'(-731,_A,볳,3.296122921935652,'ꐇ',귨,304),-1.2433159344266105,苺,"쾜陪𐕁",0.13667606339998004,谎鬡,'𐌦䏸⧅쟜뒈胩碒턌鷐还Ὄ⼱治軹꩔'(졠,_B)).

(There's a lot of chinese characters in atoms, functors and strings because I haven't weighted the codepoint distribution in anyway and chinese characters simply cover a lot of the Unicode range)

[dcnorris]

If this kind of RNG state threading is generally desirable, might this make DCGs the 'natural' formulation for sampling?

[bakaq]

Possibly. There are many things that I think are more naturally expressed threaded like this, but any kind of state (like PRNG state) is a good candidate.

[dcnorris]

Are you thinking monadically, then?

[bakaq]

Kinda, although I think linear types may be a better intuition for this in Prolog.