tobias-kohn / pyfoppl-2

Framework for First Order Probabilistic Programming Languages

PyFOPPL-2

A Framework for First Order Probabilistic Programming Languages

This is an improved version of the PyFOPPL. It is work in progress.

Usage

Write your FOPPL-model and save in a file, say, my_model.foppl. The file should either be in the root folder of your project, or in a folder called foppl-src or foppl-models. A simple model might look as follows:

(let [x (sample (normal 1.0 5.0))
      y (+ x 1)]
  (observe (normal y 2.0) 7.0)
  y)

You will find various models in the examples-folder in this project.

Once you have written your FOPPL-model, you can import it as a graphical model in your Python code like this:

from foppl import imports

# Import your model here:
from my_model import model

state = model.gen_prior_samples()
log_pdf = model.gen_pdf(state)
print(log_pdf)

You can get a visual representation of the graphical model if you have the Python packages networkx and matplotlib installed (preferably also graphviz).

model.display_graph()

The file example.py shows how you might import a FOPPL mode, print out the generated graphical model, or generate samples from it.

Options

Options.debug: bool:
When set to True, PyFOPPL will print out additional debug information. On the one hand, the output of the vertices will include addition information. On the other hand, when running gen_prior_samples(), or gen_pdf(), respectively, it will print out the performed computations as they happen.

Options.log_file: str:
This lets you specify a possible log-file. When given a filename as a string, the graph will print out the entire model and the generated code for gen_prior_samples() as well as gen_pdf() to the specified logfile.

In order to take effect during the import of any models, the options should be set before the actual import. You can, later on, deactivate the debug-option, if you do not need the runtime-output.

from foppl import Options
Options.debug = True

import my_model  # Import your model here!

Options.debug = False
state = my_model.model.gen_prior_samples()
...

License

This project is released unter the MIT-license. See LICENSE.

Papers

Discontinuous Hamiltonian Monte Carlo for Probabilistic Programs

Contributors

• Tobias Kohn
• Bradley Gram-Hansen
• Yuan Zhou
• Frank Wood
• Honseok Yang

pyfoppl-2's Issues

The following models do not compile

(let [x1 (sample (normal 0 1))
      x2 (sample (categorical [0.1 0.2 0.7]))
      y1 7]
  (if (> x1 0)
    (if (> x2 1)
      (observe (normal x1 1) y1)
      (observe (normal (+ x1 x2) 2) y1))
    (observe (normal x2 1) y1) )
  [x1 x2])

Error produced is:

    File "/Users/bradley/Documents/Projects/Pyfo_project/Pyfo/pyfo/utils/state.py", line 227, in _log_pdf
    return self._gen_logpdf(state)
  File "/Users/bradley/Documents/Projects/Pyfo_project/Pyfo/pyfo/pyfoppl/foppl/foppl_model.py", line 172, in gen_pdf
    node.update_pdf(state)
  File "/Users/bradley/Documents/Projects/Pyfo_project/Pyfo/pyfo/pyfoppl/foppl/graphs.py", line 421, in update_pdf
    log_pdf = distr.log_pdf(state[self.name])
  File "/Users/bradley/Documents/Projects/Pyfo_project/Pyfo/pyfo/distributions/distribution.py", line 185, in log_pdf
    return torch.sum(self.batch_log_pdf(x, *args, **kwargs))
  File "/Users/bradley/Documents/Projects/Pyfo_project/Pyfo/pyfo/distributions/categorical.py", line 120, in batch_log_pdf
    boolean_mask = torch_zeros_like(logits.data).scatter_(-1, x.data.long(), 1)
RuntimeError: Invalid index in scatter at /Users/bradley/pytorch/aten/src/TH/generic/THTensorMath.c:636

The second model is the following:

(let [x1 (sample (normal 0 1))
      x2 (sample (normal 0 1))
      cov [[1 0] [0 1]]
      y [1 1]]
  (if (> x1 0)
    (if (> x2 0)
      (observe (mvn [1 1] cov) y)
      (observe (mvn [1 -1] cov) y))
    (if (> x2 0)
      (observe (mvn [-1 1] cov) y)
      (observe (mvn [-1 -1] cov) y)))
  [x1 x2])

Error produced is:

  File "/Users/bradley/Documents/Projects/Pyfo_project/Pyfo/pyfo/utils/state.py", line 42, in __init__
    self._state_init = cls.gen_prior_samples()
  File "/Users/bradley/Documents/Projects/Pyfo_project/Pyfo/pyfo/pyfoppl/foppl/foppl_model.py", line 158, in gen_prior_samples
    node.update(state)
  File "/Users/bradley/Documents/Projects/Pyfo_project/Pyfo/pyfo/pyfoppl/foppl/graphs.py", line 406, in update
    result = self.evaluate_observation(state)
  File "<string>", line 1, in <lambda>
KeyError: 'data_30004'

NN model

Okay, the first issue!

Let's try the nn model

(def latent-dim 2)

(def hidden-dim 10)

(def output-dim 5)

; (require '[clojure.core.matrix :as mat :refer [mmul add mul div sub]])

(defn append-gaussian [_ v]
  (conj v (sample (normal 0.0 1.0))))

(defn make-latent-vector [_]
  (loop latent-dim [] append-gaussian))

(defn make-hidden-vector [_]
  (loop hidden-dim [] append-gaussian))

(defn make-output-vector [_]
  (loop output-dim [] append-gaussian))

(defn append-latent-vector [_ M]
  (conj M (make-latent-vector)))

(defn append-hidden-vector [_ M]
  (conj M (make-hidden-vector)))

(defn append-output-vector [_ M]
  (conj M (make-output-vector)))

(defn relu [v]
  (mul (mat/ge v 0.0) v))

(defn sigmoid [v]
  (div 1.0 (add 1.0 (mat/exp (sub 0.0 v)))))

(defn append-flip [i v p]
  (conj v (sample (flip (nth p i)))))

(let [z (make-latent-vector)

      ;; first: hidden layer
      W (loop hidden-dim [] append-latent-vector)
      b (make-hidden-vector)
      h (relu (add (mmul W z) b))

      ;; output layer
      V (loop output-dim [] append-hidden-vector)
      c (make-output-vector)]
  (loop output-dim [] append-flip (sigmoid (add (mmul V h) c))))

Implementing bijectors + transformations into the compiler

This is less of an issue right now and just a heads up.

So this has very, very recently just been implemented in pytorch:
Bijectors and Transform distributions in pytorch
and here for additional details:
link

I will have a play around with these later and see how we can implement them in our current framework, with caching and without caching.

Is following extracting the .data from behind the scenes?

From the following gmm model:

(def N 10)

(defn sample-components [_ zs pi]
  (let [z (sample (categorical pi))]
    (conj zs z)))

(defn observe-data [n _ ys zs mus]
  (let [y (get ys n)
        z (get zs n)
        mu (get mus z)]
    (observe (normal mu 1) y)
    nil))

(let [ys      (vector -2.0  -2.5  -1.7  -1.9  -2.2
                      1.5  2.2  3  1.2  2.8)
      pi [0.5 0.5]
      zs  (loop N (vector) sample-components pi)
      mu1 (sample (normal 0 100))   ; std = 10
      mu2 (sample (normal 0 100))
      mus (vector mu1 mu2)]
  (loop N nil observe-data ys zs mus)
  (vector mus zs))

we get the code output to be:

Now generating posterior python code 
# data_30001
state['data_30001'] = [-2.0, -2.5, -1.7, -1.9, -2.2, 1.5, 2.2, 3, 1.2, 2.8]
# x30002
log_pdf += dist.Categorical([0.5, 0.5]).log_pdf(state['x30002'])
# x30003
log_pdf += dist.Categorical([0.5, 0.5]).log_pdf(state['x30003'])
# x30004
log_pdf += dist.Categorical([0.5, 0.5]).log_pdf(state['x30004'])
# x30005
log_pdf += dist.Categorical([0.5, 0.5]).log_pdf(state['x30005'])
# x30006
log_pdf += dist.Categorical([0.5, 0.5]).log_pdf(state['x30006'])
# x30007
log_pdf += dist.Categorical([0.5, 0.5]).log_pdf(state['x30007'])
# x30008
log_pdf += dist.Categorical([0.5, 0.5]).log_pdf(state['x30008'])
# x30009
log_pdf += dist.Categorical([0.5, 0.5]).log_pdf(state['x30009'])
# x30010
log_pdf += dist.Categorical([0.5, 0.5]).log_pdf(state['x30010'])
# x30011
log_pdf += dist.Categorical([0.5, 0.5]).log_pdf(state['x30011'])
# x30012
log_pdf += dist.Normal(0, 100).log_pdf(state['x30012'])
# x30013
log_pdf += dist.Normal(0, 100).log_pdf(state['x30013'])
# y30015
log_pdf += dist.Normal(([state['x30012'], state['x30013']] if not ((state['x30012']-state['x30013'])>=0) else [state['x30013'], state['x30012']])[runtime.index(state['x30002'])], 1).log_pdf(-2.0)
# y30016
log_pdf += dist.Normal(([state['x30012'], state['x30013']] if not ((state['x30012']-state['x30013'])>=0) else [state['x30013'], state['x30012']])[runtime.index(state['x30003'])], 1).log_pdf(-2.5)
# y30017
log_pdf += dist.Normal(([state['x30012'], state['x30013']] if not ((state['x30012']-state['x30013'])>=0) else [state['x30013'], state['x30012']])[runtime.index(state['x30004'])], 1).log_pdf(-1.7)
# y30018
log_pdf += dist.Normal(([state['x30012'], state['x30013']] if not ((state['x30012']-state['x30013'])>=0) else [state['x30013'], state['x30012']])[runtime.index(state['x30005'])], 1).log_pdf(-1.9)
# y30019
log_pdf += dist.Normal(([state['x30012'], state['x30013']] if not ((state['x30012']-state['x30013'])>=0) else [state['x30013'], state['x30012']])[runtime.index(state['x30006'])], 1).log_pdf(-2.2)
# y30020
log_pdf += dist.Normal(([state['x30012'], state['x30013']] if not ((state['x30012']-state['x30013'])>=0) else [state['x30013'], state['x30012']])[runtime.index(state['x30007'])], 1).log_pdf(1.5)
# y30021
log_pdf += dist.Normal(([state['x30012'], state['x30013']] if not ((state['x30012']-state['x30013'])>=0) else [state['x30013'], state['x30012']])[runtime.index(state['x30008'])], 1).log_pdf(2.2)
# y30022
log_pdf += dist.Normal(([state['x30012'], state['x30013']] if not ((state['x30012']-state['x30013'])>=0) else [state['x30013'], state['x30012']])[runtime.index(state['x30009'])], 1).log_pdf(3)
# y30023
log_pdf += dist.Normal(([state['x30012'], state['x30013']] if not ((state['x30012']-state['x30013'])>=0) else [state['x30013'], state['x30012']])[runtime.index(state['x30010'])], 1).log_pdf(1.2)
# y30024
log_pdf += dist.Normal(([state['x30012'], state['x30013']] if not ((state['x30012']-state['x30013'])>=0) else [state['x30013'], state['x30012']])[runtime.index(state['x30011'])], 1).log_pdf(2.8)
...

(state['x...']-state['x...']) are missing the conditional_suffix = .data[0] ,is something happening behind the scenes, because those are definitely Variables and no error is being induced?

An adaptation of gmm_a which breaks symmetry

Let's look a model which breaks the symmetry.

(def N 10)

(defn sample-components [_ zs pi]
  (let [z (sample (categorical pi))]
    (conj zs z)))

(defn get-ordered-mu [mu1 mu2]
  (if (< mu1 mu2)
    (vector mu1 mu2)
    (vector mu2 mu1)))

(defn observe-data [n _ ys zs mus]
  (let [y (get ys n)
        z (get zs n)
        mu (get mus z)]
    (observe (normal mu 1) y)
    nil))

(let [ys      (vector -2.0  -2.5  -1.7  -1.9  -2.2
                      1.5  2.2  3  1.2  2.8)
      pi [0.5 0.5]
      zs  (loop N (vector) sample-components pi)
      mu1 (sample (normal 0 100))   ; std = 10
      mu2 (sample (normal 0 100))
      mus (get-ordered-mu mu1 mu2)]
  (loop N nil observe-data ys zs mus)
  (vector mus zs))

I could compile it correctly but met a problem while running inference. Error message as following

ERROR in y30008:
  log_pdf += dist.Normal([state['x30005'], state['x30006']] if not ((state['x30005']-state['x30006'])>=0) else [state['x30006'], state['x30005']][runtime.index(state['x30002'])], 1).log_pdf(-2.0)

Is that something wrong with slicing or is it easy to adapt it?

All GMM models

Here we list a collection of all gmm models, from easy to hard. Might write loop differently. Will update after talking to Tobias.

gmm_a
RVs: zs, mus

(defn sample-components [_ zs pi]
  (let [z (sample (categorical pi))]
    (conj zs z)))

(defn observe-data [n _ ys zs mus]
  (let [y (get ys n)
        z (get zs n)
        mu (get mus z)]
    (observe (normal mu 1) y)
    nil))

(let [ys      (vector -2.0  -2.5  -1.7  -1.9  -2.2
                      1.5  2.2  3  1.2  2.8)
      pi [0.5 0.5]
      zs  (loop 10 (vector) sample-components pi)
      mus (vector (sample (normal 0 100))   ; std = 10
                  (sample (normal 0 100)))]
  (loop 10 nil observe-data ys zs mus)
  (vector mus zs))

gmm_b
RVs: pi, zs, mus

(defn sample-components [_ zs pi]
  (let [z (sample (categorical pi))]
    (conj zs z)))

(defn observe-data [n _ ys zs mus]
  (let [y (get ys n)
        z (get zs n)
        mu (get mus z)]
    (observe (normal mu 1) y)
    nil))

(let [ys      (vector -2.0  -2.5  -1.7  -1.9  -2.2
                      1.5  2.2  3  1.2  2.8)
      k 2
      pi (sample (dirichlet [1.0 1.0]))
      zs  (loop 10 (vector) sample-components pi)
      mus (vector (sample (normal 0 100))   ; std = 10
                  (sample (normal 0 100)))]
  (loop 10 nil observe-data ys zs mus)
  (vector pi zs mus))

gmm_c
RVs: pi, zs, mus, taus

;; fix pi, unknown prior mean for the center of each cluster

(defn sample-components [_ zs pi]
  (let [z (sample (categorical pi))]
    (conj zs z)))

(defn observe-data [n _ ys zs mus taus]
  (let [y (get ys n)
        z (get zs n)
        mu (get mus z)
        tau (get taus z)]
    (observe (normal mu (/ 1 tau)) y)
    nil))

(let [ys      (vector -2.0  -2.5  -1.7  -1.9  -2.2
                      1.5  2.2  3  1.2  2.8)
      k 2
      pi (sample (dirichlet [1.0 1.0]))
      zs  (loop 10 (vector) sample-components pi)
      mus (vector (sample (normal 0 100))   ; std = 10
                  (sample (normal 0 100)))
      taus (vector (sample (gamma 1 1))
                  (sample (gamma 1 1)))]
  (loop 10 nil observe-data ys zs mus taus)
  (vector pi zs mus taus ))

Getting a string back from dependent conditions

In branching BHMC I have to know which 'x***' corresponds to which 'cond_'. I thought that accessing vertex.dependent_conditions would return the str name, which is what I need, but it actually returns the set from condition in self.conditions. But, you cannot access anything in this set, you can't even convert it a list and extract pieces of information as element 0, is just the whole condition.
Is there another way to extract dependent_conditions as just a string? I have a list of the 'cond_' but at the moment you wouldn't be able to find which variable depends on it.

The reason I require this information, is that when I am looking to see if a discontinuity has been crossed, I need both the x and the cond relating to that x, to extract a boolean from the state when the new, state has been generated. I was trying to create a hash map of {vertex.name: vertex.dependent_condition}, which I would then use within the inference.

How to return two arguments in code_distributions.py

@classmethod
  def binomial(cls, args: list):
      cls.__check_arg_count('dirichlet', 2, args)
      arg = args[1]
      if isinstance(arg, SequenceType):
          return ListType(FloatType, arg.size)
      else:
          return FloatType()

I have had added arg = args[1] to return the size of the list of the probability of success for each trial. But, I also need to return arg0 = arg[0] as this is the number of trials. This int value I would be put into the 'Support_size' attribute.

Problem:

(let [x (sample (binomial  2 [0.3, 0.5))]
  x)

This is a binomial with 2 trials and success prob of first trial is 0.3 and the success of the second trail is 0.5.

Currently I am getting the following errors:

(let [x (sample (binomial  2 [0.3, 0.5))]
  x)

Traceback (most recent call last):
  File "/Users/bradley/Documents/Projects/Pyfo_project/Pyfo/pyfo/unittests/models/binomial/binomial_test.py", line 32, in <module>
    main()
  File "/Users/bradley/Documents/Projects/Pyfo_project/Pyfo/pyfo/unittests/models/binomial/binomial_test.py", line 30, in main
    base_continous_test()
  File "/Users/bradley/Documents/Projects/Pyfo_project/Pyfo/pyfo/unittests/models/binomial/binomial_test.py", line 16, in base_continous_test
    import bin as test
  File "<frozen importlib._bootstrap>", line 971, in _find_and_load
  File "<frozen importlib._bootstrap>", line 955, in _find_and_load_unlocked
  File "<frozen importlib._bootstrap>", line 665, in _load_unlocked
  File "/Users/bradley/Documents/Projects/Pyfo_project/Pyfo/pyfo/pyfoppl/foppl/imports.py", line 28, in exec_module
    compile_module(module, input_text)
  File "/Users/bradley/Documents/Projects/Pyfo_project/Pyfo/pyfo/pyfoppl/foppl/imports.py", line 16, in compile_module
    graph, expr = compile(input_text)
  File "/Users/bradley/Documents/Projects/Pyfo_project/Pyfo/pyfo/pyfoppl/foppl/compilers.py", line 915, in compile
    ast = foppl_parser.parse(source)
  File "/Users/bradley/Documents/Projects/Pyfo_project/Pyfo/pyfo/pyfoppl/foppl/foppl_parser.py", line 359, in parse
    return parser.parse(source)
  File "/Users/bradley/Documents/Projects/Pyfo_project/Pyfo/pyfo/pyfoppl/foppl/foppl_parser.py", line 330, in parse
    return self._parsers[head].parse(form)
  File "/Users/bradley/Documents/Projects/Pyfo_project/Pyfo/pyfo/pyfoppl/foppl/foppl_parser.py", line 201, in parse
    items = [self._parse(f) for f in form.tail]
  File "/Users/bradley/Documents/Projects/Pyfo_project/Pyfo/pyfo/pyfoppl/foppl/foppl_parser.py", line 201, in <listcomp>
    items = [self._parse(f) for f in form.tail]
  File "/Users/bradley/Documents/Projects/Pyfo_project/Pyfo/pyfo/pyfoppl/foppl/foppl_parser.py", line 35, in _parse
    return self.parent.parse(form)
  File "/Users/bradley/Documents/Projects/Pyfo_project/Pyfo/pyfo/pyfoppl/foppl/foppl_parser.py", line 330, in parse
    return self._parsers[head].parse(form)
  File "/Users/bradley/Documents/Projects/Pyfo_project/Pyfo/pyfo/pyfoppl/foppl/foppl_parser.py", line 275, in parse
    _, assignments = self._parse_bindings(form[1])
  File "/Users/bradley/Documents/Projects/Pyfo_project/Pyfo/pyfo/pyfoppl/foppl/foppl_parser.py", line 137, in _parse_bindings
    raise SyntaxError('let bindings must be a vector with an even number of elements')
SyntaxError: let bindings must be a vector with an even number of elements

I have updated all foppl_distribution.py to reflect the recently added distributions, but only in the pyfo code as I know that you don't use it and I didn't want to push any changes until this issue is resolved.