Merge pull request #65 from gesinecauer/modifying_pastis_poisson

Modifying pastis poisson

Author: NelleV

README.rst

@@ -14,17 +14,21 @@ For Pastis:
 - pandas
 - iced
 
-Additional dependencies for the diploid version:
+Additional dependencies for new features (diploid inference
+multiscale optimization, etc):
 - python (>= 3.6)
 - autograd (>= 1.3)
 
-All of these dependencies can be installed easily using conda:
+Most of these dependencies can be installed at once using conda:
 `http://conda.pydata.org/miniconda.html <http://conda.pydata.org/miniconda.html>`_
 
 Once conda is installed, just type the following::
 
   conda install numpy scipy scikit-learn pandas
 
+Or, to include the new features::
+
+    conda install numpy scipy scikit-learn pandas autograd
 
 `iced` can be installed via::
 
@@ -34,17 +38,26 @@ Install PASTIS
 --------------
 
 This package uses distutils, which is the default way of installing
-python modules. To install in your home directory, use::
+python modules.
 
-  python setup.py install --user
+To install in your home directory, use::
 
-To install for all users on Unix/Linux::
+    python setup.py install --user
 
-  python setup.py build
-  sudo python setup.py install
+or using pip::
 
+    pip install --user pastis
+
+To install for all users on Unix/Linux::
+
+    python setup.py build
+    sudo python setup.py install
 
 or using pip::
 
-    pip install --user pastis
+  pip install pastis
+
+This will install a python package ``pastis``, and five programs:
+``pastis-mds``, ``pastis-nmds``, ``pastis-pm1``, ``pastis-pm2``, and
+``pastis-poisson``. Calling any of those five programs will display the help.
 

doc/install.rst

@@ -2,13 +2,8 @@
 Installation
 ================================================================================
 
-pastis
-=======
-
-This package uses distutils, which is the default way of installing
-python modules.
-
-The dependencies are:
+Dependencies
+============
 
 - python (>= 2.7)
 - setuptools
@@ -17,17 +12,31 @@ The dependencies are:
 - scikit-learn (>= 0.13)
 - iced
 
-For the diploid version, additional dependencies are required:
+Additional dependencies for new features (diploid inference
+multiscale optimization, etc):
 - python (>= 3.6)
 - autograd (>= 1.3)
 
-Most of these dependencies can be installed at once using `Anaconda
-<http://docs.continuum.io/anaconda/install.html>`_
+Most of these dependencies can be installed at once using conda:
+`http://conda.pydata.org/miniconda.html <http://conda.pydata.org/miniconda.html>`_
+
+Once conda is installed, just type the following::
+
+  conda install numpy scipy scikit-learn pandas
+
+Or, to include the new features::
+
+    conda install numpy scipy scikit-learn pandas autograd
 
 `iced` can be installed using pip::
 
   pip install --user iced
 
+Install PASTIS
+==============
+
+This package uses distutils, which is the default way of installing
+python modules.
 
 To install in your home directory, use::
 
@@ -46,7 +55,7 @@ or using pip::
 
   pip install pastis
 
-This will install a python package ``pastis``, and four programs ``pastis-mds``,
-``pastis-nmds``, ``pastis-pm1`` and ``pastis-pm2``. Calling any of those four
-programs will display the help.
+This will install a python package ``pastis``, and five programs:
+``pastis-mds``, ``pastis-nmds``, ``pastis-pm1``, ``pastis-pm2``, and
+``pastis-poisson``. Calling any of those five programs will display the help.
 

pastis/optimization/callbacks.py

@@ -2,7 +2,7 @@
 from timeit import default_timer as timer
 from datetime import timedelta
 import os
-from .utils_diploid import find_beads_to_remove
+from .utils_poisson import find_beads_to_remove
 from .multiscale_optimization import decrease_lengths_res, decrease_struct_res
 
 

pastis/optimization/constraints.py

@@ -7,7 +7,7 @@
 import autograd.numpy as ag_np
 from autograd.builtins import SequenceBox
 from .multiscale_optimization import decrease_struct_res, decrease_lengths_res
-from .utils_diploid import find_beads_to_remove
+from .utils_poisson import find_beads_to_remove
 
 
 class Constraints(object):
@@ -389,7 +389,7 @@ def distance_between_homologs(structures, lengths, ploidy, mixture_coefs=None,
 
     """
 
-    from .utils_diploid import _format_structures
+    from .utils_poisson import _format_structures
 
     structures = _format_structures(
         structures=structures, lengths=lengths, ploidy=ploidy,

pastis/optimization/counts.py

@@ -8,9 +8,8 @@
 
 from iced.io import write_counts
 
-from .utils_diploid import _constraint_dis_indices
-from .utils_diploid import find_beads_to_remove
-from .utils_diploid import find_beads_to_remove
+from .utils_poisson import _constraint_dis_indices
+from .utils_poisson import find_beads_to_remove
 
 from .multiscale_optimization import decrease_lengths_res
 from .multiscale_optimization import decrease_counts_res
@@ -37,7 +36,6 @@ def ambiguate_counts(counts, lengths, ploidy, exclude_zeros=None):
         Aggregated and ambiguated contact counts matrix.
     """
 
-
     lengths = np.array(lengths)
     n = lengths.sum()
 
@@ -193,7 +191,6 @@ def _check_counts_matrix(counts, lengths, ploidy, exclude_zeros=True,
     """Check counts dimensions, reformat, & excise selected chromosomes.
     """
 
-
     if chrom_subset_index is not None and len(chrom_subset_index) / max(counts.shape) not in (1, 2):
         raise ValueError("chrom_subset_index size (%d) does not fit counts"
                          " shape (%d, %d)." %

pastis/optimization/multiscale_optimization.py

@@ -437,7 +437,7 @@ def get_multiscale_variances_from_struct(structures, lengths, multiscale_factor,
         low-resolution bead.
     """
 
-    from .utils_diploid import _format_structures
+    from .utils_poisson import _format_structures
 
     if multiscale_factor == 1:
         return None

pastis/optimization/pastis_algorithms.py

@@ -45,7 +45,7 @@ def infer(counts_raw, lengths, ploidy, outdir='', alpha=None, seed=0,
           callback_freq=None, callback_function=None, reorienter=None,
           alpha_true=None, struct_true=None, input_weight=None,
           exclude_zeros=False, null=False, mixture_coefs=None, verbose=True):
-    """Infer 3D structures with PASTIS-diploid.
+    """Infer 3D structures with PASTIS via Poisson model.
 
     Optimize 3D structure from Hi-C contact counts data for diploid
     organisms. Optionally perform multiscale optimization during inference.
@@ -416,21 +416,20 @@ def infer(counts_raw, lengths, ploidy, outdir='', alpha=None, seed=0,
         return struct_, infer_var
 
 
-def pastis_diploid(counts, lengths, ploidy, outdir='', chromosomes=None,
-           chrom_subset=None,
-           alpha=None, seed=0, normalize=True, filter_threshold=0.04,
-           alpha_init=-3., max_alpha_loop=20, multiscale_rounds=1,
-           use_multiscale_variance=True, max_iter=10000000000, factr=10000000.,
-           pgtol=1e-05, alpha_factr=1000000000000., bcc_lambda=0.,
-           hsc_lambda=0., hsc_r=None, hsc_min_beads=5, callback_function=None,
-           print_freq=100, history_freq=100, save_freq=None, piecewise=False,
-           piecewise_step=None, piecewise_chrom=None, piecewise_min_beads=5,
-           piecewise_fix_homo=False, piecewise_opt_orient=True,
-           alpha_true=None,
-           struct_true=None, init='msd', input_weight=None,
-           exclude_zeros=False,
-           null=False, mixture_coefs=None, verbose=True):
-    """Infer 3D structures with PASTIS.
+def pastis_poisson(counts, lengths, ploidy, outdir='', chromosomes=None,
+                   chrom_subset=None, alpha=None, seed=0, normalize=True,
+                   filter_threshold=0.04, alpha_init=-3., max_alpha_loop=20,
+                   multiscale_rounds=1, use_multiscale_variance=True,
+                   max_iter=10000000000, factr=10000000., pgtol=1e-05,
+                   alpha_factr=1000000000000., bcc_lambda=0., hsc_lambda=0.,
+                   hsc_r=None, hsc_min_beads=5, callback_function=None,
+                   print_freq=100, history_freq=100, save_freq=None,
+                   piecewise=False, piecewise_step=None, piecewise_chrom=None,
+                   piecewise_min_beads=5, piecewise_fix_homo=False,
+                   piecewise_opt_orient=True, alpha_true=None, struct_true=None,
+                   init='msd', input_weight=None, exclude_zeros=False,
+                   null=False, mixture_coefs=None, verbose=True):
+    """Infer 3D structures with PASTIS via Poisson model.
 
     Infer 3D structure from Hi-C contact counts data for haploid or diploid
     organisms.

pastis/optimization/utils_diploid.py renamed to pastis/optimization/utils_poisson.py

@@ -1,8 +1,6 @@
 from __future__ import print_function
 
 import numpy as np
-from sklearn.metrics import euclidean_distances
-from scipy import sparse
 import sys
 
 
@@ -190,143 +188,3 @@ def _struct_replace_nan(struct, lengths, kind='linear', random_state=None):
             warn('The following chromosomes were all NaN: ' + ' '.join(nan_chroms))
 
         return(interpolated_struct)
-
-
-class ConstantDispersion(object):
-    def __init__(self, coef=7):
-        self.coef = coef
-
-    def fit(self, X, y):
-        return self
-
-    def predict(self, X):
-        return self.coef * np.ones(X.shape)
-
-    def derivate(self, X):
-        return np.zeros(X.shape)
-
-
-def compute_wish_distances(counts, alpha=-3., beta=1., bias=None):
-    """
-    Computes wish distances from a counts matrix
-
-    Parameters
-    ----------
-    counts : ndarray
-        Interaction counts matrix
-
-    alpha : float, optional, default: -3
-        Coefficient of the power law
-
-    beta : float, optional, default: 1
-        Scaling factor
-
-    Returns
-    -------
-    wish_distances
-    """
-    if beta == 0:
-        raise ValueError("beta cannot be equal to 0.")
-    counts = counts.copy()
-    if sparse.issparse(counts):
-        if not sparse.isspmatrix_coo(counts):
-            counts = counts.tocoo()
-        if bias is not None:
-            bias = bias.flatten()
-            counts.data /= bias[counts.row] * bias[counts.col]
-        wish_distances = counts / beta
-        wish_distances.data[wish_distances.data != 0] **= 1. / alpha
-        return wish_distances
-    else:
-        wish_distances = counts.copy() / beta
-        wish_distances[wish_distances != 0] **= 1. / alpha
-
-        return wish_distances
-
-
-def eval_no_f(x, user_data=None):
-    """
-    Evaluate the object function (no objective function).
-    """
-    return 0.
-
-
-def eval_grad_no_f(X, user_data=None):
-    m, n, counts, alpha, beta, d = user_data
-    return np.zeros(n * m)
-
-
-def eval_g_no(x, user_data=None):
-    return np.array([0, 0.])
-
-
-def eval_jac_g_no(x, flag, user_data=None):
-    if flag:
-        return 0, 0
-    return np.array([0., 0.])
-
-
-def eval_g(x, user_data=None):
-    """
-    Computes the constraints
-    """
-
-    m, n, wish_dist, alpha, beta, d = user_data
-
-    x = x.reshape((m, n))
-    dis = euclidean_distances(x)
-    dis = dis ** 2
-    mask = np.invert(np.tri(m, dtype=np.bool))
-    g = np.concatenate([dis[mask].flatten(), ((x - d) ** 2).sum(axis=1)])
-    return g
-
-
-def eval_jac_g(x, flag, user_data=None):
-    """
-    Computes the jacobian for the constraints mentionned in duan-et-al
-    """
-
-    m, n, wish_dist, alpha, beta, d = user_data
-
-    if flag:
-        ncon = m * (m - 1) / 2
-        row = np.arange(ncon).repeat(2 * n)
-
-        tmp = np.arange(n).repeat(m - 1)
-        tmp = tmp.reshape((n, m - 1))
-        tmp = tmp.T
-        tmp1 = np.arange(n, m * n)
-        tmp1 = tmp1.reshape((m - 1, n))
-        tmp = np.concatenate((tmp, tmp1), axis=1)
-        tmp1 = tmp.copy()
-        for it in range(m):
-            tmp += n
-            tmp = tmp[:-1]
-            tmp1 = np.concatenate((tmp1, tmp))
-
-        # The second part of the jacobian is the restrictions on the distances
-        # to the origin and/or the distances to the SPB/nucleolus center
-        row_2 = np.arange(m).repeat(n)
-        col = np.arange(n * m)
-
-        col = np.concatenate([tmp1.flatten(), col])
-        row = np.concatenate([row, row_2])
-        return row.flatten(), col.flatten()
-    else:
-        x = x.reshape((m, n))
-        tmp = x.repeat(m, axis=0).reshape((m, m, n))
-        dif = tmp - tmp.transpose(1, 0, 2)
-        mask = np.invert(np.tri(m, dtype=np.bool))
-        dif = dif[mask]
-        jac = 2 * np.concatenate((dif, - dif), axis=1).flatten()
-
-        # The second part of the jacobian is the restrictions on the distances
-        # to the origin and/or the distances to the SPB/nucleolus center
-        jac2 = 2 * (x - d).flatten()
-        return np.concatenate([jac, jac2]).flatten()
-
-
-def eval_h(x, lagrange, obj_factor, flag, user_data=None):
-    """
-    """
-    return False

pastis/script/pastis-diploid renamed to pastis/script/pastis-poisson

@@ -1,6 +1,6 @@
 #! /usr/bin/env python
 
-from pastis.optimization.pastis_algorithms import pastis_diploid
+from pastis.optimization.pastis_algorithms import pastis_poisson
 import argparse
 
 
@@ -105,4 +105,4 @@ parser.add_argument('--exclude_zeros', default=False,
 parser.add_argument('--null', dest='null', default=False,
                     action='store_true', help=argparse.SUPPRESS)
 
-pastis_diploid(**vars(parser.parse_args()))
+pastis_poisson(**vars(parser.parse_args()))