# \[ANN\] DPMMSubClusters.jl - Fast, Distributed, Scaleable inference for Dirichlet Process Mixture Models

**URL:** <https://discourse.julialang.org/t/ann-dpmmsubclusters-jl-fast-distributed-scaleable-inference-for-dirichlet-process-mixture-models/27895>\
**Category:** Package Announcements\
**Created:** [August 23, 2019, 2:38pm UTC](https://discourse.julialang.org/t/ann-dpmmsubclusters-jl-fast-distributed-scaleable-inference-for-dirichlet-process-mixture-models/27895 "2019-08-23T14:38:46Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![Dinari](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dinari/32/10951_2.png) [@Dinari](https://discourse.julialang.org/u/Dinari)\
**Post date:** [August 23, 2019, 2:38pm UTC](https://discourse.julialang.org/t/ann-dpmmsubclusters-jl-fast-distributed-scaleable-inference-for-dirichlet-process-mixture-models/27895/1 "2019-08-23T14:38:46Z")

</div>

## DPMMSubClusters.jl

Package provides an easy, fast and scalable way to perform inference in Dirichlet Process Mixture Models.

> **[GitHub - dinarior/DPMMSubClusters.jl: Distributed MCMC Inference in Dirichlet Process...](https://github.com/dinarior/DPMMSubClusters.jl)**
>
> Distributed MCMC Inference in Dirichlet Process Mixture Models

Developed from the code of:

[Distributed MCMC Inference in Dirichlet Process Mixture Models Using Julia](https://www.cs.bgu.ac.il/~dinari/papers/dpmm_hpml2019.pdf) by Dinari et al.

Which is based on the algorithm from:

[Parallel Sampling of DP Mixture Models using Sub-Clusters Splits](http://people.csail.mit.edu/jchang7/pubs/publications/chang13_NIPS.pdf) by Chang and Fisher.

The package currently supports Gaussian and Multinomial priors, however adding your own is very easy, and more will come in future releases.

 ![dpmm](https://global.discourse-cdn.com/julialang/original/3X/0/5/05d4fcb6c77f07bc704c8ae99115d0d9023270a3.gif)

The package is faster than any other available open code (as far as I know), including solutions for DPMM in Matlab, Python and Julia, even when running only on 1 process.

Few examples:  
[2d Gaussian with ploting](https://nbviewer.jupyter.org/github/dinarior/DPMMSubClusters.jl/blob/master/examples/2d_gaussian/gaussian_2d.ipynb)  
[Image segmentation](https://nbviewer.jupyter.org/github/dinarior/DPMMSubClusters.jl/blob/master/examples/image_seg/dpgmm-superpixels.ipynb)  
[Running-Saving-Loading-Rerunning model](https://nbviewer.jupyter.org/github/dinarior/DPMMSubClusters.jl/blob/master/examples/save_load_model/save_load_example.ipynb)

### Installation

The package (latest version) has the following dependencies:

- CatViews
- Clustering (`0.13.3`)
- Distributions
- JLD2
- NPZ
- SpecialFunctions
- StatsBase
- LinearAlgebra
- Distributed
- DistributedArrays
- Random

To install, simply:

```julia-auto
] add DPMMSubClusters

```

### Usage:

This package is aimed for distributed parallel computing, while working with no workers is possible. Adding more workers, distributed across different machines, is encouraged for increased performance.

It is recommended to use `BLAS.set_num_threads(1)`. When working with larger datasets increasing the amount of workers will do the trick, `BLAS` multi threading might disturb the multiprocessing, resulting in slower inference.

For all the workers to recognize the package, you must start with `@everywhere using DPMMSubClusters`. If you require to set the seed (using the `seed` kwarg), add `@everywhere using Random` as well.

While being very verstile in the setting and configuration, there are 2 modes which you can work with, either the _Basic_, which will use mostly predefined configuration, and will take the data as an argument, or _Advanced_ use, which allows more configuration, loading data from file, and saving the model, or running from a saved checkpoint.

### Basic

In order to run in the basic mode, use the function:

```julia-auto
labels, clusters, weights = fit(all_data::AbstractArray{Float32,2},local_hyper_params::distribution_hyper_params,α_param::Float32;
        iters::Int64 = 100, init_clusters::Int64 = 1,seed = nothing, verbose = true, save_model = false, burnout = 20, gt = nothing)

```

Or, if opting for the default Gaussian weak prior:

```julia-auto
labels, clusters, weights = fit(all_data::AbstractArray{Float32,2},α_param::Float32;
        iters::Int64 = 100, init_clusters::Int64 = 1,seed = nothing, verbose = true, save_model = false,burnout = 20, gt = nothing)

```

\* note that while we dispatch on `Float32`, other numbers will work as well, and will be cast if needed.

#### Args and Kwargs:

- all\_data - The data, should be `DxN`.
- local\_hyper\_params - The prior you plan to use, can be either Multinomial, or `NIW` (example below on how to create one)
- α\_param - Concetration parameter
- iters - Number of iterations
- seed - Random seed, can also be set seperatly. note that if seting seperatly you must set it on all workers.
- verbose - Printing status on every iteration.
- save\_model - If true, will save a checkpoint every 25 iterations, note that if you opt for saving, I recommend the advanced mode.
- burnout - How many iteration before allowing clusters to split/merge, reducing this number will result in faster inference, but with higher variance between the different runs.
- gt - Ground Truth, if supplied will perform `NMI` and `VI` tests on every iteration.

#### Return values:

`fit` will return the following:

```julia-auto
labels, cluster_params, weights, iteration_time_history, nmi_score_history,likelihood_history, cluster_count_history

```

Note that `weights` does not sum up to `1`, but to `1` minus the weight of the non-instanisated components.

### Advanced

In this mode you are required to supply a params file, example for one is the file `global_params.jl`.  
It includes all the configurable params. Running it is as simple as:

```julia-auto
dp = dp_parallel(model_params::String; verbose = true, save_model = true, burnout = 5, gt = nothing)

```

Will return:

```julia-auto
dp, iteration_time_history , nmi_score_history, liklihood_history, cluster_count_history

```

The returned value `dp` is a data structure:

```julia-auto
mutable struct dp_parallel_sampling
    model_hyperparams::model_hyper_params
    group::local_group
end

```

In which contains the `local_group`, another structure:

```julia-auto
mutable struct local_group
    model_hyperparams::model_hyper_params
    points::AbstractArray{Float64,2}
    labels::AbstractArray{Int64,1}
    labels_subcluster::AbstractArray{Int64,1}
    local_clusters::Vector{local_cluster}
    weights::Vector{Float64}
end

```

Note that for data loading the package use `NPZ` , which utilize python _numpy_ files. Thus the data files must be _pythonic_, and be of the shape `NxD`.

## Additional Functions

Additional function exposed to the user include:

- `run_model_from_checkpoint(file_name)` : Used to restart a saved run, file\_name must point to a valid checkpoint file created during a run of the model. Note that the params files used for running the model initialy must still be available and in the same location, this is true for the data as well.
- `calculate_posterior(model)` : Calculate the posterior of a model, returned from `dp_parallel`.
- `generate_gaussian_data(N::Int64, D::Int64, K::Int64, MixtureVar::Number)`: Randomly generates gaussian data, `N` points, of dimension `D` from `K` clusters, with `MixtureVar` variance between mixture componenets means. return value is `points, labels, cluster_means, cluster_covariance`.
- `generate_mnmm_data(N::Int64, D::Int64, K::Int64, trials::Int64)`: Similar to above, just for multinomial data, the return value is `points, labels, clusters`

## Toy Example

```julia-auto
using DPMMSubClusters

#Generate 10k samples of 2D gaussian data, sampled from 6 random gaussians)
x,y,clusters = generate_gaussian_data(10000,2,6,100.0)

#NIW Hyper Params:
# struct niw_hyperparams <: distribution_hyper_params
# κ::Float32
# m::AbstractArray{Float32}
# ν::Float32
# ψ::AbstractArray{Float32}
# end
hyper_params = DPMMSubClusters.niw_hyperparams(1.0,
           zeros(2),
           5,
           [1 0;0 1])

##Run with hyper params
ret_values= fit(x,hyper_params,10.0, iters = 100)

##Run without hyper params,faster burnout and gt
ret_values= fit(x,10.0, iters = 100,burnout = 10, gt = y)

labels = ret_values[1]

```

### Performance Tips

As mentioned above, it is recommended to use `BLAS.set_num_threads(1)`.

The performance increase is not linear with the processes, and on small data sets of lower dimensions adding more processes might even reduce performance.

On that note - Any optimization contributions are very welcomed 🙂

### Misc

For any questions: dinari@post.bgu.ac.il  
Also available here and on Julia slack.

Contributions, feature requests, suggestion etc.. are welcomed.

If you use this code for your work, please cite the following:

```julia-auto
@inproceedings{Dinari:CCGrid:2019,
  title={Distributed {MCMC} Inference in {Dirichlet} Process Mixture Models Using {Julia}},
  author={Dinari, Or and Angel, Yu and Freifeld, Oren and Fisher III, John W},
  booktitle={International Symposium on Cluster, Cloud and Grid Computing (CCGRID) Workshop on High Performance Machine Learning Workshop},
  year={2019}
}

```