text:
from typing import List, Callable, NamedTuple, Tuple
import sys

import tensorflow as tf
import tensorflow_probability as tfp
from tensorflow.train import AdamOptimizer
import numpy as np

from rgpy.util import tf_xavier_init
from rgpy.rbms import RBM
from rgpy import samplers

class RGRBM(RBM):
    def __init__(self,
                 energy_psi_of_v: Callable[[List[tf.Tensor]], tf.Tensor],
                 energy_theta_of_ve: Callable[[List[tf.Tensor], List[tf.Tensor]], tf.Tensor],
                 n_visible: int,
                 n_environment: int,
                 n_hidden: int,
                 lmbda=0,
                 regularization="L1",
                 learning_rate=0.01,
                 xavier_const=1.0,
                 use_tqdm=True,
                 # DEPRECATED:
                 tqdm=None,
                 binary=[0, 1]):
        """Initializes an RBM to implement the RSMI algorithm

        This assumes that the energy functions characterizing the (assumed)
        Boltzmann distributions P(V) and P(V,E) are provided.

        Visible bias is 0 since we care only about probability of hiddens
        conditioned on the vnisibles.

        Args:
            energy_theta_of_ve: a function that returns the energy of a state of
                visible and environment spins. Parameterized by Theta-RBM.
                Return shape [?] where ? is # of samples of v, e
            energy_psi_of_v: a function that returns the energy of a state of
                visible spins. Parametrized by Psi-RBM
                Return shape [?] where ? is # of samples of v
        """

        super(RGRBM, self).__init__(
            n_visible,
            n_hidden,
            lmbda=lmbda,
            regularization=regularization,
            learning_rate=learning_rate,
            xavier_const=xavier_const,
            use_tqdm=use_tqdm,
            # DEPRECATED:
            tqdm=tqdm,
            binary=binary
        )
        self.mutual_info = -1.
        self.n_environment = n_environment

        self.v = tf.placeholder(tf.float32, [None, self.n_visible])
        self.e = tf.placeholder(tf.float32, [None, self.n_environment])

        # 2 args: [v, e], well what's the point of the lambda function then?x
        self.energy_theta_of_ve = lambda v, e:(
            energy_theta_of_ve(v, e))

        self.energy_psi_of_v = lambda v:(
            energy_psi_of_v(v))

        # The sampler will perform the internal mc average over samples of
        # energy_theta_lambda_of_vh with clamped h
        #self.internal_mc_sampler=samplers.GenericLatticeKernel(n_spins=self.n_visible)

        # Tensorflow: Set up the computational graph
        self.optimizer = AdamOptimizer(learning_rate=self.learning_rate, epsilon=1.0)
        self.run_update_params = self.update_params(self.v, self.e)
        self._post_init()

    def _post_init(self):
        """ Initializes the tf Variables.
        """
        init = tf.global_variables_initializer()
        self.sess = tf.Session()
        self.sess.run(init)

    def _initialize_vars(self):
        # Through ? refers to the number of elements in a given batch

        # The local field (due to visible spins at each hidden site j)
        # Return shape [?, n_hidden]
        self._local_field_lambda_j=lambda v: (tf.matmul(v, self.W)) + self.hbias

        # The energy at each hidden site j = h_j* local_field at j
        # Return shape [?, n_hidden]
        self._energy_lambda_j=lambda v, h: self._local_field_lambda_j(v) * h

        # The total energy of a configuration of visible & hidden units
        # Return shape [?]
        self.energy_lambda_of_vh=lambda v, h: tf.reduce_sum(self._energy_lambda_j(v, h), axis=1)

        # The energy describing the marginal distribution of v of Boltzmann form
        # Return shape [?]
        self.energy_lambda_of_v=lambda v: (
            tf.reduce_sum(tf.log(1 + tf.exp(self._local_field_lambda_j(v))), axis=1))

        # Return shape [?]
        self.energy_psi_lambda_of_vh=lambda v, h:(
            self.energy_psi_of_v(v)
            + self.energy_lambda_of_vh(v, h)
            - self.energy_lambda_of_v(v))

        # Return shape [?]
        self.energy_theta_lambda_of_veh=lambda v, e, h: (
            self.energy_theta_of_ve(v, e)
            + self.energy_lambda_of_vh(v, h)
            - self.energy_lambda_of_v(v))

        # Return shape [?]
        self.delta_energy_psi_theta=lambda v, e: (
            self.energy_theta_of_ve(v, e) - self.energy_psi_of_v(v))

        # The gradient of energy_lambda_of_vh
        # First element: grads wrt hbias; at index j: h_j
        # Second element: grads wrt W: at index i, j: v_i * h_j
        # We neglect vbias, since this has no influence on P(h|v)
        # Return shape [[?, n_hidden],[?, n_visible, n_hidden]]
        self.delta_of_vh=lambda v, h: (
            [h, (tf.expand_dims(h, 1) * tf.expand_dims(v, 2))])

        # The gradient of energy_psi_lambda_of_vh
        # This is akin to the gradient of energy for P(h|v) of Boltzmann form
        # Return shape [[?, n_hidden],[?, n_visible, n_hidden]]
        def _grad_energy_psi_lambda_of_vh(v, h):
            h_expect = tf.nn.sigmoid(self._local_field_lambda_j(v))
            h_eff = h - h_expect
            return self.delta_of_vh(v, h_eff)

        self.grad_energy_psi_lambda_of_vh=_grad_energy_psi_lambda_of_vh

        # Grad (p_lambda(h|v)) = p_lambda(h|v) * (this term)
        # Return shape [[?, n_hidden],[?, n_visible, n_hidden]]
        def _grad_energy_lambda_of_v_not_h(v, h):
            p_lambda_h_j_given_v=tf.nn.sigmoid(self._energy_lambda_j(v, h))
            h_eff=(1 - p_lambda_h_j_given_v) * h
            return self.delta_of_vh(v, h_eff)

        self.grad_energy_lambda_of_v_not_h=_grad_energy_lambda_of_v_not_h

        self.compute_hidden = self.get_hidden(self.x)

    def gen_internal_mc_v_samples(self,
                                  E: List[tf.Tensor],
                                  h: List[tf.Tensor]):
        """
        Generates internal MC samples of V given E and H(V').

        Args:
            e: state of environmental spins
            h: state of hidden spins

        Returns:
            v_samples List[List[tf.Tensor]]: list of samples of visible spins
        """
        mc_energy_fn=lambda v: tf.reshape(self.energy_psi_lambda_of_vh(
            tf.reshape(v, [1, self.n_visible]),
            tf.reshape(h, [1, self.n_hidden])),
                                          [])

        samples = samplers.generic_generator_graph(
            mc_energy_fn,
            self.n_visible,
            n_burnin_steps=126,
            n_results_per_chain=2,
            n_chains=1,
            n_steps_between_results=126
        )

        # self.internal_mc_sampler.set_energy_fn(mc_energy_fn)

        # init_state = tf.constant(np.random.choice(a=[-1., 1.],
        #                                           size=(1, self.n_visible)),
        #                          dtype=tf.float32)
        # all_samples = tfp.mcmc.sample_chain(
        #     num_results=2,
        #     num_burnin_steps=126,
        #     num_steps_between_results=126,
        #     current_state=init_state,
        #     parallel_iterations=1,
        #     kernel=self.internal_mc_sampler)

        return samples
    def calc_internal_mc_expectation(
            self,
            e: List[tf.Tensor],
            h: List[tf.Tensor],
            expectation_fns: List[Callable[[List[tf.Tensor]], tf.Tensor]]
    ) -> List[tf.Tensor]:

        """
        Generates internal MC samples of V given E and H(V').
        Calculates expectations over these samples from functions in
        expectation_fns

        Args:
            e: state of environmental spins
            h: state of hidden spins

            expectation_fn (`V -> tf.Tensor`): an expectation function
                whose expectation we wish to calculate where `V` is a visible
                state, `E` an environmental and `H` a hidden state.
                `E` and `H` are clamped to the values e, h
        Returns:
            expectation_vals (List[tf.Tensor]): the expectation values
                corresponding to each of the functions in expectation_fns.
                shape = tf.shape(expectation_fns)
        """
        v_samples = self.gen_internal_mc_v_samples(e, h)

        return samplers.expectations(
            v_samples,
            expectation_fns)

    def expect_delta_energy(self,
                            e: List[tf.Tensor],
                            h: List[tf.Tensor]
    ):
        """Expectation value of the delta energy
        with respect to Boltzmann distribution with energy
        energy_psi_lambda_of_vh = E_{\\Psi, \\Lambda}(v, h)

        Args:
            e: state of environmental spins
            h: state of hidden spins

        Returns:
            < \\Delta E_{\\Psi, \\Theta}(v, e) >_h (h fixed)
                <=> \\Delta E_{\\Psi, \\Theta, \\Lambda}(v, e, h)
        """
        return self.calc_internal_mc_expectation(
            e,
            h,
            [lambda _v: self.delta_energy_psi_theta(_v, e)])

    def grad_mutual_info_expectations(self,
                                      e: List[tf.Tensor],
                                      h: List[tf.Tensor],
                                      debug=False):
        """Expectation value of the delta energy
        with respect to Boltzmann distribution with energy
        energy_psi_lambda_of_vh = E_{\\Psi, \\Lambda}(v, h)

        Args:
            e: state of environmental spins
            h: state of hidden spins

        Returns:
            < \\Delta E_{\\Psi, \\Theta}(v, e) >_h (h fixed)
                <=> \\Delta E_{\\Psi, \\Theta, \\Lambda}(v, e, h)
        """
        n_samples=2
        e_resized=tf.reshape(tf.tile(e[0], [n_samples]), [n_samples, self.n_environment])
        h_resized=tf.reshape(tf.tile(h[0], [n_samples]), [n_samples, self.n_hidden])

        _delta_energy_psi_theta=lambda _v: (
            self.delta_energy_psi_theta(_v, e_resized))

        _grad_energy_psi_lambda_of_vh=lambda _v: (
            self.grad_energy_psi_lambda_of_vh(_v, h_resized))

        _prod_delta_and_grad =lambda _v: (
            _grad_energy_psi_lambda_of_vh(_v))
        #            tf.reshape(_delta_energy_psi_theta(_v), [1])
        #           * _grad_energy_psi_lambda_of_vh(_v))

        (
            [delta_energy,
             [grad_h, grad_W],
             [prod_h, prod_W]]
        ) =self.calc_internal_mc_expectation(
            e, h,
            [_delta_energy_psi_theta,
             _grad_energy_psi_lambda_of_vh,
             _prod_delta_and_grad]
        )

        correct_h = lambda _h: tf.reduce_sum(_h, axis=0 )
        correct_W = lambda _W: tf.reduce_sum(_W, axis= 0)

        grad_h = correct_h(grad_h)
        grad_W = correct_W(grad_W)
        prod_h = correct_h(prod_h)
        prod_W = correct_W(prod_W)

        if debug:
            (
                [delta_energy,
                 [grad_h, grad_W],
                 [prod_h, prod_W]]
            ) = (
                [tf.Print(delta_energy, [tf.shape(delta_energy), delta_energy], "Delta E: "),
                 [tf.Print(grad_h, [tf.shape(grad_h), grad_h], "\ngrad h"),
                  tf.Print(grad_W, [tf.shape(grad_W), grad_W], "\ngrad W", summarize=10)],
                 [tf.Print(prod_h, [tf.shape(prod_h), prod_h], "\nprod h"),
                  tf.Print(prod_W, [tf.shape(prod_W), prod_W], "\n prod W", summarize=10)]]
            )
        return [delta_energy,
                [grad_h, grad_W],
                [prod_h, prod_W]]

    def grad_mutual_info_proxy(self,
                               v_prime_batch: List[List[tf.Tensor]],
                               e_batch: List[List[tf.Tensor]],
                               debug=False):
        """
        Approximates the mutual information across samples of X (the entire
        configuration).

        Indices should match, so e_samples[i] is the environmental block
        corresponding to v_samples[i]

        Args:
            v_prime_batch: batch of visible samples from external mc
                Shape: [batch_size, n_visible]
            e_batch: batch of environemtnal samples from external mc
                Shape: [batch_size, n_environment]
            batch_size:

        Returns:
            The gradient of a mutualinformation proxy to be used for
            differentiation in SGD.
            < \\Delta E_{\\Psi, \\Theta}(v, e) >_h (h fixed)
                where h = h(v) from v_samples, e = e from e_samples
                and internal v MC sampled
        """
        batch_size = tf.shape(v_prime_batch)[0]

        # Shape: [batch_size, n_hidden]
        h_batch = self.get_hidden(v_prime_batch)

        # Shape: [[batch_size, n_hidden], [batch_size, n_visible, n_hidden]]
        delta_v_prime_not_h=(
            lambda v, h: self.grad_energy_lambda_of_v_not_h(v, h))

        grad_mutual_info_per_term = (
            lambda delta_cond, delta_energy, prod, delta_not: (
                -delta_energy * (delta_not + delta_cond) + prod))

        def compute_term(a: tf.Tensor):
            """
            Computes the ath term of the external sum in grad_mutual_info_proxy
            TODO: Parallelize this; i.e. allow joint calculation for entire batches
            """
            correct_shape = lambda x, x_len: tf.reshape(x, [-1, x_len])
            v_prime= correct_shape(v_prime_batch[a],self.n_visible)
            e= correct_shape(e_batch[a],self.n_environment)
            h= correct_shape(h_batch[a],self.n_hidden)

            (
                [delta_energy,
                 [h_grad, W_grad],
                 [h_prod, W_prod]]
            ) =self.grad_mutual_info_expectations(e, h)

            [h_grad_not, W_grad_not] = delta_v_prime_not_h(v_prime, h)

            h_term=tf.reshape(
                grad_mutual_info_per_term(h_grad, delta_energy, h_grad, h_grad_not),
                [1])

            W_term=tf.reshape(
                grad_mutual_info_per_term(W_grad, delta_energy, W_prod, W_grad_not),
                [4, 1])
            mi_term = tf.reshape(-delta_energy, [])

            if debug:
                [h_term, W_term, mi_term] = (
                    [tf.Print(h_term, [tf.shape(h_term), h_term], "\nh_term"),
                     tf.Print(W_term, [tf.shape(W_term), W_term], "\nW term", summarize=100),
                     tf.Print(mi_term, [tf.shape(mi_term), mi_term], "\n Mutual info proxy")]
                )

            return [h_term, W_term, mi_term]

        i = tf.constant(0)
        h_sum_init = tf.zeros([self.n_hidden], dtype=tf.float32)
        W_sum_init = tf.zeros([self.n_visible, self.n_hidden], dtype=tf.float32)
        mi_sum_init = tf.constant(0., dtype=tf.float32)

        def body(_i, _h_sum, _W_sum, _mi_sum):
            [_h_update, _W_update, _mi_update] = compute_term(_i)
            return [_i + 1,
                    _h_sum + _h_update,
                    _W_sum + _W_update,
                    _mi_sum + _mi_update]

        cond = lambda _i, _h_sum, _W_sum, _mi_sum: _i < batch_size

        [_, h_sum, W_sum, mi_sum] = tf.while_loop(cond,
                                                  body,
                                                  [i, h_sum_init, W_sum_init, mi_sum_init])

        return [h_sum, W_sum, mi_sum]

    def update_params(self, v, e):
        [*grads, mutual_info] = self.grad_mutual_info_proxy(v, e)
        params = [self.hbias, self.W]
        update_op = self.optimizer.apply_gradients(zip(grads, params), self.global_step)
        return update_op, mutual_info

    def partial_fit(self, batch_ve):
        batch_v, batch_e = batch_ve
        update_op, mutual_info= self.run_update_params
        _, mutual_info_= self.sess.run([update_op, mutual_info],
                                       feed_dict={self.v: batch_v, self.e: batch_e})
        return mutual_info_

    def mutual_info_proxy(self):
        """
        Approximates a mutual information proxy across samples of X (the entire
        configuration).

        Indices should match, so e_samples[i] is the environmental block
        corresponding to v_samples[i]

        Returns:
            An approximation of the mutual information proxy
                < \\Delta E_{\\Psi, \\Theta}(v, e) >_h (h fixed)
                where h = h(v) from v_samples, e = e from e_samples
                and internal v MC sampled
                Shape= []
        """
        h_samples = self.get_hidden(self.v_samples)
        tf_zip = lambda a, b: tf.stack([a, b], 1)
        double_map_fn = lambda map_fn: (lambda ab: map_fn(a[0], a[1]))

        expect_delta_energy = tf.map_fn(
            double_map_fn(self.expect_delta_energy),
            tf_zip(self.e_samples, h_samples))
        return T.mean(expect_delta_energies)

    def fit(self,
            data_ve,
            n_epochs=10,
            batch_size=10,
            print_filters_dir="./",
            shuffle=True,
            verbose=True):
        """
        Trains the Lambda-RBM. Performs SGD where the derivative is given by
        self.grad_mutual_info_proxy.

        Args:
            `data_ve`: training data [ve_samples]
            `n_epochs`: number of epochs to train over
            `batch_size`: number of training samples per batch (external mc-avg)
            `shuffle`: whether to shuffle training data before fitting
            `verbose`: whether to display information about training
        """
        assert n_epochs > 0

        data_ve = self.convert_to_binary(data_ve)

        n_data = min(data_ve.shape[0], 20000)

        if batch_size > 0:
            n_batches = n_data // batch_size + (0 if n_data % batch_size == 0 else 1)
        else:
            n_batches = 1

        batch_spacing = n_data / n_batches

        errs = []

        for e in range(n_epochs):
            if verbose and not self._use_tqdm:
                print('Epoch: {:d}'.format(e))

            epoch_errs = np.zeros((n_batches,))
            epoch_errs_ptr = 0

            if shuffle:
                inds = np.arange(n_data)
                np.random.shuffle(inds)
                data_ve = data_ve[inds]

            r_batches = range(n_batches)

            if verbose and self._use_tqdm:
                r_batches = self._tqdm(r_batches, desc='Epoch: {:d}'.format(e), ascii=True, file=sys.stdout)

            for b in r_batches:
                batch_ve = data_ve[b * batch_size:(b + 1) * batch_size]

                # TODO: avoid splitting this up here
                batch_v = batch_ve[:, :self.n_visible]
                batch_e = batch_ve[:, self.n_visible:]

                self.mutual_info = self.partial_fit([batch_v, batch_e])
                epoch_errs[epoch_errs_ptr] = -self.mutual_info # <=> batch_error
                epoch_errs_ptr += 1

            if verbose:
                err_mean = epoch_errs.mean()
                if self._use_tqdm:
                    self._tqdm.write('Train error: {:.4f}'.format(err_mean))
                    self._tqdm.write('')
                else:
                    print('Train error: {:.4f}'.format(err_mean))
                    print('')
                    sys.stdout.flush()

            errs = np.hstack([errs, epoch_errs])

        return errs

max_stars_repo_name: jqhoogland/rgpy

max_issues_repo_name: jqhoogland/rgpy

max_forks_repo_name: jqhoogland/rgpy


text:
# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*-
# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4
#
# MDAnalysis --- https://www.mdanalysis.org
# Copyright (c) 2006-2017 The MDAnalysis Development Team and contributors
# (see the file AUTHORS for the full list of names)
#
# Released under the GNU Public Licence, v2 or any higher version
#
# Please cite your use of MDAnalysis in published work:
#
# R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler,
# D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein.
# MDAnalysis: A Python package for the rapid analysis of molecular dynamics
# simulations. In S. Benthall and S. Rostrup editors, Proceedings of the 15th
# Python in Science Conference, pages 102-109, Austin, TX, 2016. SciPy.
# doi: 10.25080/majora-629e541a-00e
#
# N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein.
# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations.
# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787
#

"""Correlations utilities --- :mod:`MDAnalysis.lib.correlations`
=================================================================================

:Authors: Paul Smith & Mateusz Bieniek
:Year: 2020
:Copyright: GNU Public License v2

.. versionadded:: 1.0.0

This module is primarily for internal use by other analysis modules. It
provides functionality for calculating the time autocorrelation function
of a binary variable (i.e one that is either true or false at each
frame for a given atom/molecule/set of molecules). This module includes
functions for calculating both the time continuous autocorrelation and
the intermittent autocorrelation. The function :func:`autocorrelation`
calculates the continuous autocorrelation only. The data may be
pre-processed using the function :func:`intermittency` in order to
acount for intermittency before passing the results to
:func:`autocorrelation`.

This module is inspired by seemingly disparate analyses that rely on the same
underlying calculation, including the survival probability of water around
proteins [Araya-Secchi2014]_, hydrogen bond lifetimes [Gowers2015]_, [Araya-Secchi2014]_,
and the rate of cholesterol flip-flop in lipid bilayers [Gu2019]_.

.. seeAlso::

    Analysis tools that make use of modules:

        * :class:`MDAnalysis.analysis.waterdynamics.SurvivalProbability`
            Calculates the continuous or intermittent survival probability
            of an atom group in a region of interest.

        * :class:`MDAnalysis.analysis.hbonds.hbond_analysis`
            Calculates the continuous or intermittent hydrogen bond
            lifetime.

.. rubric:: References

.. [Gu2019] Gu, R.-X.; Baoukina, S.; Tieleman, D. P. (2019)
            Cholesterol Flip-Flop in Heterogeneous Membranes.
            J. Chem. Theory Comput. 15 (3), 2064–2070.
            https://doi.org/10.1021/acs.jctc.8b00933.

"""

from __future__ import absolute_import, division
import numpy as np
from copy import deepcopy

def autocorrelation(list_of_sets, tau_max, window_step=1):
    r"""Implementation of a discrete autocorrelation function.

    The autocorrelation of a property :math:`x` from a time :math:`t=t_0` to :math:`t=t_0 + \tau`
    is given by:

    .. math::
        C(\tau) = \langle \frac{ x(t_0)x(t_0 +\tau) }{ x(t_0)x(t_0) } \rangle

    where :math:`x` may represent any property of a particle, such as velocity or
    potential energy.

    This function is an implementation of a special case of the time
    autocorrelation function in which the property under consideration can
    be encoded with indicator variables, :math:`0` and :math:`1`, to represent the binary
    state of said property. This special case is often referred to as the
    survival probability (:math:`S(\tau)`). As an example, in calculating the survival
    probability of water molecules within 5 Å of a protein, each water
    molecule will either be within this cutoff range (:math:`1`) or not (:math:`0`). The
    total number of water molecules within the cutoff at time :math:`t_0` will be
    given by :math:`N(t_0)`. Other cases include the Hydrogen Bond Lifetime as
    well as the translocation rate of cholesterol across a bilayer.

    The survival probability of a property of a set of particles is
    given by:

    .. math::
        S(\tau) =  \langle \frac{ N(t_0, t_0 + \tau )} { N(t_0) }\rangle

    where :math:`N(t0)` is the number of particles at time :math:`t_0` for which the feature
    is observed, and :math:`N(t0, t_0 + \tau)` is the number of particles for which
    this feature is present at every frame from :math:`t_0` to :math:`N(t0, t_0 + \tau)`.
    The angular brackets represent an average over all time origins, :math:`t_0`.

    See [Araya-Secchi2014]_ for a description survival probability.

    Parameters
    ----------
    list_of_sets : list
      List of sets. Each set corresponds to data from a single frame. Each element in a set
      may be, for example, an atom id or a tuple of atoms ids. In the case of calculating the
      survival probability of water around a protein, these atom ids in a given set will be
      those of the atoms which are within a cutoff distance of the protein at a given frame.
    tau_max : int
      The last tau (lag time, inclusive) for which to calculate the autocorrelation. e.g if tau_max = 20,
      the survival probability will be calculated over 20 frames.
    window_step : int, optional
      The step size for t0 to perform autocorrelation. Ideally, window_step will be larger than
       tau_max to ensure independence of each window for which the calculation is performed.
       Default is 1.

    Returns
    --------
    tau_timeseries : list of int
        the values of tau for which the autocorrelation was calculated
    timeseries : list of int
        the autocorelation values for each of the tau values
    timeseries_data : list of list of int
        the raw data from which the autocorrelation is computed, i.e :math:`S(\tau)` at each window.
        This allows the time dependant evolution of :math:`S(\tau)` to be investigated.

    .. versionadded:: 0.19.2
    """

    # check types
    if (type(list_of_sets) != list and len(list_of_sets) != 0) or type(list_of_sets[0]) != set:
        raise TypeError("list_of_sets must be a one-dimensional list of sets")  # pragma: no cover

    # Check dimensions of parameters
    if len(list_of_sets) < tau_max:
        raise ValueError("tau_max cannot be greater than the length of list_of_sets") # pragma: no cover

    tau_timeseries = list(range(1, tau_max + 1))
    timeseries_data = [[] for _ in range(tau_max)]

    # calculate autocorrelation
    for t in range(0, len(list_of_sets), window_step):
        Nt = len(list_of_sets[t])

        if Nt == 0:
            continue

        for tau in tau_timeseries:
            if t + tau >= len(list_of_sets):
                break

            # continuous: IDs that survive from t to t + tau and at every frame in between
            Ntau = len(set.intersection(*list_of_sets[t:t + tau + 1]))
            timeseries_data[tau - 1].append(Ntau / float(Nt))

    timeseries = [np.mean(x) for x in timeseries_data]

    # at time 0 the value has to be one
    tau_timeseries.insert(0, 0)
    timeseries.insert(0, 1)

    return tau_timeseries, timeseries, timeseries_data

def correct_intermittency(list_of_sets, intermittency):
    r"""Preprocess data to allow intermittent behaviour prior to calling :func:`autocorrelation`.

    Survival probabilty may be calculated either with a strict continuous requirement or
    a less strict intermittency. If calculating the survival probability water around a
    protein for example, in the former case the water must be within a cutoff distance
    of the protein at every frame from :math:`t_0` to :math:`t_0 + \tau` in order for it to be considered
    present at :math:`t_0 + \tau`. In the intermittent case, the water molecule is allowed to
    leave the region of interest for up to a specified consecutive number of frames whilst still
    being considered present at :math:`t_0 + \tau`.

    This function pre-processes data, such as the atom ids of water molecules within a cutoff
    distance of a protein at each frame, in order to allow for intermittent behaviour, with a
    single pass over the data.

    For example, if an atom is absent for a number of frames equal or smaller than the parameter
    :attr:`intermittency`, then this absence will be removed and thus the atom is considered to have
    not left.
    e.g 7,A,A,7 with `intermittency=2` will be replaced by 7,7,7,7, where A=absence.

    The returned data can be used as input to the function :func:`autocorrelation` in order
    to calculate the survival probability with a given intermittency.

    See [Gowers2015]_ for a description of
    intermittency in the calculation of hydrogen bond lifetimes.

    # TODO - is intermittency consitent with list of sets of sets? (hydrogen bonds)

    Parameters
    ----------
    list_of_sets: list
        In the simple case of e.g survival probability, a list of sets of atom ids present at each frame, where a
        single set contains atom ids at a given frame, e.g [{0, 1}, {0}, {0}, {0, 1}]
    intermittency : int
        The maximum gap allowed. The default `intermittency=0` means that if the datapoint is missing at any frame, no
        changes are made to the data. With the value of `intermittency=2`, all datapoints missing for up to two
        consecutive frames will be instead be considered present.

    Returns
    -------
    list_of_sets: list
        returns a new list with the IDs with added IDs which disappeared for <= :attr:`intermittency`.
        e.g If [{0, 1}, {0}, {0}, {0, 1}] is a list of sets of atom ids present at each frame and `intermittency=2`,
        both atoms will be considered present throughout and thus the returned list of sets will be
        [{0, 1}, {0, 1}, {0, 1}, {0, 1}].

    """
    if intermittency == 0:
        return list_of_sets

    list_of_sets = deepcopy(list_of_sets)

    for i, ids in enumerate(list_of_sets):
        # initially update each frame as seen 0 ago (now)
        seen_frames_ago = {i: 0 for i in ids}
        for j in range(1, intermittency + 2):
            for atomid in seen_frames_ago.keys():
                # no more frames
                if i + j >= len(list_of_sets):
                    continue

                # if the atom is absent now
                if not atomid in list_of_sets[i + j]:
                    # increase its absence counter
                    seen_frames_ago[atomid] += 1
                    continue

                # the atom is found
                if seen_frames_ago[atomid] == 0:
                    # the atom was present in the last frame
                    continue

                # it was absent more times than allowed
                if seen_frames_ago[atomid] > intermittency:
                    continue

                # the atom was absent but returned (within <= intermittency_value)
                # add it to the frames where it was absent.
                # Introduce the corrections.
                for k in range(seen_frames_ago[atomid], 0, -1):
                    list_of_sets[i + j - k].add(atomid)

                seen_frames_ago[atomid] = 0
    return list_of_sets

ext: py

max_stars_repo_name: dtklinh/GBRDE

max_issues_repo_name: dtklinh/GBRDE

max_forks_repo_name: dtklinh/GBRDE

reason: import numpy


text:
# -*- coding: utf-8 -*-
#
# Dappled tiling
# by S. Kaji
# 17/Jun/2016
#
# A python implementation of the algorithm discussed in the paper
# "Dappled tiling"
# By S. Kaji and .....
#

import numpy as np
import random

# size of the grid
M = 12
N = 10

# conditions: H[i] is the limit of the length of the horizontal strip with tile i
H = [2,5]
V = [5,2]

# (quick dirty) global variables
# for book-keeping the number of consecutive tiles
VDanger = np.ones((M+1), dtype=np.int)
HDanger = np.ones((M+1), dtype=np.int)

# choose a tile different from t
def chooseTile(t):
    return 1-t

# set barrier
def setbarrier(f,n,m):
    for i in range(n+1):
        f[i,0] = -1
    for j in range(m+1):
        f[0,j] = -1
    return f

# give a randomised tiling
def randtiling(n,m):
    f = np.zeros((n+1,m+1), dtype=np.int)
    for i in range(1,n+1):
        for j in range(1,m+1):
            f[i,j] = random.randint(0,1)
    return setbarrier(f,n,m)

# give an all-one tiling
def onetiling(n,m):
    f = np.ones((n+1,m+1), dtype=np.int)
    return setbarrier(f,n,m)

# give a diagonal tiling
def diagtiling(n,m):
    f = np.ones((n+1,m+1), dtype=np.int)
    for i in range(1,n+1):
        for j in range(i,min(i+int(m/2),m+1)):
            f[i,j] = 0
    return setbarrier(f,n,m)

# compute the number of consecutive tiles of lower weights
def compute_danger(f,i,j):
    if f[i,j]==f[i-1,j]:
        VDanger[j] = VDanger[j]+1
    else:
        VDanger[j] = 1
    if f[i,j]==f[i,j-1]:
        HDanger[j] = HDanger[j-1]+1
    else:
        HDanger[j] = 1

# print tiling patterns
def printtiling(f):
    print(f[1:N+1,1:M+1])

# convert any tiling to a dappled one
def dappled(f):
    for w in range(1,N+M+1):
        for i in range(max(1,w-M),min(w,N+1)):
            j = w-i
            prev_vd=VDanger[j]
            compute_danger(f,i,j)
            if V[f[i,j]] < VDanger[j] or H[f[i,j]] < HDanger[j]:
                f[i,j]=chooseTile(f[i,j])
                VDanger[j]=prev_vd
                compute_danger(f,i,j)
                if V[f[i,j]] < VDanger[j]:
                    f[i,j]=f[i-1,j-1]
                    f[i-1,j]=f[i,j-1]=chooseTile(f[i,j])
                    HDanger[j]=VDanger[j]=VDanger[j-1]=1
                    if j < M:
                        if f[i-1,j]==f[i-1,j+1]:
                            HDanger[j+1]=2
                elif H[f[i,j]] < HDanger[j]:
                    f[i,j]=f[i-1,j-1]
                    f[i-1,j]=f[i,j-1]=chooseTile(f[i,j])
                    VDanger[j]=HDanger[j]=VDanger[j-1]=1
                    if j < M:
                        if f[i-1,j]==f[i-1,j+1]:
                            HDanger[j+1]=2

## main
if __name__ == "__main__":
    print((M,N))
    f = randtiling(N,M)
#    f = diagtiling(N,M)
#    f = onetiling(N,M)
    print("Initial tiling")
    printtiling(f)

    dappled(f)

    print("Dappled tiling")
    printtiling(f)
    exit

max_stars_repo_name: shizuo-kaji/DappledTiling

max_issues_repo_name: shizuo-kaji/DappledTiling

max_forks_repo_name: shizuo-kaji/DappledTiling


text:
"""Multi-object Tracking evaluation code."""
import argparse
import json
import time
from functools import partial
from multiprocessing import Pool
from typing import Callable, Dict, List, Tuple, Union

import motmetrics as mm
import numpy as np

from ..common.io import open_write_text
from ..common.logger import logger
from ..common.parallel import NPROC
from ..common.typing import NDArrayF64, NDArrayI32, NDArrayU8
from ..label.io import group_and_sort, load, load_label_config
from ..label.transforms import box2d_to_bbox
from ..label.typing import Category, Config, Frame, Label
from ..label.utils import (
    check_crowd,
    check_ignored,
    get_leaf_categories,
    get_parent_categories,
)
from .result import AVERAGE, OVERALL, Result, Scores, ScoresList
from .utils import handle_inconsistent_length, label_ids_to_int

# Video = TypeVar("Video", List[Frame])
VidFunc = Callable[
    [List[Frame], List[Frame], List[str], float, float, bool],
    List[mm.MOTAccumulator],
]

METRIC_MAPS = {
    "idf1": "IDF1",
    "mota": "MOTA",
    "motp": "MOTP",
    "num_false_positives": "FP",
    "num_misses": "FN",
    "num_switches": "IDSw",
    "mostly_tracked": "MT",
    "partially_tracked": "PT",
    "mostly_lost": "ML",
    "num_fragmentations": "FM",
}
METRIC_TO_AVERAGE = set(["IDF1", "MOTA", "MOTP"])

class TrackResult(Result):
    """The class for bounding box tracking evaluation results."""

    mMOTA: float
    mMOTP: float
    mIDF1: float
    MOTA: List[Dict[str, float]]
    MOTP: List[Dict[str, float]]
    IDF1: List[Dict[str, float]]
    FP: List[Dict[str, int]]
    FN: List[Dict[str, int]]
    IDSw: List[Dict[str, int]]
    MT: List[Dict[str, int]]
    PT: List[Dict[str, int]]
    ML: List[Dict[str, int]]
    FM: List[Dict[str, int]]

    # pylint: disable=useless-super-delegation
    def __eq__(self, other: "Result") -> bool:  # type: ignore
        """Check whether two instances are equal."""
        return super().__eq__(other)

def parse_objects(
    objects: List[Label], classes: List[str], ignore_unknown_cats: bool = False
) -> Tuple[NDArrayF64, NDArrayI32, NDArrayI32, NDArrayF64]:
    """Parse objects under Scalabel formats."""
    bboxes, labels, ids, ignore_bboxes = [], [], [], []
    for obj in objects:
        box_2d = obj.box2d
        if box_2d is None:
            continue
        bbox = box2d_to_bbox(box_2d)
        category = obj.category
        if category in classes:
            if check_crowd(obj) or check_ignored(obj):
                ignore_bboxes.append(bbox)
            else:
                bboxes.append(bbox)
                labels.append(classes.index(category))
                ids.append(obj.id)
        else:
            if not ignore_unknown_cats:
                raise KeyError(f"Unknown category: {category}")
    bboxes_arr: NDArrayF64 = np.array(bboxes, dtype=np.float64)
    labels_arr: NDArrayI32 = np.array(labels, dtype=np.int32)
    ids_arr: NDArrayI32 = np.array(ids, dtype=np.int32)
    ignore_bboxes_arr: NDArrayF64 = np.array(ignore_bboxes, dtype=np.float64)
    return (bboxes_arr, labels_arr, ids_arr, ignore_bboxes_arr)

def intersection_over_area(preds: NDArrayF64, gts: NDArrayF64) -> NDArrayF64:
    """Returns the intersection over the area of the predicted box."""
    out: NDArrayF64 = np.zeros((len(preds), len(gts)), dtype=np.float64)
    for i, p in enumerate(preds):
        for j, g in enumerate(gts):
            w = min(p[0] + p[2], g[0] + g[2]) - max(p[0], g[0])
            h = min(p[1] + p[3], g[1] + g[3]) - max(p[1], g[1])
            out[i][j] = max(w, 0) * max(h, 0) / float(p[2] * p[3])
    return out

def acc_single_video_mot(
    gts: List[Frame],
    results: List[Frame],
    classes: List[str],
    iou_thr: float = 0.5,
    ignore_iof_thr: float = 0.5,
    ignore_unknown_cats: bool = False,
) -> List[mm.MOTAccumulator]:
    """Accumulate results for one video."""
    assert len(gts) == len(results)

    def get_frame_index(frame: Frame) -> int:
        return frame.frameIndex if frame.frameIndex is not None else 0

    num_classes = len(classes)
    gts = sorted(gts, key=get_frame_index)
    results = sorted(results, key=get_frame_index)
    accs = [mm.MOTAccumulator(auto_id=True) for _ in range(num_classes)]

    label_ids_to_int(gts)

    for gt, result in zip(gts, results):
        assert gt.frameIndex == result.frameIndex
        gt_bboxes, gt_labels, gt_ids, gt_ignores = parse_objects(
            gt.labels if gt.labels is not None else [],
            classes,
            ignore_unknown_cats,
        )
        pred_bboxes, pred_labels, pred_ids, _ = parse_objects(
            result.labels if result.labels is not None else [],
            classes,
            ignore_unknown_cats,
        )
        for i in range(num_classes):
            gt_inds, pred_inds = gt_labels == i, pred_labels == i
            gt_bboxes_c, gt_ids_c = gt_bboxes[gt_inds], gt_ids[gt_inds]
            pred_bboxes_c, pred_ids_c = (
                pred_bboxes[pred_inds],
                pred_ids[pred_inds],
            )
            if gt_bboxes_c.shape[0] == 0 and pred_bboxes_c.shape[0] != 0:
                distances = np.full((0, pred_bboxes_c.shape[0]), np.nan)
            elif gt_bboxes_c.shape[0] != 0 and pred_bboxes_c.shape[0] == 0:
                distances = np.full((gt_bboxes_c.shape[0], 0), np.nan)
            else:
                distances = mm.distances.iou_matrix(
                    gt_bboxes_c, pred_bboxes_c, max_iou=1 - iou_thr
                )
            if gt_ignores.shape[0] > 0:
                # 1. assign gt and preds
                fps: NDArrayU8 = np.ones(pred_bboxes_c.shape[0], dtype=bool)
                le, ri = mm.lap.linear_sum_assignment(distances)
                for m, n in zip(le, ri):
                    if np.isfinite(distances[m, n]):
                        fps[n] = False
                # 2. ignore by iof
                iofs = intersection_over_area(pred_bboxes_c, gt_ignores)
                ignores: bool = np.greater(iofs, ignore_iof_thr).any(axis=1)
                # 3. filter preds
                valid_inds = np.logical_not(np.logical_and(fps, ignores))
                pred_ids_c = pred_ids_c[valid_inds]
                distances = distances[:, valid_inds]
            if distances.shape != (0, 0):
                accs[i].update(gt_ids_c, pred_ids_c, distances)
    return accs

def aggregate_accs(
    video_accs: List[List[mm.MOTAccumulator]],
    classes: List[Category],
    super_classes: Dict[str, List[Category]],
) -> Tuple[List[str], List[List[str]], List[List[mm.MOTAccumulator]]]:
    """Aggregate the results of the entire dataset."""
    # accs for each class
    class_names: List[str] = [c.name for c in classes]
    metric_names: List[List[str]] = [[] for _ in classes]
    class_accs: List[List[mm.MOTAccumulator]] = [[] for _ in classes]
    for video_ind, _accs in enumerate(video_accs):
        for cls_ind, acc in enumerate(_accs):
            if (
                len(acc._events["Type"])  # pylint: disable=protected-access
                == 0
            ):
                continue
            name = f"{classes[cls_ind].name}_{video_ind}"
            class_accs[cls_ind].append(acc)
            metric_names[cls_ind].append(name)

    # super categories (if any)
    for super_cls, cls in super_classes.items():
        class_names.append(super_cls)
        metric_names.append(
            [n for c in cls for n in metric_names[classes.index(c)]]
        )
        class_accs.append(
            [a for c in cls for a in class_accs[classes.index(c)]]
        )

    # overall
    class_names.append(OVERALL)
    metric_names.append(
        [n for name in metric_names[: len(classes)] for n in name]
    )
    class_accs.append([a for acc in class_accs[: len(classes)] for a in acc])

    return class_names, metric_names, class_accs

def evaluate_single_class(
    names: List[str], accs: List[mm.MOTAccumulator]
) -> Dict[str, Union[int, float]]:
    """Evaluate results for one class.

    Args:
        names (list[str]): list of metric names
        accs (list[pymotmetrics.MOTAAccumulator]): list of accumulators
    Return:
        flat_dict (dict[str, int | float]):
            the dict that maps metric to score for the given class
    """
    mh = mm.metrics.create()
    summary = mh.compute_many(
        accs, names=names, metrics=METRIC_MAPS.keys(), generate_overall=True
    )
    flat_dict: Dict[str, Union[int, float]] = {
        METRIC_MAPS[k]: v["OVERALL"] for k, v in summary.to_dict().items()
    }
    if np.isnan(flat_dict["MOTP"]):
        num_dets = mh.compute_many(
            accs,
            names=names,
            metrics=["num_detections"],
            generate_overall=True,
        )
        sum_motp = (summary["motp"] * num_dets["num_detections"]).sum()
        motp = mm.math_util.quiet_divide(
            sum_motp, num_dets["num_detections"]["OVERALL"]
        )
        flat_dict["MOTP"] = float(1 - motp)
    else:
        flat_dict["MOTP"] = 1 - flat_dict["MOTP"]

    for metric, score in flat_dict.items():
        if isinstance(score, float):
            flat_dict[metric] = 100 * score
    return flat_dict

def compute_average(
    flat_dicts: List[Dict[str, Union[int, float]]],
    metrics: List[str],
    classes: List[Category],
) -> Dict[str, Union[int, float]]:
    """Calculate the AVERAGE scores."""
    ave_dict: Dict[str, Union[int, float]] = {}
    for metric in metrics:
        dtype = type(flat_dicts[-1][metric])
        v: NDArrayF64 = np.array(
            [flat_dicts[i][metric] for i in range(len(classes))],
            dtype=np.float64,
        )
        v = np.nan_to_num(v, nan=0, posinf=0, neginf=0)
        if dtype == int:
            value = int(v.sum())  # type: Union[int, float]
        elif dtype == float:
            value = float(v.mean())
        else:
            raise TypeError()
        ave_dict[metric] = value
    return ave_dict

def generate_results(
    flat_dicts: List[Dict[str, Union[int, float]]],
    class_names: List[str],
    metrics: List[str],
    classes: List[Category],
    super_classes: Dict[str, List[Category]],
) -> TrackResult:
    """Compute summary metrics for evaluation results."""
    ave_dict = compute_average(flat_dicts, metrics, classes)
    class_names.insert(len(flat_dicts) - 1, AVERAGE)
    flat_dicts.insert(len(flat_dicts) - 1, ave_dict)

    nested_dict: Dict[str, Scores] = {
        metric: {
            class_name: res_dict_[metric]
            for class_name, res_dict_ in zip(class_names, flat_dicts)
        }
        for metric in metrics
    }

    basic_set = [c.name for c in classes]
    super_set = list(super_classes.keys())
    hyper_set = [AVERAGE, OVERALL]
    class_name_sets = [basic_set, hyper_set]
    if [name for name in class_names if name in super_classes]:
        class_name_sets.insert(1, super_set)

    res_dict: Dict[str, Union[int, float, ScoresList]] = {
        metric: [
            {
                class_name: score
                for class_name, score in scores.items()
                if class_name in class_name_set
            }
            for class_name_set in class_name_sets
        ]
        for metric, scores in nested_dict.items()
    }
    res_dict.update(
        {"m" + metric: ave_dict[metric] for metric in METRIC_TO_AVERAGE}
    )
    return TrackResult(**res_dict)

def evaluate_track(
    acc_single_video: VidFunc,
    gts: List[List[Frame]],
    results: List[List[Frame]],
    config: Config,
    iou_thr: float = 0.5,
    ignore_iof_thr: float = 0.5,
    ignore_unknown_cats: bool = False,
    nproc: int = NPROC,
) -> TrackResult:
    """Evaluate CLEAR MOT metrics for a Scalabel format dataset.

    Args:
        acc_single_video: Function for calculating metrics over a single video.
        gts: the ground truth annotations in Scalabel format
        results: the prediction results in Scalabel format.
        config: Config object
        iou_thr: Minimum IoU for a bounding box to be considered a positive.
        ignore_iof_thr: Min. Intersection over foreground with ignore regions.
        ignore_unknown_cats: if False, raise KeyError when trying to evaluate
            unknown categories.
        nproc: processes number for loading files

    Returns:
        TrackResult: evaluation results.
    """
    logger.info("Tracking evaluation with CLEAR MOT metrics.")
    t = time.time()
    results = handle_inconsistent_length(gts, results)
    assert len(gts) == len(results)

    classes = get_leaf_categories(config.categories)
    super_classes = get_parent_categories(config.categories)

    logger.info("evaluating...")
    class_names = [c.name for c in classes]
    if nproc > 1:
        with Pool(nproc) as pool:
            video_accs = pool.starmap(
                partial(
                    acc_single_video,
                    classes=class_names,
                    ignore_iof_thr=ignore_iof_thr,
                    ignore_unknown_cats=ignore_unknown_cats,
                ),
                zip(gts, results),
            )
    else:
        video_accs = [
            acc_single_video(
                gt,
                result,
                class_names,
                iou_thr,
                ignore_iof_thr,
                ignore_unknown_cats,
            )
            for gt, result in zip(gts, results)
        ]

    class_names, metric_names, class_accs = aggregate_accs(
        video_accs, classes, super_classes
    )

    logger.info("accumulating...")
    if nproc > 1:
        with Pool(nproc) as pool:
            flat_dicts = pool.starmap(
                evaluate_single_class, zip(metric_names, class_accs)
            )
    else:
        flat_dicts = [
            evaluate_single_class(names, accs)
            for names, accs in zip(metric_names, class_accs)
        ]

    metrics = list(METRIC_MAPS.values())
    result = generate_results(
        flat_dicts, class_names, metrics, classes, super_classes
    )
    t = time.time() - t
    logger.info("evaluation finishes with %.1f s.", t)
    return result

def parse_arguments() -> argparse.Namespace:
    """Parse the arguments."""
    parser = argparse.ArgumentParser(description="MOT evaluation.")
    parser.add_argument(
        "--gt", "-g", required=True, help="path to mot ground truth"
    )
    parser.add_argument(
        "--result", "-r", required=True, help="path to mot results"
    )
    parser.add_argument(
        "--config",
        "-c",
        default=None,
        help="Path to config toml file. Contains definition of categories, "
        "and optionally attributes as well as resolution. For an example "
        "see scalabel/label/testcases/configs.toml",
    )
    parser.add_argument(
        "--out-file",
        default="",
        help="Output path for mot evaluation results.",
    )
    parser.add_argument(
        "--iou-thr",
        type=float,
        default=0.5,
        help="iou threshold for mot evaluation",
    )
    parser.add_argument(
        "--ignore-iof-thr",
        type=float,
        default=0.5,
        help="ignore iof threshold for mot evaluation",
    )
    parser.add_argument(
        "--ignore-unknown-cats",
        type=bool,
        default=False,
        help="ignore unknown categories for mot evaluation",
    )
    parser.add_argument(
        "--nproc",
        "-p",
        type=int,
        default=NPROC,
        help="number of processes for mot evaluation",
    )
    return parser.parse_args()

if __name__ == "__main__":
    args = parse_arguments()
    dataset = load(args.gt)
    gt_frames, cfg = dataset.frames, dataset.config
    if args.config is not None:
        cfg = load_label_config(args.config)
    if cfg is None:
        raise ValueError(
            "Dataset config is not specified. Please use --config"
            " to specify a config for this dataset."
        )
    eval_result = evaluate_track(
        acc_single_video_mot,
        group_and_sort(gt_frames),
        group_and_sort(load(args.result).frames),
        cfg,
        args.iou_thr,
        args.ignore_iof_thr,
        args.ignore_unknown_cats,
        args.nproc,
    )
    logger.info(eval_result)
    logger.info(eval_result.summary())
    if args.out_file:
        with open_write_text(args.out_file) as fp:
            json.dump(eval_result.dict(), fp, indent=2)

max_stars_repo_name: Anthuang/scalabel

max_issues_repo_name: Anthuang/scalabel

max_forks_repo_name: Anthuang/scalabel


text:
import unittest
import warnings
import os
import string
import random
import getpass
import sys

import atpy
import numpy as np
np.seterr(all='ignore')

# Size of the test table
shape = (100, )
shape_vector = (100, 10)

# SQL connection parameters
username = getpass.getuser()
password = "C2#uwQsk"

def random_int_array(dtype, shape):
    n = np.product(shape)
    n = n * np.iinfo(dtype).bits / 8
    s = "".join(chr(random.randrange(0, 256)) for i in xrange(n))
    return np.fromstring(s, dtype=dtype).reshape(shape)

def random_float_array(dtype, shape):
    n = np.product(shape)
    if dtype == np.float32:
        n = n * 4
    else:
        n = n * 8
    s = "".join(chr(random.randrange(0, 256)) for i in xrange(n))
    array = np.fromstring(s, dtype=dtype)
    if np.sum(np.isnan(array)):
        array[np.isnan(array)] = random_float_array(dtype, array[np.isnan(array)].shape)
    return array.reshape(shape)

numpy_types = [np.bool, np.int8, np.int16, np.int32, np.int64, np.uint8, \
               np.uint16, np.uint32, np.uint64, np.float32, np.float64, \
               np.dtype('|S100')]

def random_generic(dtype, name, shape):

    if 'int' in name:
        values = random_int_array(dtype, shape)
    elif 'float' in name:
        values = random_float_array(dtype, shape)
    elif 'bool' in name:
        values = np.random.random(shape) > 0.5
    else:
        values = np.zeros(shape, dtype=dtype)
        for i in range(shape[0]):
            s = ""
            for k in range(dtype.itemsize):
                s += random.choice(string.letters + string.digits)
            values[i] = s

    return values

def generate_simple_table(dtype, shape):

    table = atpy.Table(name='atpy_test')

    try:
        name = dtype.__name__
    except:
        name = dtype.type.__name__

    values = random_generic(dtype, name, shape)

    table.add_column('col_'+name, values, dtype=dtype)

    return table

class ColumnsDefaultTestCase():

    def test_uint8(self):
        self.generic_test(np.uint8)

    def test_uint16(self):
        self.generic_test(np.uint16)

    def test_uint32(self):
        self.generic_test(np.uint32)

    def test_uint64(self):
        self.generic_test(np.uint64)

    def test_int8(self):
        self.generic_test(np.int8)

    def test_int16(self):
        self.generic_test(np.int16)

    def test_int32(self):
        self.generic_test(np.int32)

    def test_int64(self):
        self.generic_test(np.int64)

    def test_float32(self):
        self.generic_test(np.float32)

    def test_float64(self):
        self.generic_test(np.float64)

    def test_string(self):
        self.generic_test(np.dtype('|S100'))

class EmptyColumnsTestCase(unittest.TestCase, ColumnsDefaultTestCase):

    def generic_test(self, dtype):
        try:
            t = atpy.Table()
            t.add_empty_column('a', dtype, shape=shape)
            t.add_empty_column('b', dtype)
            t.add_empty_column('c', dtype)
        except:
            self.fail(sys.exc_info()[1])

class EmptyVectorColumnsTestCase(unittest.TestCase, ColumnsDefaultTestCase):

    def generic_test(self, dtype):
        try:
            t = atpy.Table()
            t.add_empty_column('a', dtype, shape=shape_vector)
            t.add_empty_column('b', dtype)
            t.add_empty_column('c', dtype, shape=shape_vector)
        except:
            self.fail(sys.exc_info()[1])

class DefaultTestCase():

    def assertAlmostEqualSig(test, first, second, significant=7, msg=None):
        ratio = first / second
        if np.abs(ratio - 1.) > 10.**(-significant+1):
            raise unittest.TestCase.failureException(msg or '%r != %r within %r significant digits' % (first, second, significant))

    def integer_test(self, dtype):
        colname = 'col_%s' % dtype.__name__
        self.writeread(dtype)
        before, after = self.table_orig.data[colname], self.table_new.data[colname]
        self.assertEqual(before.shape, after.shape)
        if before.ndim == 1:
            for i in range(before.shape[0]):
                self.assertEqual(before[i], after[i])
        else:
            for i in range(before.shape[0]):
                for j in range(before.shape[1]):
                    self.assertEqual(before[i, j], after[i, j])

    def test_uint8(self):
        self.integer_test(np.uint8)

    def test_uint16(self):
        self.integer_test(np.uint16)

    def test_uint32(self):
        self.integer_test(np.uint32)

    def test_uint64(self):
        self.integer_test(np.uint64)

    def test_int8(self):
        self.integer_test(np.int8)

    def test_int16(self):
        self.integer_test(np.int16)

    def test_int32(self):
        self.integer_test(np.int32)

    def test_int64(self):
        self.integer_test(np.int64)

    def float_test(self, dtype, significant=None):
        colname = 'col_%s' % dtype.__name__
        self.writeread(dtype)
        before, after = self.table_orig.data[colname], self.table_new.data[colname]
        self.assertEqual(before.shape, after.shape)
        self.assertEqual(before.dtype.type, after.dtype.type)
        if before.ndim == 1:
            for i in range(before.shape[0]):
                if(np.isnan(before[i])):
                    self.failUnless(np.isnan(after[i]))
                elif(np.isinf(before[i])):
                    self.failUnless(np.isinf(after[i]))
                else:
                    if significant:
                        self.assertAlmostEqualSig(before[i], after[i], significant=significant)
                    else:
                        self.assertEqual(before[i], after[i])
        else:
            for i in range(before.shape[0]):
                for j in range(before.shape[1]):
                    if(np.isnan(before[i, j])):
                        self.failUnless(np.isnan(after[i, j]))
                    elif(np.isinf(before[i, j])):
                        self.failUnless(np.isinf(after[i, j]))
                    else:
                        if significant:
                            self.assertAlmostEqualSig(before[i, j], after[i, j], significant=significant)
                        else:
                            self.assertEqual(before[i, j], after[i, j])

    def test_float32(self):
        if self.format == 'mysql':
            self.float_test(np.float32, significant=6)
        elif self.format == 'postgres':
            self.float_test(np.float32, significant=6)
        else:
            self.float_test(np.float32)

    def test_float64(self):
        if self.format == 'mysql':
            self.float_test(np.float64, significant=15)
        elif self.format == 'postgres':
            self.float_test(np.float64, significant=12)
        else:
            self.float_test(np.float64)

    def test_string(self):
        self.writeread(np.dtype('|S100'))
        before, after = self.table_orig.data['col_string_'], self.table_new.data['col_string_']
        self.assertEqual(before.shape, after.shape)
        self.assertEqual(before.dtype.type, after.dtype.type)
        for i in range(len(self.table_orig)):
            self.assertEqual(before[i], after[i])

try:
    import pyfits

    class FITSTestCase(unittest.TestCase, DefaultTestCase):

        format = 'fits'

        test_uint64 = None # unsupported

        def writeread(self, dtype):

            self.table_orig = generate_simple_table(dtype, shape)
            self.table_orig.write('test_atpy.fits', verbose=False, overwrite=True)
            self.table_new = atpy.Table('test_atpy.fits', verbose=False)
            os.remove('test_atpy.fits')

    class FITSTestCaseVector(unittest.TestCase, DefaultTestCase):

        format = 'fits'

        test_string = None # unsupported
        test_uint64 = None # unsupported

        def writeread(self, dtype):

            self.table_orig = generate_simple_table(dtype, shape_vector)
            self.table_orig.write('test_atpy.fits', verbose=False, overwrite=True)
            self.table_new = atpy.Table('test_atpy.fits', verbose=False)
            os.remove('test_atpy.fits')

except:
    print "pyfits is not installed - skipping FITS tests"
    pass

try:
    import h5py

    class HDF5TestCase(unittest.TestCase, DefaultTestCase):

        format = 'hdf5'

        def writeread(self, dtype):

            self.table_orig = generate_simple_table(dtype, shape)
            self.table_orig.write('test_atpy.hdf5', verbose=False, overwrite=True)
            self.table_new = atpy.Table('test_atpy.hdf5', verbose=False)
            os.remove('test_atpy.hdf5')

    class HDF5TestCaseVector(unittest.TestCase, DefaultTestCase):

        format = 'hdf5'

        test_string = None # unsupported

        def writeread(self, dtype):

            self.table_orig = generate_simple_table(dtype, shape_vector)
            self.table_orig.write('test_atpy.hdf5', verbose=False, overwrite=True)
            self.table_new = atpy.Table('test_atpy.hdf5', verbose=False)
            os.remove('test_atpy.hdf5')

except:
    print "h5py is not installed - skipping HDF5 tests"
    pass

try:
    import vo

    class VOTestCase(unittest.TestCase, DefaultTestCase):

        format = 'vo'

        test_uint64 = None # unsupported

        def writeread(self, dtype):

            self.table_orig = generate_simple_table(dtype, shape)
            self.table_orig.write('test_atpy.xml', verbose=False, overwrite=True)
            self.table_new = atpy.Table('test_atpy.xml', verbose=False)
            os.remove('test_atpy.xml')

    class VOTestCaseVector(unittest.TestCase, DefaultTestCase):

        format = 'vo'

        test_string = None # unsupported
        test_uint64 = None # unsupported

        def writeread(self, dtype):

            self.table_orig = generate_simple_table(dtype, shape_vector)
            self.table_orig.write('test_atpy.xml', verbose=False, overwrite=True)
            self.table_new = atpy.Table('test_atpy.xml', verbose=False)
            os.remove('test_atpy.xml')

except:
    print "vo is not installed - skipping MySQL tests"
    pass

class IPACTestCase(unittest.TestCase, DefaultTestCase):

    format = 'ipac'

    def writeread(self, dtype):

        self.table_orig = generate_simple_table(dtype, shape)
        self.table_orig.write('test_atpy.tbl', verbose=False, overwrite=True)
        self.table_new = atpy.Table('test_atpy.tbl', verbose=False)
        os.remove('test_atpy.tbl')

try:
    import sqlite3

    class SQLiteTestCase(unittest.TestCase, DefaultTestCase):

        format = 'sqlite'

        test_uint64 = None # unsupported

        def writeread(self, dtype):

            self.table_orig = generate_simple_table(dtype, shape)
            self.table_orig.write('sqlite', 'test_atpy.db', verbose=False, overwrite=True)
            self.table_new = atpy.Table('sqlite', 'test_atpy.db', verbose=False)
            os.remove('test_atpy.db')

    class SQLiteTestCaseQuery(unittest.TestCase, DefaultTestCase):

        format = 'sqlite'

        test_uint64 = None # unsupported

        def writeread(self, dtype):

            self.table_orig = generate_simple_table(dtype, shape)
            self.table_orig.write('sqlite', 'test_atpy.db', verbose=False, overwrite=True)
            self.table_new = atpy.Table('sqlite', 'test_atpy.db', verbose=False, query='select * from atpy_test')
            os.remove('test_atpy.db')

except:
    print "sqlite3 is not installed - skipping SQLite tests"
    pass

try:
    import MySQLdb

    class MySQLTestCase(unittest.TestCase, DefaultTestCase):

        format = 'mysql'

        def writeread(self, dtype):

            self.table_orig = generate_simple_table(dtype, shape)
            self.table_orig.write('mysql', db='python', overwrite=True, verbose=False, user=username, passwd=password)
            self.table_new = atpy.Table('mysql', db='python', verbose=False, user=username, passwd=password, table='atpy_test')

    class MySQLTestCaseQuery(unittest.TestCase, DefaultTestCase):

        format = 'mysql'

        def writeread(self, dtype):

            self.table_orig = generate_simple_table(dtype, shape)
            self.table_orig.write('mysql', db='python', overwrite=True, verbose=False, user=username, passwd=password)
            self.table_new = atpy.Table('mysql', db='python', verbose=False, user=username, passwd=password, query='select * from atpy_test')

except:
    print "MySQLdb is not installed - skipping MySQL tests"
    pass

try:
    import pgdb

    class PostGreSQLTestCase(unittest.TestCase, DefaultTestCase):

        format = 'postgres'

        test_uint64 = None # unsupported

        def writeread(self, dtype):

            self.table_orig = generate_simple_table(dtype, shape)
            self.table_orig.write('postgres', database='python', overwrite=True, verbose=False, user=username, password=password)
            self.table_new = atpy.Table('postgres', database='python', verbose=False, user=username, password=password, table='atpy_test')

    class PostGreSQLTestCaseQuery(unittest.TestCase, DefaultTestCase):

        format = 'postgres'

        test_uint64 = None # unsupported

        def writeread(self, dtype):

            self.table_orig = generate_simple_table(dtype, shape)
            self.table_orig.write('postgres', database='python', overwrite=True, verbose=False, user=username, password=password)
            self.table_new = atpy.Table('postgres', database='python', verbose=False, user=username, password=password, query='select * from atpy_test')

except:
    print "pgdb is not installed - skipping PostGreSQL tests"
    pass

if __name__ == '__main__':
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        unittest.main()

max_stars_repo_name: keflavich/atpy

max_issues_repo_name: keflavich/atpy

max_forks_repo_name: keflavich/atpy


text:
#import matplotlib
#matplotlib.use('Agg')

import numpy as np
from model_selection import cross_validate
from visualization import graphics
from service_package import utils
import argparse as ap

import warnings
warnings.filterwarnings("ignore")

def main(test_size_ratio,num_shuffles,stab_gap,
         X_train_pathway,y_train_pathway,
         feat_names_pathway,example_names_pathway,
         cv_folds,alpha_start,alpha_end,num_alphas,
         smoothness_power_start,smoothness_power_end,
         smoothness_base,num_feat_filtered,data_contamination,
         max_backfit_iter=20,num_jobs_factor=2,
         X_test_pathway=None,y_test_pathway=None):

    print("Data loading")
    [X_train,y_train],[feat_names_full,example_names],[X_test,y_test]=\
    utils.load_data(X_train_pathway,y_train_pathway,feat_names_pathway,
                    example_names_pathway,data_contamination,
                    X_test_pathway,y_test_pathway,num_feat_filtered)

    print("Parameter grid seeding")
    smooth_factors_list=utils.get_smoothness_grid(smoothness_power_start,
                                                  smoothness_power_end,
                                                  smooth_base=smoothness_base)
    alphas_list=utils.get_alphas_grid(alpha_start,alpha_end,num_alphas)

    print("Grid search of parameters via cross-validation and stability selection")
    gridsearch=cross_validate.ModelSelector(smooth_factor_grid=smooth_factors_list,
                                            alpha_grid=alphas_list,
                                            X_data_training=X_train,
                                            y_training=y_train,
                                            max_backfit_iter=max_backfit_iter,
                                            num_jobs_factor=num_jobs_factor,
                                            cv_folds=cv_folds)
    stable_feat_indices,pred_f_3d,y_2d_training,y_2d_test=\
    gridsearch.run_stability_selection(test_size_ratio,num_shuffles,stab_gap)

    num_bins_hist=15
    print("Visualization of results")
    plots_instance=graphics.BuildPlots(num_shuffles,
                                       test_size_ratio,
                                       feat_names_full,
                                       stable_feat_indices,
                                       y_train,X_train,num_bins_hist)
    plots_instance.plot_stability_selection_results(gridsearch.stability_values,
                                                    y_2d_training,y_2d_test,
                                                    gridsearch.diagnostics_ss_array)
    plots_instance.plot_splines(pred_f_3d,gridsearch.best_parameters_cv_array)

    print("Done!")

if __name__ == "__main__":
    __spec__ = None

    parser = ap.ArgumentParser(description='Run SpAM model')
    parser.add_argument('train_data_name',
                        action="store",
                        type=str,
                        help="Name of the training data tab-delimited file")
    parser.add_argument('train_label_name',
                        action="store",
                        type=str,
                        help="Name of the training data label tab-delimited file")
    parser.add_argument('features_file_name',
                        action="store",
                        type=str,
                        help="Name of the features tab-delimited file")
    parser.add_argument('examples_file_name',
                        action="store",
                        type=str,
                        help="Name of the examples tab-delimited file")
    parser.add_argument('-test_size_ratio',
                        action="store",
                        type=float,
                        default=0.2,
                        help="Ratio of a test dataset, must be in (0,1] interval, default is 0.2")
    parser.add_argument('-stab_gap',
                        action="store",
                        type=float,
                        default=0.2,
                        help="Stability selection gap threshold, must be in (0,1] interval, default is 0.2")
    parser.add_argument('-num_shuffles',
                        action="store",
                        type=int,
                        default=10,
                        help="Number of random stratified shuffles, default is 10")
    parser.add_argument('-cv_folds',
                        action="store",
                        type=int,
                        default=5,
                        help="Number of random cross-validation shuffles, default is 5")
    parser.add_argument('-alpha_bounds',
                        action="store",
                        nargs='+',
                        type=float,
                        default=[0.1,0.9],
                        help="Range boundary values of alpha hyperparameter, default are [0.1,0.9]")
    parser.add_argument('-smoothess_power_bounds',
                        action="store",
                        nargs='+',
                        type=int,
                        default=[0,8],
                        help="Range boundary values of spline smoothnes power values, default are [0,8]")
    parser.add_argument('-num_alphas',
                        action="store",
                        type=int,
                        default=5,
                        help="Number of alphas in alpha list, default is 5")
    parser.add_argument('-smoothness_base',
                        action="store",
                        type=int,
                        default=2,
                        help="Power base for spline smoothness values, default is 2")
    parser.add_argument('-num_feat_filtered',
                        action="store",
                        type=int,
                        default=100,
                        help="Number of preselected features, default is 100")
    parser.add_argument('-data_contamination',
                        action="store",
                        type=float,
                        default=0.1,
                        help="Ratio of outliers in data, default is 0.1")

    args = parser.parse_args()
    X_train_file_name=args.train_data_name
    y_train_file_name=args.train_label_name
    features_file_name=args.features_file_name
    examples_file_name=args.examples_file_name
    test_size_ratio=args.test_size_ratio
    stab_gap=args.stab_gap
    num_shuffles=args.num_shuffles
    cv_folds=args.cv_folds
    alpha_bounds=args.alpha_bounds
    num_alphas=args.num_alphas
    smoothess_power_bounds=args.smoothess_power_bounds
    smoothness_base=args.smoothness_base
    num_feat_filtered=args.num_feat_filtered
    data_contamination=args.data_contamination

    input_data_folder_path='dataset/input_data/'
    X_train_pathway=input_data_folder_path+X_train_file_name
    y_train_pathway=input_data_folder_path+y_train_file_name
    feat_names_pathway=input_data_folder_path+features_file_name
    example_names_pathway=input_data_folder_path+examples_file_name

    [alpha_start,alpha_end]=alpha_bounds
    [smoothness_power_start,smoothness_power_end]=smoothess_power_bounds

    main(test_size_ratio,num_shuffles,stab_gap,
         X_train_pathway,y_train_pathway,
         feat_names_pathway,example_names_pathway,
         cv_folds,
         alpha_start,alpha_end,num_alphas,
         smoothness_power_start,smoothness_power_end,
         smoothness_base,num_feat_filtered,
         data_contamination)

max_stars_repo_name: imansultan/dim

max_issues_repo_name: imansultan/dim

max_forks_repo_name: imansultan/dim


text:
# Copyright (C) 2020 Intel Corporation
#
# SPDX-License-Identifier: MIT

from tools.test import *
import os
from copy import copy
import jsonpickle
import numpy as np

class ModelHandler:
    def __init__(self):
        # Setup device
        self.device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')
        torch.backends.cudnn.benchmark = True

        base_dir = os.path.abspath(os.environ.get("MODEL_PATH",
            "/opt/nuclio/SiamMask/experiments/siammask_sharp"))
        class configPath:
            config = os.path.join(base_dir, "config_davis.json")

        self.config = load_config(configPath)
        from custom import Custom
        siammask = Custom(anchors=self.config['anchors'])
        self.siammask = load_pretrain(siammask, os.path.join(base_dir, "SiamMask_DAVIS.pth"))
        self.siammask.eval().to(self.device)

    def encode_state(self, state):
        state['net.zf'] = state['net'].zf
        state.pop('net', None)
        state.pop('mask', None)

        for k,v in state.items():
            state[k] = jsonpickle.encode(v)

        return state

    def decode_state(self, state):
        for k,v in state.items():
            state[k] = jsonpickle.decode(v)

        state['net'] = copy(self.siammask)
        state['net'].zf = state['net.zf']
        del state['net.zf']

        return state

    def infer(self, image, shape, state):
        image = np.array(image)
        if state is None: # init tracking
            # xtl, ytl, xbr, ybr = shape
            xtl, ytl, xbr, ybr = [0,0,1,1]
            target_pos = np.array([(xtl + xbr) / 2, (ytl + ybr) / 2])
            target_sz = np.array([xbr - xtl, ybr - ytl])
            siammask = copy(self.siammask) # don't modify self.siammask
            state = siamese_init(image, target_pos, target_sz, siammask,
                self.config['hp'], device=self.device)
            state = self.encode_state(state)
        else: # track
            state = self.decode_state(state)
            state = siamese_track(state, image, mask_enable=True,
                refine_enable=True, device=self.device)
            shape = state['ploygon'].flatten().tolist()
            state = self.encode_state(state)

        return {"shape": shape, "state": state}

max_stars_repo_name: Linestro/happy-FaaS

max_stars_repo_licenses: Intel MIT

max_issues_repo_name: Linestro/happy-FaaS

max_issues_repo_licenses: Intel MIT

max_forks_repo_name: Linestro/happy-FaaS

max_forks_repo_licenses: Intel MIT


text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-

from argparse import Namespace
from copy import deepcopy

import numpy as np
from scipy.special import gamma, gammainc

from ..models.transforms import logsfr_ratios_to_masses
from ..sources.constants import cosmo
from .corner import quantile

__all__ = ["params_to_sfh", "parametric_pset",
           "parametric_cmf", "parametric_mwa", "parametric_sfr",
           "compute_mass_formed",
           "ratios_to_sfrs", "sfh_quantiles",
           "sfh_to_cmf", "nonpar_mwa", "nonpar_recent_sfr"]

def params_to_sfh(params, time=None, agebins=None):

    parametric = (time is not None)

    if parametric:
        taus, tages, masses = params["tau"], params["tage"], params["mass"]
        sfhs = []
        cmfs = []
        for tau, tage, mass in zip(taus, tages, masses):
            sfpar = dict(tau=tau, tage=tage, mass=mass, sfh=params["sfh"])
            sfhs.append(parametric_sfr(times=time, tavg=0, **sfpar))
            cmfs.append(parametric_cmf(times=time, **sfpar))
        lookback = time.max() - time
        sfhs = np.array(sfhs)
        cmfs = np.array(cmfs)

    else:
        logmass = params["logmass"]
        logsfr_ratios = params["logsfr_ratios"]
        sfhs = np.array([ratios_to_sfrs(logm, sr, agebins)
                        for logm, sr in zip(logmass, logsfr_ratios)])
        cmfs = sfh_to_cmf(sfhs, agebins)
        lookback = 10**(agebins-9)

    return lookback, sfhs, cmfs

def parametric_pset(logmass=None, **sfh):
    """Convert a dicionary of FSPS parametric SFH parameters into a
    namespace, making sure they are at least 1d vectors

    :param sfh: dicionary
        FSPS parameteric SFH parameters

    :returns pset:
        A Namespace instance with attributes giving the SFH parameters
    """
    # TODO: make multiple psets if any of the SFH parameters have np.size() > 1

    vectors = ["mass", "sf_start", "tage", "tau", "const", "fburst", "tburst", "sf_trunc", "sf_slope"]
    pset = Namespace(mass=1.0, sfh=4., sf_start=0, tage=1,
                     tau=1.,
                     const=0.,
                     fburst=0., tburst=1.,
                     sf_trunc=0, sf_slope=0.)

    if logmass:
        sfh["mass"] = 10**logmass
    for k in vars(pset).keys():
        if k in sfh:
            setattr(pset, k, sfh[k])
    # vectorize
    for k in vectors:
        setattr(pset, k, np.atleast_1d(getattr(pset, k)))

    return pset

def sfh_quantiles(tvec, bins, sfrs, weights=None, q=[16, 50, 84]):
    """Compute quantiles of a binned SFH

    :param tvec: shape (ntime,)
        Vector of lookback times onto which the SFH will be interpolated.

    :param bins: shape (nsamples, nbin, 2)
        The age bins, in linear untis, same units as tvec

    :param sfrs: shape (nsamples, nbin)
        The SFR in each bin

    :returns sfh_q: shape(ntime, nq)
        The quantiles of the SFHs at each lookback time in `tvec`
    """
    tt = bins.reshape(bins.shape[0], -1)
    ss = np.array([sfrs, sfrs]).transpose(1, 2, 0).reshape(bins.shape[0], -1)
    sf = np.array([np.interp(tvec, t, s, left=0, right=0) for t, s in zip(tt, ss)])
    if weights is not None:
        qq = quantile(sf.T, q=np.array(q)/100., weights=weights)
    else:
        qq = np.percentile(sf, axis=0, q=q)
    return qq

def parametric_sfr(times=None, tavg=1e-3, tage=1, **sfh):
    """Return the SFR (Msun/yr) for the given parameters of a parametric SFH,
    optionally averaging over some timescale.

    :param times: (optional, ndarray)
        If given, a set of *lookback* times where you want to calculate the sfr,
        same units as `tau` and `tage`

    :param tavg: (optional, float, default: 1e-3)
        If non-zero, average the SFR over the last `tavg` Gyr. This can help
        capture bursts.  If zero, the instantaneous SFR will be returned.

    :param sfh: optional keywords
        FSPS parametric SFH parametrs, e.g. sfh, tage, tau, sf_trunc

    :returns sfr:
        SFR in M_sun/year either for the lookback times given by `times` or at
        lookback time 0 if no times are given.  The SFR will either be
        instaneous or averaged over the last `tavg` Gyr.
    """
    if times is None:
        times = np.array(tage)

    pset = parametric_pset(tage=tage, **sfh)
    sfr, mass = compute_mass_formed(tage - times, pset)
    if tavg > 0:
        _, meps = compute_mass_formed((tage - times) - tavg, pset)
        sfr = (mass - meps) / (tavg * 1e9)
    return sfr

def parametric_cmf(times=None, tage=1., **sfh):
    """Return the cumulative formed mass for the given parameters of a
    parametric SFH.

    :param times: (optional, ndarray)
        If given, a set of *lookback* times (relative to `tage`) where you want
        to calculate the formed mass, in Gyr.  If not given, the formed mass
        will be computed for loockback time of 0.

    :param sfh: optional keywords
        FSPS parametric SFH parametrs, e.g. sfh, tage, tau, sf_trunc

    :returns mass: (ndarray)
        Mass formed up to the supplied lookback time, in units of M_sun.
        Same shape as `times`
    """
    if times is None:
        times = np.array(sfh["tage"])

    pset = parametric_pset(tage=tage**sfh)
    _, mass = compute_mass_formed(tage - times, pset)
    return mass

def parametric_mwa_numerical(tau=4, tage=13.7, power=1, n=1000):
    """Compute Mass-weighted age

    :param power: (optional, default: 1)
        Use 0 for exponential decline, and 1 for te^{-t} (delayed exponential decline)
    """
    p = power + 1
    t = np.linspace(0, tage, n)
    tavg = np.trapz((t**p)*np.exp(-t/tau), t) / np.trapz(t**(power) * np.exp(-t/tau), t)
    return tage - tavg

def parametric_mwa(tau=4, tage=13.7, power=1):
    """Compute Mass-weighted age. This is done analytically

    :param power: (optional, default: 1)
        Use 0 for exponential decline, and 1 for te^{-t} (delayed exponential decline)
    """
    tt = tage / tau
    mwt = gammainc(power+2, tt) * gamma(power+2) / gammainc(power+1, tt) * tau
    return tage - mwt

def ratios_to_sfrs(logmass, logsfr_ratios, agebins):
    """scalar
    """
    masses = logsfr_ratios_to_masses(np.squeeze(logmass),
                                     np.squeeze(logsfr_ratios),
                                     agebins)
    dt = (10**agebins[:, 1] - 10**agebins[:, 0])
    sfrs = masses / dt
    return sfrs

def nonpar_recent_sfr(logmass, logsfr_ratios, agebins, sfr_period=0.1):
    """vectorized
    """
    masses = [logsfr_ratios_to_masses(np.squeeze(logm), np.squeeze(sr), agebins)
              for logm, sr in zip(logmass, logsfr_ratios)]
    masses = np.array(masses)
    ages = 10**(agebins - 9)
    # fractional coverage of the bin by the sfr period
    ft = np.clip((sfr_period - ages[:, 0]) / (ages[:, 1] - ages[:, 0]), 0., 1)
    mformed = (ft * masses).sum(axis=-1)
    return mformed / (sfr_period * 1e9)

def nonpar_mwa(logmass, logsfr_ratios, agebins):
    """mass-weighted age, vectorized
    """
    sfrs = np.array([ratios_to_sfrs(logm, sr, agebins)
                     for logm, sr in zip(logmass, logsfr_ratios)])
    ages = 10**(agebins)
    dtsq = (ages[:, 1]**2 - ages[:, 0]**2) / 2
    mwa = [(dtsq * sfr).sum() / 10**logm
           for sfr, logm in zip(sfrs, logmass)]
    return np.array(mwa) / 1e9

def sfh_to_cmf(sfrs, agebins):
    sfrs = np.atleast_2d(sfrs)
    dt = (10**agebins[:, 1] - 10**agebins[:, 0])
    masses = (sfrs * dt)[..., ::-1]
    cmfs = masses.cumsum(axis=-1)
    cmfs /= cmfs[..., -1][..., None]
    zshape = list(cmfs.shape[:-1]) + [1]
    zeros = np.zeros(zshape)
    cmfs = np.append(zeros, cmfs, axis=-1)
    ages = 10**(np.array(agebins) - 9)
    ages = np.array(ages[:, 0].tolist() + [ages[-1, 1]])
    return ages, np.squeeze(cmfs[..., ::-1])

def compute_mass_formed(times, pset):
    """Compute the SFR and stellar mass formed in a parametric SFH,
    as a function of (forward) time.

    The linear portion of the Simha SFH (sfh=5) is defined as:
        psi(t) = psi_trunc + psi_trunc * sf_slope * (t - sf_trunc)
    where psi_trunc is the SFR of the delay-tau SFH at time sf_trunc

    :param times: ndarray of shape (nt,)
        Forward time in Gyr. Use times = pset.tage - t_lookback to
        convert from lookback times

    :param pset: Namespace instance
        The FSPS SFH parameters, assumed to be scalar or 1 element 1-d arrays.
        Usually the output of parametric_pset()

    :returns sfr: ndarray of shape (nt,)
        The instaneous SFR in M_sun/yr at each of `times`

    :returns mfromed: ndarray of shape (nt,)
        The total stellar mass formed from t=0 to `times`, in unist of M_sun
    """
    # TODO: use broadcasting to deal with multiple sfhs?

    # subtract sf_start
    tmass = pset.tage - pset.sf_start  # the age at which the mass is *specified*
    tprime = times - pset.sf_start     # the ages at which sfr and formed mass are requested

    if pset.sfh == 3:
        raise NotImplementedError("This method does not support tabular SFH")

    if tmass < 0:
        raise ValueError("SF never started (tage - sf_start < 0)")

    if (pset.const + pset.fburst) > 1:
        raise ValueError("Constant and burst fractions combine to be > 1")

    if (pset.sfh == 0):
        # SSPs
        mfrac = 1.0 * (tprime > tmass)
        sfr = np.zeros_like(tprime)  # actually the sfr is infinity

    elif pset.sfh > 0:
        # Compute tau model component, for SFH=1,4,5
        #
        # Integration limits are from 0 to tmax and 0 to tprime, where
        #   - tmass is the tage, and
        #   - tprime is the given `time`,
        #   - ttrunc is where the delay-tau truncates
        ttrunc, tend = np.max(tprime), tmass
        if (pset.sf_trunc > 0) and (pset.sf_trunc > pset.sf_start):
            # Otherwise we integrate tau model to sf_trunc - sf_start
            ttrunc = pset.sf_trunc - pset.sf_start
            tend = min(tmass, ttrunc)

        # Now integrate to get mass formed by Tprime and by Tmax, dealing with
        # truncation that happens after sf_start but before Tmax and/or Tprime.
        power = 1 + int(pset.sfh > 3)
        total_mass_tau = pset.tau * gammainc(power, tend / pset.tau)

        tt = np.clip(tprime, 0, ttrunc)
        mass_tau = (tprime > 0.) * pset.tau * gammainc(power, tt / pset.tau)
        # The SFR at Tprime (unnormalized)
        sfr_tau = (tprime > 0.) * (tprime <= ttrunc) * (tprime / pset.tau)**(power-1.) * np.exp(-tprime / pset.tau)
        # fraction of tau component mass formed by tprime
        mfrac_tau = mass_tau / total_mass_tau

    # Add the constant and burst portions, for SFH=1,4.
    if ((pset.sfh == 1) or (pset.sfh == 4)):
        # Fraction of the burst mass formed by Tprime
        tburst = (pset.tburst - pset.sf_start)
        mfrac_burst = 1.0 * (tprime > tburst)
        # SFR from constant portion at Tprime (integrates to 1 at tmax)
        sfr_const = (tprime > 0) * 1.0 / tmass
        # fraction of constant mass formed by tprime
        mfrac_const = np.clip(tprime, 0, ttrunc) * sfr_const

        # Add formed mass fractions for each component, weighted by component fractions.
        # Fraction of the constant mass formed by Tprime is just Tprime/Tmax
        # TODO : The FSPS source does not include the tburst < tmass logic....
        mfrac = ((1. - pset.const - pset.fburst * (tburst < tmass)) * mfrac_tau +
                 pset.const * mfrac_const +
                 pset.fburst * mfrac_burst)

        # N.B. for Tprime = tburst, sfr is infinite, but we ignore that case.
        sfr = ((1. - pset.const - pset.fburst) * sfr_tau / total_mass_tau + pset.const * sfr_const)
        # We've truncated
        sfr *= (tprime <= ttrunc)

    # Add the linear portion, for Simha, SFH=5.
    # This is the integral of sfr_trunc*(1 - m * (T - Ttrunc)) from Ttrunc to Tz
    elif (pset.sfh == 5):
        #raise NotImplementedError

        m = -pset.sf_slope
        if (m > 0):
            # find time at which SFR=0, if m>0
            Tz = ttrunc + 1.0 / m
        else:
            # m <= 0 will never reach SFR=0
            Tz = np.max(tprime)

        # Logic for Linear portion
        if (ttrunc < 0):
            # Truncation does not occur during the SFH.
            total_mass_linear = 0.
            mass_linear = 0.
            sfr = sfr_tau / total_mass_tau
        else:
            # Truncation does occur, integrate linear to zero crossing or tage.
            Thi = min(Tz, tmass)
            sfr_trunc = (ttrunc/pset.tau) * np.exp(-ttrunc / pset.tau)
            total_mass_linear = (Thi > ttrunc) * sfr_trunc * linear_mass(Thi, ttrunc, m)
            mass_linear = (tprime > ttrunc) * sfr_trunc * linear_mass(tprime, ttrunc, m)
            mass_linear[tprime > Tz] = sfr_trunc * linear_mass(Tz, ttrunc, m)
            # SFR in linear portion
            sfr = sfr_trunc * (1 - m * (tprime - ttrunc)) / (total_mass_tau + total_mass_linear)
            sfr *= ((tprime > ttrunc) & (tprime <= Tz))
            # add portion for tau
            sfr[tprime <= ttrunc] = sfr_tau[tprime <= ttrunc] / (total_mass_tau + total_mass_linear)

        mfrac = (mass_tau + mass_linear) / (total_mass_tau + total_mass_linear)

    return pset.mass * sfr/1e9, pset.mass * mfrac

def linear_mass(t, ttrunc, m):
    """Integrate (1-m*(a-ttrunc)) da from a=ttrunc to a=t
    """
    tt = t - ttrunc
    return ((tt + ttrunc * m * tt) - m/2. * (t**2 - ttrunc**2))

default_sfh = dict(mass=1.0, sfh=4, tage=1., tau=2.,
                   sf_start=0., sf_trunc=0., fburst=0., const=0., tburst=0.5)

def show_par_sfh(times, label="", axes=[], tavg=0.01, tol=1e-3, **params):
    sfh = deepcopy(default_sfh)
    sfh.update(params)
    pset = parametric_pset(**sfh)
    sfr, mass = compute_mass_formed(times, pset)
    sf_label = r"SFR"
    if tavg > 0:
        sfr = parametric_sfr(times, **sfh)
        sf_label = r"$\langle {{\rm SFR}}\rangle_{{{}}}$".format(tavg)

    ax = axes[0]
    ax.plot(times, sfr)
    ax.set_ylabel(sf_label)

    ax = axes[1]
    ax.plot(times, mass)
    ax.set_ylabel("M(<t)")
    ax.axhline(pset.mass, linestyle=":", color='k', linewidth=1.0, label="mass at tage")

    [ax.axvline(pset.tage, linestyle="--", color='firebrick', label="tage") for ax in axes]
    if pset.sf_trunc > 0:
        [ax.axvline(pset.sf_trunc, linestyle="--", color='seagreen', label="sf_trunc") for ax in axes]
    if pset.fburst > 0:
        [ax.axvline(pset.tburst, linestyle="--", color='darkorange', label="tburst") for ax in axes]
    if pset.sfh > 5:
        [ax.axvline(pset.sf_trunc, linestyle="--", label="tburst") for ax in axes]

    mr = np.trapz(sfr, times) * 1e9 / mass.max()
    print("{}: {:0.6f}".format(label, mr))
    if (pset.sfh > 0) & (pset.fburst == 0):
        assert np.abs(mr - 1) < tol
    return mr

if __name__ == "__main__":

    # time array
    times = np.linspace(0, 5, 1000)

    # --- test cases ---
    ncases = 14
    import matplotlib.pyplot as pl
    pl.ion()
    fig, axes = pl.subplots(ncases + 1, 2, sharex="row",
                            figsize=(8.5, ncases * 1.5))
    i = 0

    # default
    show_par_sfh(times, label="default", axes=axes[i])

    # sf start
    i += 1
    show_par_sfh(times, label="test SF start", axes=axes[i],
                 sf_start=1.0, tage=2.0)

    # with const
    i += 1
    show_par_sfh(times, label="with const", axes=axes[i],
                 const=0.5)

    # const w/ sf_start
    i += 1
    show_par_sfh(times, label="const w/ sf_start", axes=axes[i],
                 const=0.5, sf_start=1.0, tage=2.0)

    # const w/ sf_start & trunc
    i += 1
    show_par_sfh(times, label="const w/ sf_start & trunc", axes=axes[i],
                 const=0.5, sf_trunc=4., sf_start=1.0, tage=2.0)

    # pure const w/ sf_start & trunc
    i += 1
    show_par_sfh(times, label="pure const w/ sf_start & trunc", tol=5e-3, axes=axes[i],
                 const=1., sf_trunc=4., sf_start=1.0, tage=2.0)

    # burst before tage
    i += 1
    show_par_sfh(times, label="burst before tage", axes=axes[i],
                 fburst=0.5, tburst=0.5)

    # burst after tage
    i += 1
    show_par_sfh(times, label="burst after tage", axes=axes[i],
                 fburst=0.5, tburst=1.5)

    # burst at sf_start
    i += 1
    show_par_sfh(times, label="burst at sf_start", axes=axes[i],
                 fburst=0.5, tburst=1.5, sf_start=1.5, tage=3)

    # pure exp
    i += 1
    show_par_sfh(times, label="pure EXP", tol=1e-2, axes=axes[i],
                 sfh=1)

    # SSP
    i += 1
    show_par_sfh(times, label="SSP", axes=axes[i],
                 sfh=0)

    # positive slope sf trunc before tage
    show_par_sfh(times, label="pos quench before tage", axes=axes[-4],
                 sfh=5, tage=3, sf_trunc=2, tau=0.5, sf_slope=1)

    # negative slope sf trunc before tage
    show_par_sfh(times, label="neg quench before tage", axes=axes[-3],
                 sfh=5, tage=2, sf_trunc=1.5, tau=3, sf_slope=-1)

    # positive slope sf trunc after tage
    show_par_sfh(times, label="pos quench after tage", axes=axes[-2],
                 sfh=5, sf_trunc=3, tau=1, sf_slope=1)

    # negative slope sf trunc after tage
    show_par_sfh(times, label="neg quench after tage", axes=axes[-1],
                 sfh=5, sf_trunc=3, tau=3, sf_slope=-1)

max_stars_repo_name: wpb-astro/prospector

max_issues_repo_name: wpb-astro/prospector

max_forks_repo_name: wpb-astro/prospector

reason: import numpy,from scipy


text:
import gzip
import numpy as np

def vectorized_result(j):
    """
    Return a (10, 1)-dimensional unit vector with a 1.0 in the jth position and
    zeroes elsewhere. This is used to convert a digit (0...9) into a corresponding
    desired output from the neural network.
    """

    e = np.zeros((10, 1))
    e[j] = 1.0
    return e

def load_data():
    """
    Return a tuple containing `(training_data, validation_data, test_data)`.

    In particular, `training_data` is a list containing 50'000 2-tuples `(x, y)`.
    `x` is a (784, 1)-dimensional numpy.ndarray containing the input image. `y` is a
    (10, 1)-dimensional numpy.ndarray representing the unit vector corresponding
    to the correct digit for `x`.

    `validation_data` and `test_data` are lists containing 10'000 items of the
    same format as `training_data`.
    """

    with gzip.open('../data/train-images-idx3-ubyte.gz', 'rb') as f:
        train_images = f.read()
    with gzip.open('../data/train-labels-idx1-ubyte.gz', 'rb') as f:
        train_labels = f.read()
    with gzip.open('../data/test-images-idx3-ubyte.gz', 'rb') as f:
        test_images = f.read()
    with gzip.open('../data/test-labels-idx1-ubyte.gz', 'rb') as f:
        test_labels = f.read()

    train_images = train_images[16:]
    training_inputs = [
        np.frombuffer(train_images[i:i+784], np.uint8).reshape((784, 1)) / 255
            for i in map(lambda x: 784*x, range(50000))
    ]
    validation_inputs = [
        np.frombuffer(train_images[i:i+784], np.uint8).reshape((784, 1)) / 255
            for i in map(lambda x: 784*x, range(50000, 60000))
    ]

    train_labels = train_labels[8:]
    training_results = [
        vectorized_result(train_labels[i])
            for i in range(50000)
    ]
    validation_results = [
        vectorized_result(train_labels[i])
            for i in range(50000, 60000)
    ]

    test_images = test_images[16:]
    test_inputs = [
        np.frombuffer(test_images[i:i+784], np.uint8).reshape((784, 1)) / 255
            for i in map(lambda x: 784*x, range(10000))
    ]

    test_labels = test_labels[8:]
    test_results = [
        vectorized_result(test_labels[i])
            for i in range(10000)
    ]

    training_data = list(zip(training_inputs, training_results))
    validation_data = list(zip(validation_inputs, validation_results))
    test_data = list(zip(test_inputs, test_results))

    return training_data, validation_data, test_data


text:
import biom
import networkx as nx
import os
from mako.scripts.neo4biom import Biom2Neo
from mako.scripts.io import IoDriver

def main():
    filepath = "C:/Users/username/demo/"

    biom_filepaths = [filepath + "11766.biom",
                      filepath + "11888.biom",
                      filepath + "11947.biom",
                      filepath + "12021.biom",
                      filepath + "12716.biom"]

    network_filepaths = [filepath + "11766.graphml",
                         filepath + "11888.graphml",
                         filepath + "11947.graphml",
                         filepath + "12021.graphml",
                         filepath + "12716.graphml"]

    driver = Biom2Neo(uri='neo4j://localhost:7688',
                      user='neo4j',
                      password='test',
                      filepath=filepath,
                      encrypted=False)

    for file in biom_filepaths:
        name = os.path.basename(file)
        name = name.split(".")[0]
        biomtab = biom.load_table(file)
        driver.convert_biom(biomfile=biomtab, exp_id=name, obs=True)
    driver.close()

    driver = IoDriver(uri='neo4j://localhost:7688',
                      user='neo4j',
                      password='test',
                      filepath=filepath,
                      encrypted=False)

    for file in network_filepaths:
        name = os.path.basename(file)
        name = name.split(".")[0]
        net = nx.read_graphml(file)
        driver.convert_networkx(network=net, network_id=name)
    driver.close()

if __name__ == '__main__':
    main()
    print("Completed running database setup.")

reason: import networkx


text:
import numpy as np

from .points import get_fractions, interpolate_along_axis

def prepare(start: np.ndarray, stop: np.ndarray):
    start = normalize(start)
    stop = normalize(stop)

    diff = len(start) - len(stop)
    if diff < 0:
        start = add_points(start, -diff)
    elif diff > 0:
        stop = add_points(stop, diff)

    start = rotate(start, stop)
    return start, stop

def get_area(line):
    line = np.insert(line, 0, line[-1], axis=0)
    return np.cross(line[1:], line[:-1]).sum() / 2

def normalize(line):
    if line.shape[1] != 2:
        raise ValueError('In case of closed lines only 2D points are supported.')

    if get_area(line) > 0:
        line = line[::-1]

    return line

def rotate(line, reference):
    min_, best_offset = np.inf, 0

    for offset in range(len(line)):
        deviation = np.linalg.norm(reference - np.roll(line, offset, axis=0))

        if deviation < min_:
            min_ = deviation
            best_offset = offset

    if best_offset:
        line = np.roll(line, best_offset, axis=0)

    return line

def add_points(line, n_points):
    line = np.insert(line, len(line), line[0], axis=0)
    fractions = get_fractions(line)
    new_fractions = sorted(
        list(fractions) +
        list(np.linspace(0, 1, n_points, endpoint=False) + 1 / (2 * n_points))
    )[:-1]

    return interpolate_along_axis(new_fractions, fractions, line)