Terminating with uncaught exception

[leoguelman]

I'm getting the following error message when running the code below. Any hints would be appreciated.

import warnings

import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse, transforms
import seaborn as sns; sns.set(style="ticks", color_codes=True)
import arviz as az
import pandas as pd
import numpy as np

import pystan

from jax import lax
import jax.numpy as np
from jax.random import PRNGKey, randint
from jax.scipy.special import expit

import numpyro
from numpyro.diagnostics import effective_sample_size
import numpyro.distributions as dist
from numpyro.infer import MCMC, NUTS, Predictive
numpyro.set_host_device_count(4)


def sim_data(a_bar = 1,
             a_sigma = 1.5,
             b1_bar = -2, 
             b1_sigma = 1,
             b2_bar = 2, 
             b2_sigma = 1,
             rho = -0.6, 
             seed=42): # corr between b1 and b2
    
    segment_obs = np.linspace(10, 200, 20, dtype = int)
    segment_id = np.repeat(np.linspace(1, 20, 20, dtype = int), segment_obs)
    segment_N = len(segment_obs)
    N = len(segment_id)

    a_segment = dist.Normal(a_bar, a_sigma).sample(PRNGKey(seed), (segment_N, ))
    sigmas = np.array([b1_sigma, b2_sigma])  # standard deviations
    Rho = np.array([[1, rho], [rho, 1]])  # correlation matrix
    Mu = np.array([b1_bar, b2_bar])
    # now matrix multiply to get covariance matrix
    Sigma = np.diag(sigmas) @ Rho @ np.diag(sigmas)

    beta_segment = dist.MultivariateNormal(Mu, Sigma).sample(PRNGKey(seed), (segment_N,))

    X = dist.Normal(0, 1).sample(PRNGKey(seed), (N, 2))

    lp = np.repeat(a_segment, segment_obs) + \
                   np.repeat(beta_segment[:,0], segment_obs) * X[:,0] + \
                   np.repeat(beta_segment[:,1], segment_obs) * X[:,1] 

    D = dist.Binomial(1, probs=1/(1+np.exp(-lp))).sample(PRNGKey(seed))

    df = pd.DataFrame(dict(D=D, x1=X[:,0], x2=X[:,1], segment_id = segment_id))
    
    return df, a_segment, beta_segment, Mu, Sigma, segment_N 


df_train, a_segment, beta_segment, Mu, Sigma, segment_N = sim_data(seed=922)

df_test, *_ = sim_data(seed=624)

df_train

dat = dict(segment_id=df_train.segment_id.values, 
           D = df_train.D.values,
           x1 = df_train.x1.values,
           x2 = df_train.x2.values)

model_code_stan = """

data {
    int D[2100];
    vector[2100] x1;
    vector[2100] x2;
    int segment_id[2100];
}

parameters {

    vector[20] a;
    vector[20] b1;
    vector[20] b2;
    real a_bar;
    real b1_bar;
    real b2_bar;
    real<lower=0> a_sigma;
    vector<lower=0>[2] sigma_segment;
    corr_matrix[2] Rho;
}

model {
    vector[2100] p;
    Rho ~ lkj_corr(2);
    sigma_segment ~ exponential(1);
    a_sigma ~ exponential(1);
    b1_bar ~ (-2, 1);
    b2_bar ~ normal(2, 1);
    a_bar ~ normal(1, 1.5);
    {
    vector[2] YY[20];
    vector[2] MU;
    MU = [b1_bar, b2_bar]';
    for (j in 1:20) YY[j] = [b1[j], b2[j]]';
    YY ~ multi_normal(MU, quad_form_diag(Rho, sigma_segment));
    }
    a ~ normal(a_bar, a_sigma);
    for (i in 1:2100) {
        p[i] = a[segment_id[i]] + b1[segment_id[i]] * x1[i] + b2[segment_id[i]] * x2[i];
        p[i] = inv_logit(p[i]);
    }
    D ~ binomial(1,p);
    
    }
}

"""

sm = pystan.StanModel(model_code=model_code_stan)
fit = sm.sampling(data=dat, iter=1000, chains=4)

Error Message

libc++abi.dylib: terminating with uncaught exception of type std::invalid_argument
Abort trap: 6

PyStan Version: 2.19.1.1

Python Version: 3.7.4

Operating System: MacOS

[ahartikainen]

Hard to say.

So error happens at the sampling step?

Make sure that the ints going to sampling are really ints (df["var"].values.astype(int))

[leoguelman]

Thank you for your prompt response. I've change this line to ensure var types are as expected. Still getting the same error.

df = pd.DataFrame(dict(D=D.astype(int), x1=X[:,0], x2=X[:,1], segment_id = segment_id.astype(int))) 

def sim_data(a_bar = 1,
             a_sigma = 1.5,
             b1_bar = -2, 
             b1_sigma = 1,
             b2_bar = 2, 
             b2_sigma = 1,
             rho = -0.6, 
             seed=42): # corr between b1 and b2
    
    segment_obs = np.linspace(10, 200, 20, dtype = int)
    segment_id = np.repeat(np.linspace(1, 20, 20, dtype = int), segment_obs)
    segment_N = len(segment_obs)
    N = len(segment_id)

    a_segment = dist.Normal(a_bar, a_sigma).sample(PRNGKey(seed), (segment_N, ))
    sigmas = np.array([b1_sigma, b2_sigma])  # standard deviations
    Rho = np.array([[1, rho], [rho, 1]])  # correlation matrix
    Mu = np.array([b1_bar, b2_bar])
    # now matrix multiply to get covariance matrix
    Sigma = np.diag(sigmas) @ Rho @ np.diag(sigmas)

    beta_segment = dist.MultivariateNormal(Mu, Sigma).sample(PRNGKey(seed), (segment_N,))

    X = dist.Normal(0, 1).sample(PRNGKey(seed), (N, 2))

    lp = np.repeat(a_segment, segment_obs) + \
                   np.repeat(beta_segment[:,0], segment_obs) * X[:,0] + \
                   np.repeat(beta_segment[:,1], segment_obs) * X[:,1] 

    D = dist.Binomial(1, probs=1/(1+np.exp(-lp))).sample(PRNGKey(seed))

    df = pd.DataFrame(dict(D=D.astype(int), x1=X[:,0], x2=X[:,1], segment_id = segment_id.astype(int))) 
    
    return df, a_segment, beta_segment, Mu, Sigma, segment_N 


df_train, a_segment, beta_segment, Mu, Sigma, segment_N = sim_data(seed=922)

df_test, *_ = sim_data(seed=624)

df_train



dat = dict(segment_id=df_train.segment_id.values, 
           D = df_train.D.values,
           x1 = df_train.x1.values,
           x2 = df_train.x2.values)

model_code_stan = """

data {
    int D[2100];
    vector[2100] x1;
    vector[2100] x2;
    int segment_id[2100];
}

parameters {

    vector[20] a;
    vector[20] b1;
    vector[20] b2;
    real a_bar;
    real b1_bar;
    real b2_bar;
    real<lower=0> a_sigma;
    vector<lower=0>[2] sigma_segment;
    corr_matrix[2] Rho;
}

model {
    vector[2100] p;
    Rho ~ lkj_corr(2);
    sigma_segment ~ exponential(1);
    a_sigma ~ exponential(1);
    b1_bar ~ (-2, 1);
    b2_bar ~ normal(2, 1);
    a_bar ~ normal(1, 1.5);
    {
    vector[2] YY[20];
    vector[2] MU;
    MU = [b1_bar, b2_bar]';
    for (j in 1:20) YY[j] = [b1[j], b2[j]]';
    YY ~ multi_normal(MU, quad_form_diag(Rho, sigma_segment));
    }
    a ~ normal(a_bar, a_sigma);
    for (i in 1:2100) {
        p[i] = a[segment_id[i]] + b1[segment_id[i]] * x1[i] + b2[segment_id[i]] * x2[i];
        p[i] = inv_logit(p[i]);
    }
    D ~ binomial(1,p);
    
    }
}

"""

sm = pystan.StanModel(model_code=model_code_stan)
fit = sm.sampling(data=dat, iter=1000, chains=4)

Var Types

{'segment_id': array([ 1, 1, 1, ..., 20, 20, 20], dtype=int32),
'D': array([0, 0, 0, ..., 0, 1, 0], dtype=int32),
'x1': array([-1.4905814 , 0.36862397, 0.3050782 , ..., -0.99298227,
-0.8500902 , -0.69721663], dtype=float32),
'x2': array([-0.0581735, -1.9590635, -1.340289 , ..., -1.4912671, 1.3437554,
-1.0460204], dtype=float32)}

[ahartikainen]

Maybe you need to use np.int64 and np.float64

[farr]

Is it possible that you are encountering this matplotlib bug matplotlib/matplotlib#15410 ? If so, the solution is to do anything requiring multiprocessing before you make any plots using matplotlib (the font cache file that matplotlib keeps open is interacting poorly with the fork call that is used to implement multiprocessing on the latest Mac OS X).

[tacaswell]

The better solution is given in https://bugs.python.org/issue33725#msg329923 which is to set

mp.set_start_method('forkserver')    # or spawn

'spawn' is now the default in py38.

[ahartikainen]

Thanks for the solution.

cc @AlexAndorra

[AlexAndorra]

Thanks for the ping @ahartikainen, and for the proposed solution @tacaswell -- this has been torturing me for the past couple of days! Gonna try and see if this fixes things on PyMC3 🤞

[farr]

Oh, score! This is great.

…

On Thu, May 7, 2020 at 6:01 AM Alexandre ANDORRA ***@***.***> wrote: Thanks for the ping @ahartikainen <https://github.com/ahartikainen>, and for the proposed solution @tacaswell <https://github.com/tacaswell> -- this has been torturing me for the past couple of days! Gonna try and see if this fixes things on PyMC3 🤞 — You are receiving this because you commented. Reply to this email directly, view it on GitHub <https://github.com/stan-dev/pystan/issues/691#issuecomment-625155880>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAAV7QD5Y6IPOP3FUCLJVRDRQKBF7ANCNFSM4LSTOK3A> .

[AlexAndorra]

@ahartikainen, in case this is useful, here's how we fixed the issue on PyMC3