Randomly generate typable SK-combinator calculus terms using Boltzmann samplers

sk_boltzmann.pl

:- ensure_loaded(library(apply)).
:- ensure_loaded(library(random)).

:- set_prolog_flag(optimise, true).
:- set_prolog_flag(optimise_unify, true).

% Randomly generate typable SK-combinator calculus terms using Boltzmann
% samplers.
%
% Based on the paper Random generation of closed simply-typed λ-terms: a
% synergy between logic programming and Boltzmann samplers
%
% https://arxiv.org/pdf/1612.07682.pdf
%
% Related:
% https://github.com/fredokun/arbogen
% https://github.com/maciej-bendkowski/boltzmann-brain
% https://github.com/maciej-bendkowski/lambda-sampler
%
% I haven't read it yet, but I think Paul Tarau does (among other things) what
% I do here in this paper: `https://arxiv.org/pdf/1910.01775.pdf`.

% The generating function for SK-combinator calculus terms, given the
% combinatorial class defined by `size/2`, is:
%
%           1 - sqrt(1 - 8z)
%    A(z) = ----------------
%                  2z
%
% The radius of convergence of `A(z)` is:
%
%    rho = 0.125
%
% More generally, if there are `n` combinators, the combinatorial class for the
% analogous size function is:
%
%             1 - sqrt(1 - 4nz)
%    A_n(z) = -----------------
%                    2z
%
% The radius of convergence of `A_n(z)` is:
%
%    rho = 1 / 4n
%
% See section 3 of `https://arxiv.org/pdf/1612.07682.pdf`.
size(s, 0).
size(k, 0).
size(a(X, Y), S) :- size(X, S0), size(Y, S1), S is S0 + S1 + 1.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  The Boltzmann sampler  %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

min_size(90).
boltzmann_combinator(R) :- R < 0.5027624308761950.

next(R1, R2, Size1, Size2) :-
  random(R1),
  random(R2),
  Size2 is Size1 + 1.

random_sk(X, A) :-
  random(R),
  random_sk(X, A, R, 0, Size),
  min_size(MinSize),
  Size >= MinSize,
  !.
random_sk(X, A) :- random_sk(X, A).

random_sk(X, A, R) -->
  { boltzmann_combinator(R), !, random(NewR) },
  random_combinator(X, A, NewR).
random_sk(a(X, Y), B, _) -->
  next(R1, R2),
  random_sk(X, A -> B, R1),
  random_sk(Y, A0, R2),
  { unify_with_occurs_check(A0, A) }.

random_combinator(k, A -> _B -> A, R) --> { R < 0.5, ! }.
random_combinator(s, (A -> B -> C) -> (A -> B) -> A -> C, _) --> [].

multithreaded_random_sk(X, A) :-
  prolog_flag(cpu_count, MaxThreads0),
  MaxThreads is MaxThreads0 - 1,
  Goal = random_sk(X, A),
  length(Goals, MaxThreads),
  maplist(=(Goal), Goals),
  first_solution([X, A], Goals, []).

% Pretty print an SK-combinator calculus term.
pretty(X) :-
  numbervars(X, 0, _),
  pretty(X, Xs, []),
  maplist(write, Xs),
  nl.

pretty(s) --> [s].
pretty(k) --> [k].
pretty(a(X, Y)) --> ['('], pretty(X), [' '], pretty(Y), [')'].

% Pretty print a type.
pretty_type(A) :-
  numbervars(A, 0, _),
  pretty_type(A, As, []),
  maplist(write, As),
  nl.

pretty_type(A -> B) -->
  ['('], pretty_type(A), !, [' → '], pretty_type(B), [')'].
pretty_type(A) --> [A].

main :-
  multithreaded_random_sk(X, A),
  pretty(X),
  pretty_type(A).

sk_boltzmann.py

#!/usr/bin/env python3

from decimal import Decimal, getcontext

getcontext().prec = 100

d1 = Decimal(1)
d2 = Decimal(2)
d4 = Decimal(4)
d8 = Decimal(8)

def A0(z):
    #
    #         1 - sqrt(1 - 8z)
    # A0(z) = ----------------
    #                2z
    #
    z = Decimal(z)

    dividend = d1 - (d1 - (d8 * z)).sqrt()

    divisor = d2 * z

    return dividend / divisor

# First derivative of `A0`.
def A1(z):
    #
    #         1 - 4z - sqrt(1 - 8z) 
    # A1(z) = ----------------------
    #         2 * sqrt(1 - 8z) * z^2
    #
    z = Decimal(z)
    aux = (d1 - (d8 * z)).sqrt()

    dividend = d1 - (d4 * z) - aux

    divisor = d2 * aux * (z ** d2)

    return dividend / divisor

# Expected size.
def E(x):
    x = Decimal(x)
    return x * (A1(x) / A0(x))

rho = Decimal('0.125')

def parameter(size):
    size = Decimal(size)

    divisor = d2
    guess = rho / divisor

    while True:
        expected_size = E(guess)

        if abs(size - expected_size) < Decimal('0.0001'):
            return guess

        divisor *= d2
        guess += rho / divisor if size > expected_size else -rho / divisor

def probability_combinator(x):
    x = Decimal(x)
    return 2 / A0(x)

size = 90

x = parameter(size)
print(f'min_size({size}).')

p1 = probability_combinator(x)
print(f'boltzmann_combinator(R) :- R < {p1:.16}.')