I tried to use SPPL to rn several BNs from bnrepository https://www.bnlearn.com/bnrepository/, The script I create is attached at https://github.com/shiruizhao/swift/blob/master/example/water_benchmark.py. It reported following error, could you have a look if I use the library wrong?

"python3.8/site-packages/sppl/compilers/sppl_to_python.py", line 360, in compile tree = ast.parse(source_prime) File "python3.8/ast.py", line 47, in parse return compile(source, filename, mode, flags, MemoryError

Consider the following SPML program:

X ~ Normal(0, 1);
if X > 0:
    Z := X**3
else:
    Z := 1 + X

probability(Z > 0);

Here, each sub-SPN in the overall mixture induced by if needs to be able to map Z differently.

https://github.com/probcomp/sum-product-dsl/blob/ba6e29362dec68422924cc23151cad2bce263eea/src/spn.py#L316

Faulty reasoning in disjoint-union algorithm for probabilities.

Key idea: Pr[A or B] // where A and B are in DNF = Pr[A or (B and ~A)] // disjoint union property = Pr[A] + Pr[B and ~A] // since events are disjoint

Since A and B are in DNF, write them as A = ϕ(A; X) and ϕ(A; Y) B = ϕ(B; X) and ϕ(B; Y)

To assess B and ~A, we have: [ϕ(B; X) and ϕ(B; Y)] and ~[ϕ(A; X) and ϕ(B; Y)]

= [ϕ(B; X) and ϕ(B; Y)] and [~ϕ(A; X) or ~ϕ(B; Y)]

= [ϕ(B; X) and ϕ(B; Y) and ~ϕ(A; X)] or [ϕ(B; X) and ϕ(B; Y) and ~ϕ(B; Y)]

= [ϕ(B; X) and ~ϕ(A; X) and ϕ(B; Y)] or [ϕ(B; X) and ϕ(B; Y) and ~ϕ(B; Y)]

the issue is that clauses in this DNF expression not necessarily disjoint.

Example:

ϕ(A; X) = X > 0 ϕ(A; Y) = Y < 1 ϕ(B; X) = X < 1 ϕ(B; Y) = Y < 3

Then the final expression is:

[(X < 1) and ~(X > 0) and (Y < 3)] or [(X < 1) and (Y < 3) and ~(Y < 1)] = [(X < 1) and (X ≤ 0) and (Y < 3)] or [(X < 1) and (Y < 3) and (Y ≥ 1)] = [ (X ≤ 0) and (Y < 3) ] or [ (X < 1) and (1 ≤ Y < 3) ]

which is not disjoint, and reduces to inclusion-exclusion.

X << ['a', 'b']

should render as

X << {'a', 'b'}

so as to distinguish between symbolic entities (x) and string entities ('a', 'b') in the language.

Similarly, for negation we have

~(X << ['a', 'b'])

should render as

X << UniversalSet() \ {'a', 'b'}

The example fairness-income-model.ipynb has strange scale parameters of the two Normal distributions: age and capital_gain. E.g. the first age variable is defined as

age ~= norm(loc=38.4208, scale=184.9151) 

but this might yield very weird values for an age. Here are some random examples:

In [.]: scipy.stats.norm.rvs(loc=38.4208, scale=184.9151,size=10)                                          
Out[.]: 
array([ 456.25033811,   83.38018331,  160.93861454,  129.31420506,
        -81.63928612, -249.89561325, -168.10803064,  103.88993016,
         83.14042404,   14.17379334])
In [.]: np.max(scipy.stats.norm.rvs(loc=38.4208, scale=184.9151,size=10000))                               
Out[.]: 697.0747156159017

Similarly, the first capital_gain is defined as

capital_gain ~= norm(loc=568.4105, scale=sqrt(24248365.5428))

which might yield weird values:

In [.]: scipy.stats.norm.rvs(loc=568.4105, scale=(24248365.5428),size=10)                                  
Out[.]: 
array([  4363679.84553238,  27916806.96928081,  24912266.33616957,
        49814111.74657502, -10570052.00261526, -17910371.43687703,
       -42457176.11145068,  21447098.3397193 , -13593136.75306042,
       -10651368.29633577])

The "decision model" part tests for capital_gain values that is much smaller, in the range of 4500..8300.

I tested to sqrt of the scale values, and then it make a little more sense (the model constrain the age to be > 18 so negative values can be ignored).

# age
In [.]: scipy.stats.norm.rvs(loc=38.4208, scale=math.sqrt(184.9151),size=10)                               
Out[.]: 
array([43.60703122, 30.07099217, 48.63010393, 47.4886019 , 36.10868353,
       34.54669056, -7.37485058, 35.89182812, 31.11277021, 25.77872902])

# capital_gain
In [.]: scipy.stats.norm.rvs(loc=568.4105, scale=math.sqrt(24248365.5428),size=10)                         
Out[.]: 
array([  4058.16881999,   1920.27026625,  -1806.45867856,  -2287.77208551,
        -2710.24830138, -10841.8993238 ,   3112.10957644,  -4933.77496474,
         3259.75181009,  -7738.50462773])

However, after sqrt'ing all the scale values in the model, the result is now quite different from the original encoding. Instead of

female_prior: 0.33066502427854466
female_given_no_hire: 0.33412533774074804
p_female_given_no hire / p_female: 1.046470962495416

the adjusted model give the following:

female_prior: 0.33076164956248716
female_given_no_hire: 0.21962121326454123
p_female_given_no hire / p_female:  -33.60136716120391

However, I'm not sure if this is the proper way of fixing the model. Or perhaps I have misunderstood something here...

In the syntax " x >> dict", we should check whether the keys of the dictionary are numeric or symbolic to determine which constructor to call.

It should suffice to specify the SPN using Python if-else notation and leveraging the ast module (extended docs)

For example, we can define the Indian GPA problem as:

nationality = choose({'India': 0.5, 'USA': 0.5})
score = choose({'imperfect': 0.99, 'perfect': 0.01})

if (nationality == 'India'):
    if (score == 'perfect'):
        gpa = Atomic(10)
    if (score == 'imperfect'):
        gpa = Uniform(0, 10)

if (nationality == 'USA'):
    if (score == 'perfect'):
        gpa = Atomic(4)
    if (score == 'imperfect'):
        gpa = Uniform(0, 4)

In general, any SPN will have this specific structure.

Users report that installing graphviz is a pain. We should re-active the extras in setup.py https://github.com/probcomp/sppl/blob/e77eba5408a749fe4811bb23bf00ff7f912114ea/setup.py#L70

When running spn.prob(GPA <= 0.0) on Indian GPA, error is encountered at this line https://github.com/probcomp/sum-product-dsl/blob/aafed1006b8f28708290911dbf946765139dd495/src/spn.py#L961

Note the use of 0.0 instead of 0 in the conditioning event.

It appears to be an error from Sympy.

>>> sympy.Intersection(sympy.Range(10,11,1), sympy.Interval(-sympy.oo, 0.0))
Traceback (most recent call last):
  File "/home/fsaad/sum-product-dsl/.venv/lib/python3.6/site-packages/sympy/core/compatibility.py", line 410, in as_int
    raise TypeError
TypeError

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/fsaad/sum-product-dsl/.venv/lib/python3.6/site-packages/sympy/sets/fancysets.py", line 426, in __new__
    for w in (start, stop, step)]
  File "/home/fsaad/sum-product-dsl/.venv/lib/python3.6/site-packages/sympy/sets/fancysets.py", line 426, in <listcomp>
    for w in (start, stop, step)]
  File "/home/fsaad/sum-product-dsl/.venv/lib/python3.6/site-packages/sympy/core/compatibility.py", line 416, in as_int
    raise ValueError('%s is not an integer' % (n,))
ValueError: 1.00000000000000 is not an integer

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/fsaad/sum-product-dsl/.venv/lib/python3.6/site-packages/sympy/sets/sets.py", line 1213, in __new__
    return simplify_intersection(args)
  File "/home/fsaad/sum-product-dsl/.venv/lib/python3.6/site-packages/sympy/sets/sets.py", line 1966, in simplify_intersection
    new_set = intersection_sets(s, t)
  File "/home/fsaad/sum-product-dsl/.venv/lib/python3.6/site-packages/sympy/multipledispatch/dispatcher.py", line 198, in __call__
    return func(*args, **kwargs)
  File "/home/fsaad/sum-product-dsl/.venv/lib/python3.6/site-packages/sympy/sets/handlers/intersection.py", line 104, in intersection_sets
    return intersection_sets(a, Range(start, end + 1))
  File "/home/fsaad/sum-product-dsl/.venv/lib/python3.6/site-packages/sympy/sets/fancysets.py", line 431, in __new__
    [0, 1/10, 1/5].'''))
ValueError: 
Finite arguments to Range must be integers; `imageset` can define
other cases, e.g. use `imageset(i, i/10, Range(3))` to give [0, 1/10,
1/5].

spn = (X << poisson(mu=5))
w = spn.condition(X << {1,2,3, 10})

The conditional distribution needs to be a mixture of atoms, instead of assuming full support on the bounds.

Rather than iterate through all 2^m subsets, first compute the probabilities of the m conjunctions and immediately return 0 if each conjunction in the union has probability zero, which removes the duplicate computation inside prune_events.

Is there a good way to execute SPPL model source code while having it rely on a previously existing model? For example, assume I have an SPPL model m (of type ProductSPE) which models variables a and b. I have this source code source in form of a string:

if a < 17 
    c ~= bernoulli(p=0.1)

else:
    c ~= bernoulli(p=0.9)

which relies on a being defined but which would otherwise compile fine if I run:

SPPL_Compiler(source).execute_module().model

Obviously, without knowing about m it won't compile because the compiler doesn't know about a. Is there a good way to tell the compiler about a model (like m) as a starting point? i.e. something like this:

SPPL_Compiler(m, source).execute_module().model

MWE

>>> source = """
i = 0.5
Y ~= bernoulli(p=0 if i == 0.5 else 1)
"""
>>> compiler = SPPL_Compiler(source)
>>> namespace = compiler.execute_module()
TypeError: unsupported operand type(s) for <<: 'float' and 'set'

The reason is that the generated Python is

# MODEL DEFINITION
command = Sequence(
    Sample(Y, bernoulli(p=(0 if (i << {0.5}) else 1))),
)

The offending line is here: https://github.com/probcomp/sppl/blob/8b0fe0c37ed15dd19936d13e0fa652c3b5237cac/src/compilers/sppl_to_python.py#L145

Currently sppl depends on hard-coded strict-equality versions of its dependencies.

This makes its inclusion to ongoing projects somewhat problematic.

We always have the choice of vendoring the code instead of installing it through pip, and ignoring restrictions on the setup.py, but I think it'd be best if sppl used a system similar to poetry or pip-tools, with loosely-restricted versions for use as a library, and strict, locked versions for experiment reproducibility in the papers, etc.

Is this something you have considered?

The .name attribute of a frozen rv_discrete.dist object does not correspond to the underlying class, i.e.,

>>> from scipy.stats import norm
>>> d = norm(loc=0, scale=1)
>>> d.dist.name
'norm'

>>> from scipy.stats import rv_discrete
>>> d = rv_discrete(values=((1, 2), (.5, .5))).freeze()
>>> d.dist.name
'Distribution'

This behavior will cause issue in serializing, since we rely on the name attribute to correspond to a scipy class: https://github.com/probcomp/sppl/blob/efff34fb3d3703247dd7001c36970069c5ac3825/src/compilers/spe_to_dict.py#L48

Related #121

The constructor of rv_discrete does accept a name attribute, which will need to be handled correctly when loading the scipy dist.

>>> from sppl.distributions import rv_discrete
>>> from sppl.spe import DiscreteLeaf
>>> d1 = rv_discrete(values=((1, 2), (.5, .5)))
>>> d2 = rv_discrete(values=((1, 2), (.8, .2)))
>>> l1 = Id("X") >> a
>>> l2 = Id("X") >> b
>>> l1 == l2
True

The problem is that scipy does not store the values under the kwds in the frozen rv_discrete object, that is

>>> l1.dist.kwds
{}
>>> l2.dist.kwds
{}

and so the equality logic checking passes: https://github.com/probcomp/sppl/blob/efff34fb3d3703247dd7001c36970069c5ac3825/src/spe.py#L828-L836

The kwds field is populated by scipy for all other distributions norm, poisson, etc.

Can we bump this version to >=1.7 in requirements (?) I tried a few different things here -- like installing into a fresh conda environment. For whatever reason, this old version of scipy fails to build -- but the newest version is okay.