.. _examples:

Examples
========

Example 1: Basics
-----------------

The source file can be downloaded here:
:download:`example1_basics.py <_static/example1_basics.py>`

This example illustrates the basic operations of the StrainPycon package. First
we draw random strain barcodes and frequencies, and then we generate two
measurement vectors: one without noise and one with small additive Gaussian
noise. The reconstruction, i.e., "MAP estimate", is computed for both vectors.
We also compute the misfit values, i.e., negative log-likelihoods, when the
number of strains in the reconstruction varies.

We start by importing the necessary modules and setting the random seed for
reproducibility:

.. literalinclude:: _static/example1_basics.py
    :lines: 1-4

We choose the number of SNP sites and the number of strains for our synthetic
measurements:

.. literalinclude:: _static/example1_basics.py
    :lines: 6-7

The binary strain barcode matrix has shape ``(n, m)``. We draw it from uniform
distribution. The frequency vector has shape ``(n,)`` and it has positive
entries that sum up to one. StrainPycon provides a utility function for drawing
uniformly random frequency vectors. For illustration purposes, we sort the
vector in descending order, because the reconstruction method always returns the
vector in descending order. The noiseless measurement is obtained by simply
taking a weighted sum of the binary barcodes, which can be written as a
vector-matrix multiplication. For the noisy measurement, we add independent
Gaussian random variables with standard deviation 0.01:
    
.. literalinclude:: _static/example1_basics.py
    :lines: 9-12
    
Now we are ready to test the strain identification functions in StrainPycon
package. Most of the useful functions are methods of StrainRecon class, so first
we create an object of that class:

.. literalinclude:: _static/example1_basics.py
    :lines: 14

The MAP estimates are computed as follows. Here, the correct number of strains
is assumed to be known:

.. literalinclude:: _static/example1_basics.py
    :lines: 16-17

Finally, we examine the misfits, i.e., negative log-likelihoods, when ``n`` in
the reconstruction varies between 1 and 7:

.. literalinclude:: _static/example1_basics.py
    :lines: 19-21

The matrices, frequency vectors, and misfits can be visualized as follows:

.. literalinclude:: _static/example1_basics.py
    :lines: 23-


Example 2: Uncertainty quantification
-------------------------------------

The source file can be downloaded here:
:download:`example2_uncertainty.py <_static/example2_uncertainty.py>`

In this example, we consider experimental data and visualize the posterior
statistics for the barcode matrix. Essentially, we reproduce the matrix part of
the image shown as Figure 3 in https://doi.org/10.1088/1361-6420/aad7cd . See
the corresponding section in the text for details about the measurement data.

We import the modules as above and provide the data. The number of SNP sites is
16 and there are three strains:

.. literalinclude:: _static/example2_uncertainty.py
    :lines: 1-27

Next, we compute the reconstruction with uncertainty information. We assume the
knowledge of the number of strains, and the standard deviation of noise is set
to 0.1 as in the article:

.. literalinclude:: _static/example2_uncertainty.py
    :lines: 29-30

Lastly, we create the figure:

.. literalinclude:: _static/example2_uncertainty.py
    :lines: 32-
    

Example 3: Multiple categories
------------------------------

The source file can be downloaded here:
:download:`example3_multicat.py <_static/example3_multicat.py>`

Here we look at the case where the SNP sites are not binary but can have more
than two (here, four) possible values:

.. literalinclude:: _static/example3_multicat.py
    :lines: 1-8

For convenience, the package provides a function that we can use to generate
synthetic data. The last argument (here, 0.1) is the standard deviation of the
Gaussian noise:

.. literalinclude:: _static/example3_multicat.py
    :lines: 10-11

The shape of ``meas`` is now ``(cats, m)`` instead of ``(m,)``. However, we can
use the ``compute`` method as before:

.. literalinclude:: _static/example3_multicat.py
    :lines: 13
    
Visualization:

.. literalinclude:: _static/example3_multicat.py
    :lines: 15-