Error using Ordered() in Mixture of Uniform Distributions

[bdyetton]

@aseyboldt

class Ordered(pm.distributions.transforms.ElemwiseTransform): name = “ordered” def forward(self, x): out = tt.zeros(x.shape) out = tt.inc_subtensor(out[0], x[0]) out = tt.inc_subtensor(out[1:], tt.log(x[1:] - x[:-1])) return out def forward_val(self, x, point=None): x, = pm.distributions.distribution.draw_values(, point=point) return self.forward(x) def backward(self, y): out = tt.zeros(y.shape) out = tt.inc_subtensor(out[0], y[0]) out = tt.inc_subtensor(out[1:], tt.exp(y[1:])) return tt.cumsum(out) def jacobian_det(self, y): return tt.sum(y[1:])

Hi @aseyboldt and @junpenglao,

I’m trying to use this Ordered() class to order a Mixture of Uniform Distributions.

I not sure its relevent, but the application is to find the true onset of an “spindle” event (0-25 seconds) from the estimates of a number of noisy human raters, each with differing expertise. I’m also trying to infer there expertise by considering how close each rater is to the inferred true onset.

My code is as follows:

    max_spindles_per_epoch = 6
    expected_std_for_accuracy = 0.2
    n_r = data['rater_i'].max()+1

    #set up model
    with pm.Model() as model:
        BoundedNormal = pm.Bound(pm.Normal, lower=0)
        p = pm.Dirichlet('p', a=np.ones(max_spindles_per_epoch)) #proportion of rater spindles to each real spindle. Uniform.
        gss = pm.Uniform('gss', lower=0, upper=25, shape=max_spindles_per_epoch, transform=Ordered(),
                         testval=np.linspace(0, 25, max_spindles_per_epoch))  # Real spindles
        comp_dists = gss.distribution
        comp_dists.mode = 12.5 #hack, https://discourse.pymc.io/t/how-to-use-a-densitydist-in-a-mixture/1371/2
        s = pm.Mixture('m_s', w=p, comp_dists=comp_dists)
        rater_expertise = BoundedNormal('r_E', mu=expected_std_for_accuracy, sd=0.25, shape=n_r)
        raters_spindle_location = pm.Normal('raters_spindle_location',
                                            mu=s, 
                                            sd=rater_expertise[data['rater_i']],
                                            observed=data['raters_spindle_location'])

    return model

Unfortunately, I’m getting an error with no real message returned:

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [r_E, m_s, gss, p]
Sampling 2 chains:   0%|          | 0/2000 [00:00<?, ?draws/s]forrtl: error (200): program aborting due to control-C event
Image              PC                Routine            Line        Source             
libifcoremd.dll    00007FFF10EA94C4  Unknown               Unknown  Unknown
KERNELBASE.dll     00007FFF668356FD  Unknown               Unknown  Unknown
KERNEL32.DLL       00007FFF6A1A3034  Unknown               Unknown  Unknown
ntdll.dll          00007FFF6A3A1461  Unknown               Unknown  Unknown
forrtl: error (200): program aborting due to control-C event
Image              PC                Routine            Line        Source             
libifcoremd.dll    00007FFF10EA94C4  Unknown               Unknown  Unknown
KERNELBASE.dll     00007FFF668356FD  Unknown               Unknown  Unknown
KERNEL32.DLL       00007FFF6A1A3034  Unknown               Unknown  Unknown
ntdll.dll          00007FFF6A3A1461  Unknown               Unknown  Unknown
forrtl: error (200): program aborting due to control-C event
Image              PC                Routine            Line        Source             
libifcoremd.dll    00007FFF10EA94C4  Unknown               Unknown  Unknown
KERNELBASE.dll     00007FFF668356FD  Unknown               Unknown  Unknown
KERNEL32.DLL       00007FFF6A1A3034  Unknown               Unknown  Unknown
ntdll.dll          00007FFF6A3A1461  Unknown               Unknown  Unknown
Process finished with exit code 0

What can i do to track this error, or have i missspecified my model somehow?
Thanks!!

[junpenglao]

Does your model work (can run sampling) if you remove the Ordered transformation?

[bdyetton]

Hi @junpenglao,
Thanks for getting back to me so quickly.
Yes, it samples fine, although, because of the lack of constraints, the gss variables are not informative (or at least that’s my intuition)

p.s. i also tried using the pm.Potential() to set the order constraints, but this resulted in a bad inital energy error.

[junpenglao]

Could you please post your data and model in a notebook?

[aseyboldt]

Hi
This or similar errors have been reported from Windows users a couple of times, and I have spend quite a bit of time trying to figure out what’s going on there. So far I haven’t even been able to reproduce it. And I don’t really know anything about Windows, which certainly makes debugging a lot harder. It is interesting that it seems to be related to the ordered transformation here, because so far I suspected that it might be related to multithreading in blas. If this is the same issue as before, then maybe that isn’t the case after all. A reproducible notebook would really be great. Also, what Windows version are you using? And what is the output of numpy.__config__.show()?
Could you also try to run this with sample(cores=1) so that it doesn’t use multiprocessing?
Thanks.

[junpenglao]

It might related to indexing in theano under WinOS.

[bdyetton]

Hi Guys,
I ended up moving to the non-marginalized mixture model (here) and enforcing order using pm.potential. This worked as expected

For anyone else interested in mixture models, code is below:

    max_spindles_per_epoch = 6
    n_r = data['rater_i'].max()+1
    expected_std_for_accuracy = 0.2
    n_t = data.shape[0]

    max_spindle_per_rater = data['spindle_i'].max()

    #set up model
    with pm.Model() as model:
        BoundedNormal = pm.Bound(pm.Normal, lower=0.)
        p = pm.Dirichlet('p', a=np.ones(max_spindles_per_epoch),
                         testval=[2]*max_spindles_per_epoch) #proportion of rater spindles to each real spindle.
        gss = pm.Uniform('gss', lower=0., upper=25., shape=max_spindles_per_epoch,
                         testval=[1, 2, 6, 9, 12, 24])  # Real spindles
        gsd = pm.Normal('gsd', mu=duration_av, sd=duration_sd, shape=max_spindles_per_epoch,
                         testval=[1, 1, 1, 1, 1, 1])  # Real spindles

        # Enforce ordering of markers
        switching = tt.switch(gss[1] - gss[0] < 0, -np.inf, 0)
        for i in range(1, max_spindles_per_epoch-1):
            switching += tt.switch(gss[i+1]-gss[i] < 0, -np.inf, 0)

        pm.Potential('order', switching)

        catergory_w = pm.Categorical('w', p=p, shape=n_t)

        rater_expertise = BoundedNormal('r_E', mu=expected_std_for_accuracy, sd=0.25, shape=n_r)

        raters_spindle_location = pm.Normal('s',
                                            mu=gss[catergory_w],
                                            sd=rater_expertise[data['rater_i']],
                                            observed=data['s'])

As for debugging the original error, here is a python notebook:

bdyetton/pymc3_ordered_issues

Interestingly the notebook gave more errors than pycharm did. @aseyboldt looks like you were right, its some sort of multiprocess error. Setting njobs=1 does not throw the multiprocess error, instead, it gives a bad initial energy error (i.e. model problems, not pymc3 bugs).
Given that I have moved to a different model, i dont think its worth your time debugging my old model

[bdyetton]

Oh, and Windows 10 Creator update (lastest).

numpy:

mkl_info:
    libraries = ['mkl_rt']
    library_dirs = ['C:/Users/bdyet/Anaconda3\\Library\\lib']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['C:\\Program Files (x86)\\IntelSWTools\\compilers_and_libraries_2016.4.246\\windows\\mkl', 'C:\\Prog
ram Files (x86)\\IntelSWTools\\compilers_and_libraries_2016.4.246\\windows\\mkl\\include', 'C:\\Program Files (x86)\\Int
elSWTools\\compilers_and_libraries_2016.4.246\\windows\\mkl\\lib', 'C:/Users/bdyet/Anaconda3\\Library\\include']
blas_mkl_info:
    libraries = ['mkl_rt']
    library_dirs = ['C:/Users/bdyet/Anaconda3\\Library\\lib']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['C:\\Program Files (x86)\\IntelSWTools\\compilers_and_libraries_2016.4.246\\windows\\mkl', 'C:\\Prog
ram Files (x86)\\IntelSWTools\\compilers_and_libraries_2016.4.246\\windows\\mkl\\include', 'C:\\Program Files (x86)\\Int
elSWTools\\compilers_and_libraries_2016.4.246\\windows\\mkl\\lib', 'C:/Users/bdyet/Anaconda3\\Library\\include']
blas_opt_info:
    libraries = ['mkl_rt']
    library_dirs = ['C:/Users/bdyet/Anaconda3\\Library\\lib']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['C:\\Program Files (x86)\\IntelSWTools\\compilers_and_libraries_2016.4.246\\windows\\mkl', 'C:\\Prog
ram Files (x86)\\IntelSWTools\\compilers_and_libraries_2016.4.246\\windows\\mkl\\include', 'C:\\Program Files (x86)\\Int
elSWTools\\compilers_and_libraries_2016.4.246\\windows\\mkl\\lib', 'C:/Users/bdyet/Anaconda3\\Library\\include']
lapack_mkl_info:
    libraries = ['mkl_rt']
    library_dirs = ['C:/Users/bdyet/Anaconda3\\Library\\lib']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['C:\\Program Files (x86)\\IntelSWTools\\compilers_and_libraries_2016.4.246\\windows\\mkl', 'C:\\Prog
ram Files (x86)\\IntelSWTools\\compilers_and_libraries_2016.4.246\\windows\\mkl\\include', 'C:\\Program Files (x86)\\Int
elSWTools\\compilers_and_libraries_2016.4.246\\windows\\mkl\\lib', 'C:/Users/bdyet/Anaconda3\\Library\\include']
lapack_opt_info:
    libraries = ['mkl_rt']
    library_dirs = ['C:/Users/bdyet/Anaconda3\\Library\\lib']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['C:\\Program Files (x86)\\IntelSWTools\\compilers_and_libraries_2016.4.246\\windows\\mkl', 'C:\\Prog
ram Files (x86)\\IntelSWTools\\compilers_and_libraries_2016.4.246\\windows\\mkl\\include', 'C:\\Program Files (x86)\\Int
elSWTools\\compilers_and_libraries_2016.4.246\\windows\\mkl\\lib', 'C:/Users/bdyet/Anaconda3\\Library\\include']

[junpenglao]

@bdyetton

instead, it gives a bad initial energy error

You can try setting init="adapt_diag". You cannot use the default "jitter+adapt_diag" here, as the test_value for the ordered RV needs to be ordered (and jitter would likely break that)

[bdyetton]

@junpenglao

init=“adapt_diag”

@junpenglao Brilliant, thanks!!

[leka0024]

Hi @junpenglao @aseyboldt, (feel free to split this off to a new discussion, as needed)

I was trying to use Ordered() on slightly modified version of example here: https://docs.pymc.io/notebooks/marginalized_gaussian_mixture_model.html (really all I did was play around with how many components, the true weights, means, and precisions, and added a lines dict for having the true values in the trace plot)

I tried to use this transform in three different ways, also noting that it now seems to be available directly:

with pm.Model() as model:

    w = pm.Dirichlet('w', np.ones(N))

#     mu = pm.Normal('mu', 0, 10, shape=N, transform=pm.distributions.transforms.Ordered())
#     tau = pm.Gamma('tau', 1, 1, shape=N, transform=pm.distributions.transforms.Ordered())

#     mu = pm.Normal('mu', 0, 10, shape=N, transform=Ordered())
#     tau = pm.Gamma('tau', 1, 1, shape=N, transform=Ordered())

    trafo = Composed(pm.distributions.transforms.LogOdds(), Ordered())
    mu = pm.Normal('mu', 0, 10, shape=N, transform=trafo)
    tau = pm.Gamma('tau', 1, 1, shape=N, transform=trafo)

    x = pm.NormalMixture('x', w, mu, tau=tau, observed=xi)

    trace = pm.sample(init='adapt_diag',tune=3000,random_seed=123)

pm.traceplot(trace,lines=lines_dict)

(note: N was just how many components I used in the simulated mixture data)

However, on all three attempts it gives the same error, always at the first pm.Normal line:

~/venv_py36/lib/python3.6/site-packages/theano/tensor/var.py in __getitem__(self, args)
    508         # Check if the number of dimensions isn't too large.
    509         if self.ndim < index_dim_count:
--> 510             raise IndexError('too many indices for array')
    511 
    512         # Convert an Ellipsis if provided into an appropriate number of

IndexError: too many indices for array

For the bottom two (where I defined the class Ordered() and Composed(), like in the old discussion), the offending line was:

<ipython-input-3-3c5cc1901fbd> in forward(self, x)
      4     def forward(self, x):
      5         out = tt.zeros(x.shape)
----> 6         out = tt.inc_subtensor(out[0], x[0])
      7         out = tt.inc_subtensor(out[1:], tt.log(x[1:] - x[:-1]))
      8         return out

For the built-in Ordered(), the offending line was:

~/venv_py36/lib/python3.6/site-packages/pymc3/distributions/transforms.py in forward(self, x)
    254     def forward(self, x):
    255         y = tt.zeros(x.shape)
--> 256         y = tt.inc_subtensor(y[..., 0], x[..., 0])
    257         y = tt.inc_subtensor(y[..., 1:], tt.log(x[..., 1:] - x[..., :-1]))
    258         return y

Do you have any suggestions?

The code works fine (depending how well separated the data is) without the ordering, just that I noticed often the respective colored lines don’t line up with the color of the trace data, so thought the ordering might be an issue.

Thank you!

[leka0024]

PS - I tried a much simpler thing and it seemed to work pretty well actually.

Just wrapped the mu=pm.Normal(… in a tt.sort() (was also talked about in the old discussion). That seems to have caused some improvement. Interestingly, if I also wrap the tau=pm.Gamma(… in the tt.sort(), then the performance is severely downgrade.

Regardless, the above issue is still at hand because the error occurred in the line for applying transform to mu line.

[junpenglao]

You need to set the testval, something like: testval=np.linspace(0., 1., N). Currently the default value is not handle correctly for ordered transformation.

[leka0024]

Indeed that worked! At it produced probably the best results, again if I only applied it to the means (and not the precisions). Interestingly, the testvals chosen have a dramatic impact. I first set them -1 to 1, since the actual means in my code are symmetric about zero, and that perform VERY poorly and took a long time to run. Setting to 0 to 1, as you said, was about as fast as no sorting, and better results (better than tt.sort() as well). Then tried 0 to 2, results not as good and slower again. Obviously some math and/or code I don’t understand under the hood …

[junpenglao]

Interesting… thanks for reporting back! The testval your provide is the one after transformation (hence they need to be sorted), and the sampler “sees” the one in the unconstrained space. My guess is that different testval might backward map to some large value that becomes bad initialization for NUTS.