text:
"""
A few different backtracking line search subsolvers.

BoundsEnforceLS - Only checks bounds and enforces them by one of three methods.
ArmijoGoldsteinLS -- Like above, but terminates with the ArmijoGoldsteinLS condition.

"""

import sys
import numpy as np

from openmdao.core.analysis_error import AnalysisError
from openmdao.solvers.solver import NonlinearSolver
from openmdao.recorders.recording_iteration_stack import Recording

def _print_violations(outputs, lower, upper):
    """
    Print out which variables exceed their bounds.

    Parameters
    ----------
    outputs : <Vector>
        Vector containing the outputs.
    lower : <Vector>
        Vector containing the lower bounds.
    upper : <Vector>
        Vector containing the upper bounds.
    """
    start = end = 0
    for name, val in outputs._abs_item_iter():
        end += val.size
        if upper is not None and any(val > upper[start:end]):
            print("'%s' exceeds upper bounds" % name)
            print("  Val:", val)
            print("  Upper:", upper[start:end], '\n')

        if lower is not None and any(val < lower[start:end]):
            print("'%s' exceeds lower bounds" % name)
            print("  Val:", val)
            print("  Lower:", lower[start:end], '\n')

        start = end

class LinesearchSolver(NonlinearSolver):
    """
    Base class for line search solvers.

    Attributes
    ----------
    _do_subsolve : bool
        Flag used by parent solver to tell the line search whether to solve subsystems while
        backtracking.
    _lower_bounds : ndarray or None
        Lower bounds array.
    _upper_bounds : ndarray or None
        Upper bounds array.
    """

    def __init__(self, **kwargs):
        """
        Initialize all attributes.

        Parameters
        ----------
        **kwargs : dict
            Options dictionary.
        """
        super().__init__(**kwargs)
        # Parent solver sets this to control whether to solve subsystems.
        self._do_subsolve = False
        self._lower_bounds = None
        self._upper_bounds = None

    def _declare_options(self):
        """
        Declare options before kwargs are processed in the init method.
        """
        super()._declare_options()
        opt = self.options
        opt.declare(
            'bound_enforcement', default='scalar', values=['vector', 'scalar', 'wall'],
            desc="If this is set to 'vector', the entire vector is backtracked together " +
                 "when a bound is violated. If this is set to 'scalar', only the violating " +
                 "entries are set to the bound and then the backtracking occurs on the vector " +
                 "as a whole. If this is set to 'wall', only the violating entries are set " +
                 "to the bound, and then the backtracking follows the wall - i.e., the " +
                 "violating entries do not change during the line search.")
        opt.declare('print_bound_enforce', default=False,
                    desc="Set to True to print out names and values of variables that are pulled "
                    "back to their bounds.")

    def _setup_solvers(self, system, depth):
        """
        Assign system instance, set depth, and optionally perform setup.

        Parameters
        ----------
        system : System
            pointer to the owning system.
        depth : int
            depth of the current system (already incremented).
        """
        super()._setup_solvers(system, depth)
        if system._has_bounds:
            abs2meta_out = system._var_abs2meta['output']
            start = end = 0
            for abs_name, val in system._outputs._abs_item_iter():
                end += val.size
                meta = abs2meta_out[abs_name]
                var_lower = meta['lower']
                var_upper = meta['upper']

                if var_lower is None and var_upper is None:
                    start = end
                    continue

                ref0 = meta['ref0']
                ref = meta['ref']

                if not np.isscalar(ref0):
                    ref0 = ref0.ravel()
                if not np.isscalar(ref):
                    ref = ref.ravel()

                if var_lower is not None:
                    if self._lower_bounds is None:
                        self._lower_bounds = np.full(len(system._outputs), -np.inf)
                    if not np.isscalar(var_lower):
                        var_lower = var_lower.ravel()
                    self._lower_bounds[start:end] = (var_lower - ref0) / (ref - ref0)

                if var_upper is not None:
                    if self._upper_bounds is None:
                        self._upper_bounds = np.full(len(system._outputs), np.inf)
                    if not np.isscalar(var_upper):
                        var_upper = var_upper.ravel()
                    self._upper_bounds[start:end] = (var_upper - ref0) / (ref - ref0)

                start = end
        else:
            self._lower_bounds = self._upper_bounds = None

    def _enforce_bounds(self, step, alpha):
        """
        Enforce lower/upper bounds.

        Modifies the vector of outputs and the step.

        Parameters
        ----------
        step : <Vector>
            Newton step; the backtracking is applied to this vector in-place.
        alpha : float
            Step size parameter.
        """
        system = self._system()
        if not system._has_bounds:
            return

        options = self.options
        method = options['bound_enforcement']
        lower = self._lower_bounds
        upper = self._upper_bounds

        if options['print_bound_enforce']:
            _print_violations(system._outputs, lower, upper)

        if method == 'vector':
            _enforce_bounds_vector(system._outputs, step, alpha, lower, upper)
        elif method == 'scalar':
            _enforce_bounds_scalar(system._outputs, step, alpha, lower, upper)
        elif method == 'wall':
            _enforce_bounds_wall(system._outputs, step, alpha, lower, upper)

class BoundsEnforceLS(LinesearchSolver):
    """
    Bounds enforcement only.

    Not so much a linesearch; just check the bounds and if they are violated, then pull back to a
    non-violating point and evaluate.
    """

    SOLVER = 'LS: BCHK'

    def _declare_options(self):
        """
        Declare options before kwargs are processed in the init method.
        """
        super()._declare_options()
        opt = self.options

        # Remove unused options from base options here, so that users
        # attempting to set them will get KeyErrors.
        for unused_option in ("atol", "rtol", "maxiter", "err_on_non_converge"):
            opt.undeclare(unused_option)

    def _solve(self):
        """
        Run the iterative solver.
        """
        self._iter_count = 0
        system = self._system()

        u = system._outputs
        du = system._vectors['output']['linear']

        if not system._has_bounds:
            u += du
            return

        self._run_apply()

        norm0 = self._iter_get_norm()
        if norm0 == 0.0:
            norm0 = 1.0
        self._norm0 = norm0
        u += du

        with Recording('BoundsEnforceLS', self._iter_count, self) as rec:
            self._enforce_bounds(step=du, alpha=1.0)

            self._run_apply()
            norm = self._iter_get_norm()

            # Save the norm values in the context manager so they can also be recorded.
            rec.abs = norm
            rec.rel = norm / norm0

        self._mpi_print(self._iter_count, norm, norm / norm0)

class ArmijoGoldsteinLS(LinesearchSolver):
    """
    Backtracking line search that terminates using the Armijo-Goldstein condition.

    Attributes
    ----------
    _analysis_error_raised : bool
        Flag is set to True if a subsystem raises an AnalysisError.
    """

    SOLVER = 'LS: AG'

    def __init__(self, **kwargs):
        """
        Initialize all attributes.

        Parameters
        ----------
        **kwargs : dict
            Options dictionary.
        """
        super().__init__(**kwargs)

        self._analysis_error_raised = False

    def _line_search_objective(self):
        """
        Calculate the objective function of the line search.

        Returns
        -------
        float
            Line search objective (residual norm).
        """
        return self._iter_get_norm()

    def _iter_initialize(self):
        """
        Perform any necessary pre-processing operations.

        Returns
        -------
        float
            initial error.
        float
            error at the first iteration.
        """
        system = self._system()
        self.alpha = alpha = self.options['alpha']

        u = system._outputs
        du = system._vectors['output']['linear']

        self._run_apply()
        phi0 = self._line_search_objective()
        if phi0 == 0.0:
            phi0 = 1.0
        self._phi0 = phi0
        # From definition of Newton's method one full step should drive the linearized residuals
        # to 0, hence the directional derivative is equal to the initial function value.
        self._dir_derivative = -phi0

        # Initial step length based on the input step length parameter
        u.add_scal_vec(alpha, du)

        self._enforce_bounds(step=du, alpha=alpha)

        try:
            cache = self._solver_info.save_cache()

            self._run_apply()
            phi = self._line_search_objective()

        except AnalysisError as err:
            self._solver_info.restore_cache(cache)

            if self.options['retry_on_analysis_error']:
                self._analysis_error_raised = True
            else:
                raise err

            phi = np.nan

        return phi

    def _declare_options(self):
        """
        Declare options before kwargs are processed in the init method.
        """
        super()._declare_options()
        opt = self.options
        opt['maxiter'] = 5
        opt.declare('c', default=0.1, lower=0.0, upper=1.0, desc="Slope parameter for line of "
                    "sufficient decrease. The larger the step, the more decrease is required to "
                    "terminate the line search.")
        opt.declare('rho', default=0.5, lower=0.0, upper=1.0, desc="Contraction factor.")
        opt.declare('alpha', default=1.0, lower=0.0, desc="Initial line search step.")
        opt.declare('retry_on_analysis_error', default=True,
                    desc="Backtrack and retry if an AnalysisError is raised.")
        opt.declare('method', default='Armijo', values=['Armijo', 'Goldstein'],
                    desc="Method to calculate stopping condition.")

    def _single_iteration(self):
        """
        Perform the operations in the iteration loop.
        """
        self._analysis_error_raised = False
        system = self._system()

        # Hybrid newton support.
        if self._do_subsolve and self._iter_count > 0:
            self._solver_info.append_solver()

            try:
                cache = self._solver_info.save_cache()
                self._gs_iter()
                self._run_apply()

            except AnalysisError as err:
                self._solver_info.restore_cache(cache)

                if self.options['retry_on_analysis_error']:
                    self._analysis_error_raised = True

                else:
                    raise err

            finally:
                self._solver_info.pop()

        else:
            self._run_apply()

    def _stopping_criteria(self, fval, method):
        """
        Sufficient decrease criteria for the line search.

        The initial line search objective and the step length parameter are stored in the class
        instance.

        Parameters
        ----------
        fval : float
            Current line search objective value.
        method : str, optional
            Method to calculate stopping condition. Can be "Armijo" or "Goldstein".

        Returns
        -------
        bool
            Stopping condition is satisfied.
        """
        method = method.lower()
        fval0 = self._phi0
        df_dalpha = self._dir_derivative
        c1 = self.options['c']
        alpha = self.alpha
        if method == 'armijo':
            return fval <= fval0 + c1 * alpha * df_dalpha
        elif method == 'goldstein':
            return fval0 + (1 - c1) * alpha * df_dalpha <= fval <= fval0 + c1 * alpha * df_dalpha

    def _update_step_length_parameter(self, rho):
        """
        Update the step length parameter by multiplying with the contraction factor.

        Parameters
        ----------
        rho : float
            Contraction factor
        """
        self.alpha *= rho  # update alpha

    def _solve(self):
        """
        Run the iterative solver.
        """
        options = self.options
        maxiter = options['maxiter']
        rho = options['rho']
        method = options['method']

        system = self._system()
        u = system._outputs
        du = system._vectors['output']['linear']  # Newton step

        self._iter_count = 0
        phi = self._iter_initialize()
        phi0 = self._phi0

        # Further backtracking if needed.
        while (self._iter_count < maxiter and
               (not self._stopping_criteria(phi, method) or self._analysis_error_raised)):

            with Recording('ArmijoGoldsteinLS', self._iter_count, self) as rec:

                if self._iter_count > 0:
                    alpha_old = self.alpha
                    self._update_step_length_parameter(rho)
                    # Moving on the line search with the difference of the old and new step length.
                    u.add_scal_vec(self.alpha - alpha_old, du)
                cache = self._solver_info.save_cache()

                try:
                    self._single_iteration()
                    self._iter_count += 1

                    phi = self._line_search_objective()

                    # Save the norm values in the context manager so they can also be recorded.
                    rec.abs = phi
                    rec.rel = phi / phi0

                except AnalysisError as err:
                    self._solver_info.restore_cache(cache)
                    self._iter_count += 1

                    if self.options['retry_on_analysis_error']:
                        self._analysis_error_raised = True
                        rec.abs = np.nan
                        rec.rel = np.nan

                    else:
                        raise err

            # self._mpi_print(self._iter_count, norm, norm / norm0)
            self._mpi_print(self._iter_count, phi, self.alpha)

def _enforce_bounds_vector(u, du, alpha, lower_bounds, upper_bounds):
    """
    Enforce lower/upper bounds, backtracking the entire vector together.

    This method modifies both self (u) and step (du) in-place.

    Parameters
    ----------
    u :<Vector>
        Output vector.
    du : <Vector>
        Newton step; the backtracking is applied to this vector in-place.
    alpha : float
        step size.
    lower_bounds : ndarray
        Lower bounds array.
    upper_bounds : ndarray
        Upper bounds array.
    """
    # The assumption is that alpha * du has been added to self (i.e., u)
    # just prior to this method being called. We are currently in the
    # initialization of a line search, and we're trying to ensure that
    # the u does not violate bounds in the first iteration. If it does,
    # we modify the du vector directly.

    # This is the required change in step size, relative to the du vector.
    d_alpha = 0

    # Find the largest amount a bound is violated
    # where positive means a bound is violated - i.e. the required d_alpha.
    du_arr = du.asarray()
    mask = du_arr != 0
    if mask.any():
        abs_du_mask = np.abs(du_arr[mask])
        u_mask = u.asarray()[mask]

        # Check lower bound
        if lower_bounds is not None:
            max_d_alpha = np.amax((lower_bounds[mask] - u_mask) / abs_du_mask)
            if max_d_alpha > d_alpha:
                d_alpha = max_d_alpha

        # Check upper bound
        if upper_bounds is not None:
            max_d_alpha = np.amax((u_mask - upper_bounds[mask]) / abs_du_mask)
            if max_d_alpha > d_alpha:
                d_alpha = max_d_alpha

    if d_alpha > 0:
        # d_alpha will not be negative because it was initialized to be 0
        # and we've only done max operations.
        # d_alpha will not be greater than alpha because the assumption is that
        # the original point was valid - i.e., no bounds were violated.
        # Therefore 0 <= d_alpha <= alpha.

        # We first update u to reflect the required change to du.
        u.add_scal_vec(-d_alpha, du)

        # At this point, we normalize d_alpha by alpha to figure out the relative
        # amount that the du vector has to be reduced, then apply the reduction.
        du *= 1 - d_alpha / alpha

def _enforce_bounds_scalar(u, du, alpha, lower_bounds, upper_bounds):
    """
    Enforce lower/upper bounds on each scalar separately, then backtrack as a vector.

    This method modifies both self (u) and step (du) in-place.

    Parameters
    ----------
    u :<Vector>
        Output vector.
    du : <Vector>
        Newton step; the backtracking is applied to this vector in-place.
    alpha : float
        step size.
    lower_bounds : ndarray
        Lower bounds array.
    upper_bounds : ndarray
        Upper bounds array.
    """
    # The assumption is that alpha * step has been added to this vector
    # just prior to this method being called. We are currently in the
    # initialization of a line search, and we're trying to ensure that
    # the initial step does not violate bounds. If it does, we modify
    # the step vector directly.

    # enforce bounds on step in-place.
    u_data = u.asarray()

    # If u > lower, we're just adding zero. Otherwise, we're adding
    # the step required to get up to the lower bound.
    # For du, we normalize by alpha since du eventually gets
    # multiplied by alpha.
    change_lower = 0. if lower_bounds is None else np.maximum(u_data, lower_bounds) - u_data

    # If u < upper, we're just adding zero. Otherwise, we're adding
    # the step required to get down to the upper bound, but normalized
    # by alpha since du eventually gets multiplied by alpha.
    change_upper = 0. if upper_bounds is None else np.minimum(u_data, upper_bounds) - u_data

    change = change_lower + change_upper
    u_data += change
    du += change / alpha

def _enforce_bounds_wall(u, du, alpha, lower_bounds, upper_bounds):
    """
    Enforce lower/upper bounds on each scalar separately, then backtrack along the wall.

    This method modifies both self (u) and step (du) in-place.

    Parameters
    ----------
    u :<Vector>
        Output vector.
    du : <Vector>
        Newton step; the backtracking is applied to this vector in-place.
    alpha : float
        step size.
    lower_bounds : ndarray
        Lower bounds array.
    upper_bounds : ndarray
        Upper bounds array.
    """
    # The assumption is that alpha * step has been added to this vector
    # just prior to this method being called. We are currently in the
    # initialization of a line search, and we're trying to ensure that
    # the initial step does not violate bounds. If it does, we modify
    # the step vector directly.

    # enforce bounds on step in-place.
    u_data = u.asarray()
    du_data = du.asarray()

    # If u > lower, we're just adding zero. Otherwise, we're adding
    # the step required to get up to the lower bound.
    # For du, we normalize by alpha since du eventually gets
    # multiplied by alpha.
    change_lower = 0. if lower_bounds is None else np.maximum(u_data, lower_bounds) - u_data

    # If u < upper, we're just adding zero. Otherwise, we're adding
    # the step required to get down to the upper bound, but normalized
    # by alpha since du eventually gets multiplied by alpha.
    change_upper = 0. if upper_bounds is None else np.minimum(u_data, upper_bounds) - u_data

    change = change_lower + change_upper

    u_data += change
    du_data += change / alpha

    # Now we ensure that we will backtrack along the wall during the
    # line search by setting the entries of du at the bounds to zero.
    changed_either = change.astype(bool)
    du_data[changed_either] = 0.

max_stars_repo_name: wright/OpenMDAO

max_issues_repo_name: wright/OpenMDAO

max_forks_repo_name: wright/OpenMDAO


text:
import unittest
import numpy as np
import tensorflow as tf
from timemachine import force
from timemachine import bonded_force

class PeriodicTorsionForceOpenMM(force.ConservativeForce):

    def __init__(self,
        params,
        torsion_idxs,
        param_idxs,
        precision=tf.float64):
        """
        Implements the OpenMM version of PeriodicTorsionForce
        """
        self.params = params
        self.torsion_idxs = torsion_idxs
        self.param_idxs = param_idxs

    @staticmethod
    def get_signed_angle(ci, cj, ck, cl):
        rij = ci - cj
        rkj = ck - cj
        rkl = ck - cl

        n1 = tf.cross(rij, rkj)
        n2 = tf.cross(rkj, rkl)

        lhs = tf.norm(n1, axis=-1)
        rhs = tf.norm(n2, axis=-1)
        bot = lhs * rhs

        top = tf.reduce_sum(tf.multiply(n1, n2), -1)
        cos_angles = top/bot # debug

        sign = tf.sign(tf.reduce_sum(tf.multiply(rkj, tf.cross(n1, n2)), -1))
        angle = sign * tf.acos(cos_angles)

        return angle

    def energy(self, conf):
        ci = tf.gather(conf, self.torsion_idxs[:, 0])
        cj = tf.gather(conf, self.torsion_idxs[:, 1])
        ck = tf.gather(conf, self.torsion_idxs[:, 2])
        cl = tf.gather(conf, self.torsion_idxs[:, 3])

        k0s = tf.gather(self.params, self.param_idxs[:, 0])
        k1s = tf.gather(self.params, self.param_idxs[:, 1])
        k2s = tf.gather(self.params, self.param_idxs[:, 2])
        t0s = tf.gather(self.params, self.param_idxs[:, 3])
        t1s = tf.gather(self.params, self.param_idxs[:, 4])
        t2s = tf.gather(self.params, self.param_idxs[:, 5])

        angle = self.get_signed_angle(ci, cj, ck, cl)

        e0 = k0s*(1+tf.cos(1 * angle - t0s)) # cos(a) cos(t0) + sin(a) sin(t0)
        e1 = k1s*(1+tf.cos(2 * angle - t1s)) # cos(2a) cos(t1) + sin(2a) sin(t1) -> (2*cos(a)^2 - 1) cos(t1) + (2*sin(a)*cos(a)) sin(t1s)
        e2 = k2s*(1+tf.cos(3 * angle - t2s)) # cos(3a) cos(t2) + sin(3a) sin(t2) -> (4*cos(a)^3 - 3 cos(a)) cos(t2) + (3 sin(a) - 4 * sin(a)^3) sin(t2)
        return tf.reduce_sum(e0+e1+e2, axis=-1)

class TestBondedForce(unittest.TestCase):

    def setUp(self):
        self.conformers = np.array([
            [[-0.6000563454193615, 0.376172954382274 ,-0.2487295756125901],
             [ 0.561317027011325 , 0.2066950040043141, 0.3670430960815993],
             [-1.187055522272264 ,-0.3415864358441354, 0.0871382207830652],
             [ 0.9399773448903637,-0.6888774474110431, 0.2104211949995816]],
            [[-0.6000563454193615, 0.376172954382274 ,-0.2487295756125901],
             [ 0.5613170270113252, 0.2066950040043142, 0.3670430960815993],
             [-1.187055522272264 ,-0.3415864358441354, 0.0871382207830652],
             [ 1.283345455745044 ,-0.0356257425880843,-0.2573923896494185]],
            [[-0.6000563454193615, 0.376172954382274 ,-0.2487295756125901],
             [ 0.561317027011325 , 0.2066950040043142, 0.3670430960815992],
             [-1.187055522272264 ,-0.3415864358441354, 0.0871382207830652],
             [ 1.263820400176392 , 0.7964992122869241, 0.0084568741589791]],
            [[-0.6000563454193615, 0.376172954382274 ,-0.2487295756125901],
             [ 0.5613170270113252, 0.2066950040043142, 0.3670430960815992],
             [-1.187055522272264 ,-0.3415864358441354, 0.0871382207830652],
             [ 0.8993534242298198, 1.042445571242743 , 0.7635483993060286]],
            [[-0.6000563454193615, 0.376172954382274 ,-0.2487295756125901],
             [ 0.5613170270113255, 0.2066950040043142, 0.3670430960815993],
             [-1.187055522272264 ,-0.3415864358441354, 0.0871382207830652],
             [ 0.5250337847650304, 0.476091386095139 , 1.3136545198545133]],
            [[-0.6000563454193615, 0.376172954382274 ,-0.2487295756125901],
             [ 0.5613170270113255, 0.2066950040043141, 0.3670430960815993],
             [-1.187055522272264 ,-0.3415864358441354, 0.0871382207830652],
             [ 0.485009232042489 ,-0.3818599172073237, 1.1530102055165103]],
            ], dtype=np.float64)

        self.nan_conformers = np.array([
            [[-0.6000563454193615, 0.376172954382274 ,-0.2487295756125901],
             [ 0.5613170270113252, 0.2066950040043142, 0.3670430960815993],
             [-1.187055522272264 ,-0.3415864358441354, 0.0871382207830652],
             [ 1.2278668040866427, 0.8805184219394547, 0.099391329616366 ]],
            [[-0.6000563454193615, 0.376172954382274 ,-0.2487295756125901],
             [ 0.561317027011325 , 0.206695004004314 , 0.3670430960815994],
             [-1.187055522272264 ,-0.3415864358441354, 0.0871382207830652],
             [ 0.5494071252089705,-0.5626592973923106, 0.9817919758125693]],
            ], dtype=np.float64)

    def tearDown(self):
        tf.reset_default_graph()

    def test_torsions_with_openmm(self):
        """
        Test agreement of torsions with OpenMM's implementation of torsion terms.
        """

        torsion_idxs = np.array([
            [0, 1, 2, 3],
        ], dtype=np.int32)

        params = tf.get_variable(name="parameters", dtype=tf.float64, shape=(6,), initializer=tf.constant_initializer([
            2.3, # k0
            5.4, # k1
            9.0, # k2
            0.0, # t0
            3.0, # t1
            5.8, # t2
        ]))

        param_idxs = np.array([
            [0, 1, 2, 3, 4, 5]
        ], dtype=np.int32)

        x_ph = tf.placeholder(shape=(4, 3), dtype=tf.float64)

        test_force = bonded_force.PeriodicTorsionForce(params, torsion_idxs, param_idxs)
        ref_force = PeriodicTorsionForceOpenMM(params, torsion_idxs, param_idxs)

        test_nrg = test_force.energy(x_ph)
        test_grad = test_force.gradients(x_ph)
        test_hessian = test_force.hessians(x_ph)
        test_mixed = test_force.mixed_partials(x_ph)

        ref_nrg = ref_force.energy(x_ph)
        ref_grad = ref_force.gradients(x_ph)
        ref_hessian = ref_force.hessians(x_ph)
        ref_mixed = ref_force.mixed_partials(x_ph)

        sess = tf.Session()
        sess.run(tf.initializers.global_variables())

        for conf_idx, conf in enumerate(self.conformers):

            t0, t1, t2, t3, r0, r1, r2, r3 = sess.run([
                test_nrg,
                test_grad,
                test_hessian,
                test_mixed,
                ref_nrg,
                ref_grad,
                ref_hessian,
                ref_mixed], feed_dict={x_ph: conf})

            np.testing.assert_array_almost_equal(t0, r0, decimal=14) # energy
            t1, r1 = t1.values[t1.indices], r1.values[r1.indices]
            np.testing.assert_array_almost_equal(t1, r1, decimal=13) # grad
            np.testing.assert_array_almost_equal(t2, r2, decimal=12) # hessian
            np.testing.assert_array_almost_equal(t3, r3, decimal=13) # mixed partials

        # OpenMM's vanilla implementation of the energy is non-differentiable when
        # the angle is equal to zero. The atan2 version is numerically stable since
        # it avoids taking a derivative of an arccos.
        for conf_idx, conf in enumerate(self.nan_conformers):

            t0, t1, t2, t3, r0, r1, r2, r3 = sess.run([
                test_nrg,
                test_grad,
                test_hessian,
                test_mixed,
                ref_nrg,
                ref_grad,
                ref_hessian,
                ref_mixed], feed_dict={x_ph: conf})

            np.testing.assert_array_almost_equal(t0, r0, decimal=14) # energy
            t1, r1 = t1.values[t1.indices], r1.values[r1.indices]
            assert not np.any(np.isnan(t1))
            # assert np.any(np.isnan(r1))
            assert not np.any(np.isnan(t2))
            # assert np.any(np.isnan(r2))
            assert not np.any(np.isnan(t3))
            # assert np.any(np.isnan(r3))

            # net force should be zero
            np.testing.assert_almost_equal(np.sum(t1, axis=0), [0,0,0], decimal=14)

if __name__ == "__main__":
    unittest.main()

ext: py

max_stars_repo_name: silky/timemachine

max_issues_repo_name: silky/timemachine

max_forks_repo_name: silky/timemachine

reason: import numpy


text:
from datashape.user import *
from datashape import dshape
from datetime import date, time, datetime
import numpy as np

def test_validate():
    assert validate(int, 1)
    assert validate('int', 1)
    assert validate(str, 'Alice')
    assert validate(dshape('string'), 'Alice')
    assert validate(dshape('int'), 1)
    assert validate(dshape('int')[0], 1)
    assert validate('real', 2.0)
    assert validate('2 * int', (1, 2))
    assert not validate('3 * int', (1, 2))
    assert not validate('2 * int', 2)

def test_nested_iteratables():
    assert validate('2 * 3 * int', [(1, 2, 3), (4, 5, 6)])

def test_numeric_tower():
    assert validate(np.integer, np.int32(1))
    assert validate(np.number, np.int32(1))

def test_validate_dicts():
    assert validate('{x: int, y: int}', {'x': 1, 'y': 2})
    assert not validate('{x: int, y: int}', {'x': 1, 'y': 2.0})
    assert not validate('{x: int, y: int}', {'x': 1, 'z': 2})

    assert validate('var * {x: int, y: int}', [{'x': 1, 'y': 2}])

    assert validate('var * {x: int, y: int}', [{'x': 1, 'y': 2},
                                               {'x': 3, 'y': 4}])

def test_tuples_can_be_records_too():
    assert validate('{x: int, y: real}', (1, 2.0))
    assert not validate('{x: int, y: real}', (1.0, 2))

def test_datetimes():
    assert validate('time', time(12, 0, 0))
    assert validate('date', date(1999, 1, 20))
    assert validate('datetime', datetime(1999, 1, 20, 12, 0, 0))

def test_numpy():
    assert validate('2 * int32', np.array([1, 2], dtype='int32'))

def test_issubschema():
    assert issubschema('int', 'int')
    assert not issubschema('int', 'float32')

    assert issubschema('2 * int', '2 * int')
    assert not issubschema('2 * int', '3 * int')

    # assert issubschema('float32', 'real')

def test_integration():
    assert validate('{name: string, arrived: date}',
                    {'name': 'Alice', 'arrived': date(2012, 1, 5)})

max_stars_repo_name: mrocklin/datashape

max_issues_repo_name: mrocklin/datashape

max_forks_repo_name: mrocklin/datashape


text:
"""
Copyright (c) 2018-2020 Intel Corporation

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

      http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
"""

import numpy as np
from .format_converter import FileBasedAnnotationConverter, ConverterReturn
from ..representation import PoseEstimationAnnotation
from ..utils import read_xml, check_file_existence
from ..config import PathField

LABELS_TO_COCO = {
    'nose': 0,
    'r_shoulder': 6,
    'r_elbow': 8,
    'r_wrist': 9,
    'l_shoulder': 5,
    'l_elbow': 7,
    'l_wrist': 10,
    'r_hip': 12,
    'r_knee': 14,
    'r_ankle': 16,
    'l_hip': 11,
    'l_knee': 13,
    'l_ankle': 15,
    'r_eye': 2,
    'l_eye': 1,
    'r_ear': 3,
    'l_ear': 4
}

class CVATPoseEstimationConverter(FileBasedAnnotationConverter):
    __provider__ = 'cvat_pose_estimation'
    annotation_types = (PoseEstimationAnnotation, )

    @classmethod
    def parameters(cls):
        configuration_parameters = super().parameters()
        configuration_parameters.update({
            'images_dir': PathField(
                is_directory=True, optional=True,
                description='path to dataset images, used only for content existence check'
            )
        })
        return configuration_parameters

    def configure(self):
        super().configure()
        self.images_dir = self.get_value_from_config('images_dir') or self.annotation_file.parent

    def convert(self, check_content=False, progress_callback=None, progress_interval=100, **kwargs):
        annotation = read_xml(self.annotation_file)
        meta = annotation.find('meta')
        size = int(meta.find('task').find('size').text)
        annotations = []
        content_errors = None if not check_content else []
        for image_id, image in enumerate(annotation.iter('image')):
            identifier = image.attrib['name'].split('/')[-1]
            if check_content:
                if not check_file_existence(self.images_dir / identifier):
                    content_errors.append('{}: does not exist'.format(self.images_dir / identifier))
            label = [1]
            x_vals, y_vals = np.zeros((1, len(LABELS_TO_COCO))), np.zeros((1, len(LABELS_TO_COCO)))
            visilibity = np.zeros((1, len(LABELS_TO_COCO)))
            for point in image.iter('points'):
                point_label = point.attrib['label']
                if point_label not in LABELS_TO_COCO:
                    continue
                point_id = LABELS_TO_COCO[point_label]
                point_x, point_y = point.attrib['points'].split(',')
                x_vals[0, point_id] = float(point_x)
                y_vals[0, point_id] = float(point_y)
                if int(point.attrib['occluded']):
                    continue
                visilibity[0, point_id] = 2
            annotations.append(PoseEstimationAnnotation(identifier, x_vals, y_vals, visilibity, label))

            if progress_callback is not None and image_id % progress_interval == 0:
                progress_callback(image_id * 100 / size)
        meta = {'label_map': {1: 'person'}}

        return ConverterReturn(annotations, meta, content_errors)

    @staticmethod
    def get_pose(image_annotation, num_landmarks):
        landmarks_x, landmarks_y = np.zeros(num_landmarks), np.zeros(num_landmarks)
        for point in image_annotation:
            idx = int(point.attrib['label'])
            x, y = point.attrib['points'].split(',')
            landmarks_x[idx] = float(x)
            landmarks_y[idx] = float(y)

        return landmarks_x, landmarks_y


text:
import json
import os
import cv2
import json
import skimage
import numpy as np
import pandas as pd

from mrcnn.config import Config
from mrcnn import utils

class MolesConfig(Config):
    """Derives from the base Config class and overrides values specific
    to the Moles dataset.
    """
    # Give the configuration a recognizable name
    NAME = "moles"

    # Train on 1 GPU and 8 images per GPU. We can put multiple images on each
    # GPU because the images are small. Batch size is 8 (GPUs * images/GPU).
    GPU_COUNT = 1
    IMAGES_PER_GPU = 4

    # Number of classes (including background)
    NUM_CLASSES = 1 + 3  # background + nevus + melanoma + sk

    # Uses small images for faster training. Set the limits of the small side
    # the large side, and that determines the image shape.
    IMAGE_MIN_DIM = 512
    IMAGE_MAX_DIM = 512

    # Use smaller anchors because our image and objects are small
    RPN_ANCHOR_SCALES = (32, 64, 128, 256, 512)  # anchor side in pixels

    # Reduce training ROIs per image because the images are small and have
    # few objects. Aim to allow ROI sampling to pick 33% positive ROIs.
    TRAIN_ROIS_PER_IMAGE = 64

    # Minimum probability value to accept a detected instance
    # ROIs below this threshold are skipped
    DETECTION_MIN_CONFIDENCE = 0.8

    # Number of training steps per epoch
    STEPS_PER_EPOCH = 400

    # Validation steps at end of epoch
    VALIDATION_STEPS = 5

    LOSS_WEIGHTS = {
        "rpn_class_loss": 1.,
        "rpn_bbox_loss": 1.,
        "mrcnn_class_loss": 2.,
        "mrcnn_bbox_loss": 1.,
        "mrcnn_mask_loss": 1.
    }

    MEAN_PIXEL = np.array([180.4, 150.8, 139.4])

config = MolesConfig()

class MolesDataset(utils.Dataset):
    def load_moles(self, data_path):
        # Add classes
        self.add_class("moles", 1, "malignant")
        self.add_class("moles", 2, "benign")

        if not os.path.exists(data_path):
            raise Exception(data_path + " Does not exists")

        for filename in os.listdir(os.path.join(data_path, "Descriptions")):
            #TODO Use real image ids
            data = json.load(open(os.path.join(data_path,"Descriptions",filename)))

            image_id = int(filename[-7:])

            image_filename = os.path.join(data_path, "Images", filename + ".jpeg")
            if not os.path.isfile(image_filename):
                continue

            mask_filename = os.path.join(data_path, "Segmentation", filename + "_expert.png")
            if not os.path.isfile(mask_filename):
                continue

            metadata = data["meta"]
            mole_type = metadata["clinical"]["benign_malignant"]
            if (mole_type != "benign") and (mole_type != "malignant"):
                continue

            self.add_image("moles", image_id, path=image_filename, mask_path=mask_filename,mole_type=mole_type)

    def load_image(self, image_id):
        """Load the specified image and return a [H,W,3] Numpy array.
        """
        # Load image

        image = skimage.io.imread(self.image_info[image_id]['path'])
        image = cv2.resize(image, (256,256))
        return image

    def load_mask(self, image_id):
        info  = self.image_info[image_id]
        mask_path = info["mask_path"]

        mask = cv2.imread(mask_path)
        mask = cv2.resize(mask, (config.IMAGE_MIN_DIM, config.IMAGE_MIN_DIM))
        mask = mask[:,:,0] > 11
        mask = mask[...,np.newaxis]

        class_id = self.class_names.index(info["mole_type"])
        class_ids = np.array([class_id])
        return mask.astype(np.bool), class_ids.astype(np.int32)

    def image_reference(self, image_id):
        if self.image_info[image_id]["source"] == "moles":
            return self.image_info[image_id]["mask_path"]
        else:
            super(self.__class__).image_reference(self, image_id)

class MolesDatasetFast(MolesDataset):
    def __init__(self):
        self.images = {}
        self.masks = {}
        super().__init__()

    def load_moles(self, data_path):
        # Add classes
        self.add_class("moles", 1, "malignant")
        self.add_class("moles", 2, "benign")

        if not os.path.exists(data_path):
            raise Exception(data_path + " Does not exists")

        for filename in os.listdir(os.path.join(data_path, "Descriptions")):
            data = json.load(open(os.path.join(data_path,"Descriptions",filename)))

            image_id = int(filename[-7:])

            image_filename = os.path.join(data_path, "Images", filename + ".jpeg")
            if not os.path.isfile(image_filename):
                continue

            mask_filename = os.path.join(data_path, "Segmentation", filename + "_expert.png")
            if not os.path.isfile(mask_filename):
                continue

            metadata = data["meta"]
            mole_type = metadata["clinical"]["benign_malignant"]
            if (mole_type != "benign") and (mole_type != "malignant"):
                continue

            self.add_image("moles", image_id, path=image_filename, mask_path=mask_filename,mole_type=mole_type)

    def load_image(self, image_id):
        """Load the specified image and return a [H,W,3] Numpy array.
        """
        if image_id in self.images:
            return self.images[image_id]

        # Load image
        image = skimage.io.imread(self.image_info[image_id]['path'])
        image = cv2.resize(image, (config.IMAGE_MIN_DIM, config.IMAGE_MIN_DIM))

        # Add to dict
        self.images[image_id] = image

        return image

    def load_mask(self, image_id):

        # Check if it's in the dictionary
        if image_id in self.masks:
            return self.masks[image_id]

        # If not then load mask
        mask = super().load_mask(image_id)

        # Add to dict
        self.masks[image_id] = mask

        return mask

class BalancedDataset(MolesDataset):
    def load_moles(self, data_path):
        # Add classes
        self.add_class("moles", 1, "malignant")
        self.add_class("moles", 2, "benign")

        if not os.path.exists(data_path):
            raise Exception(data_path + " Does not exists")

        cap = 3000
        b_count = 0

        for filename in os.listdir(os.path.join(data_path, "Descriptions")):
            #TODO Use real image ids
            data = json.load(open(os.path.join(data_path,"Descriptions",filename)))

            image_id = int(filename[-7:])

            image_filename = os.path.join(data_path, "Images", filename + ".jpeg")
            if not os.path.isfile(image_filename):
                continue

            mask_filename = os.path.join(data_path, "Segmentation", filename + "_expert.png")
            if not os.path.isfile(mask_filename):
                continue

            metadata = data["meta"]
            mole_type = metadata["clinical"]["benign_malignant"]
            if (mole_type == "benign") and b_count < cap:
                b_count += 1
            elif (mole_type == "malignant"):
                pass
            else:
                continue

            self.add_image("moles", image_id, path=image_filename, mask_path=mask_filename,mole_type=mole_type)

class BalancedDatasetFast(BalancedDataset):
    def __init__(self):
        self.images = {}
        self.masks = {}
        super().__init__()

    def load_image(self, image_id):
        """Load the specified image and return a [H,W,3] Numpy array.
        """
        if image_id in self.images:
            return self.images[image_id]

        # Load image
        image = skimage.io.imread(self.image_info[image_id]['path'])
        image = cv2.resize(image, (256,256))

        # Add to dict
        self.images[image_id] = image

        return image

    def load_mask(self, image_id):
        # Check if it's in the dictionary
        if image_id in self.masks:
            return self.masks[image_id]

        # If not then load mask
        mask = super().load_mask(image_id)

        # Add to dict
        self.masks[image_id] = mask

        return mask

class ISIC17Dataset(MolesDatasetFast):
     def load_moles(self, data_path):
        # Add classes
        self.add_class("moles", 1, "nevus")
        self.add_class("moles", 2, "melanoma")
        self.add_class("moles", 3, "seborrheic_keratosis")

        classes = {1:"nevus", 2:"melanoma", 3:"seborrheic_keratosis"}
        if not os.path.exists(data_path):
            raise Exception(data_path + " Does not exists")

        ground_truth_file = os.path.join(data_path, "../", "ISIC-2017_Training_Part3_GroundTruth.csv")
        ground_truth_data = pd.read_csv(ground_truth_file)

        images_path = os.path.join(data_path, "images")
        masks_path = os.path.join(data_path, "masks")

        for filename in os.listdir(images_path):
            #TODO Use real image ids
            #data = json.load(open(os.path.join(data_path,"Descriptions",filename)))

            image_id = int(filename[-11:-4])
            image_name = filename[:12]

            image_filename = os.path.join(images_path, image_name + ".jpg")
            if not os.path.isfile(image_filename):
                continue

            mask_filename = os.path.join(masks_path, image_name + "_segmentation.png")
            if not os.path.isfile(mask_filename):
                continue

            image_info = ground_truth_data[ground_truth_data["image_id"] == image_name]
            label = int((image_info["melanoma"] * 1 + image_info["seborrheic_keratosis"] * 2) + 1)

            self.add_image("moles", image_id, path=image_filename, mask_path=mask_filename,mole_type=classes[label])

class ISIC17AugDataset(MolesDatasetFast):
    def load_moles(self, data_path):
        # Add classes
        self.add_class("moles", 1, "nevus")
        self.add_class("moles", 2, "melanoma")
        self.add_class("moles", 3, "seborrheic_keratosis")

        classes = {1:"nevus", 2:"melanoma", 3:"seborrheic_keratosis"}
        if not os.path.exists(data_path):
            raise Exception(data_path + " Does not exists")

        ground_truth_file = os.path.join(data_path, "labels.csv")
        ground_truth_data = pd.read_csv(ground_truth_file)

        images_path = os.path.join(data_path, "images")
        masks_path = os.path.join(data_path, "masks")

        for filename in os.listdir(images_path):
            #TODO Use real image ids

            image_id = int(filename[-11:-4])
            image_name = filename[:12]

            image_filename = os.path.join(images_path, image_name + ".jpg")
            if not os.path.isfile(image_filename):
                continue

            mask_filename = os.path.join(masks_path, image_name + "_segmentation.png")
            if not os.path.isfile(mask_filename):
                continue

            image_info = ground_truth_data[ground_truth_data["image_id"] == image_name]
            label = int(image_info["label"])

            self.add_image("moles", image_id, path=image_filename, mask_path=mask_filename,mole_type=classes[label])


text:
import numpy as np
from error_calculation import variance_of_matrix

import logging
import sys
import traceback

class inversion:
    def __init__(self,
                 A, f):
        self.A = A
        self.f = f

    def __call__(self):
        Ap_inv = np.linalg.pinv(self.A)
        mus = self.f * Ap_inv
        mus = np.asarray(mus)[0]
        try:
            var = variance_of_matrix(self.A.T)
        except np.linalg.LinAlgError:
            exc_type, exc_value, exc_traceback = sys.exc_info()
            lines = traceback.format_exception(exc_type,
                                               exc_value,
                                               exc_traceback)
            logging.error(''.join(line for line in lines))
            logging.info('Solving for sparse matrix instead')
            from scipy.sparse.linalg import lsqr
            x = lsqr(self.A.T, self.f, calc_var=True)
            mus = x[0]
            var = x[-1]
        self.reco = mus
        self.var = var
        return mus, var

    def reco_hist(self):
        if self.unfolded:
            return self.reco
        else:
            self.__call__()
            return self.reco

reason: import numpy,from scipy


text:
import os
import argparse
import scipy.misc
import scipy.ndimage

import numpy as np

from matplotlib import colors as mcolors

voc_cls = {'background':0,
           'aeroplane': 2,# 1,
           'bus':3, #6,
           'car':8, #7,
           'cat':7, # 8,
           'dog':6, #12,
           'horse':5, #13,
           'person':1, # 15,
           'tv/monitor':4} # 20}
cls_color = {
    0:  [0, 0, 0],
    2:  [128, 0, 0],
    3:  [0, 128, 128],
    8:  [128, 128, 128],
    7:  [64, 0, 0],
    6: [64, 0, 128],
    5: [192, 0, 128],
    1: [192, 128, 128],
    4: [0, 64, 128]
}
def mask_edge_detection(mask, edge_width):
    h = mask.shape[0]
    w = mask.shape[1]

    edge_mask = np.zeros((h,w))

    for i in range(h):
        for j in range(1,w):
            j_prev = j - 1
            # horizantal #
            if not mask[i][j] == mask[i][j_prev]: # horizontal
                if mask[i][j]==1: # 0 -> 1
                    edge_mask[i][j] = 1
                    for add in range(1,edge_width):
                        if j + add < w and mask[i][j+add] == 1:
                            edge_mask[i][j+add] = 1

                else : # 1 -> 0
                    edge_mask[i][j_prev] = 1
                    for minus in range(1,edge_width):
                        if j_prev - minus >= 0 and mask[i][j_prev - minus] == 1:
                            edge_mask[i][j_prev - minus] = 1
            # vertical #
            if not i == 0 :
                i_prev = i - 1
                if not mask[i][j] == mask[i_prev][j]:
                    if mask[i][j]==1: # 0 -> 1
                        edge_mask[i][j] = 1
                        for add in range(1,edge_width):
                            if i + add < h and mask[i+add][j] == 1:
                                edge_mask[i+add][j] = 1
                    else : # 1 -> 0
                        edge_mask[i_prev][j] = 1
                        for minus in range(1,edge_width):
                            if i_prev - minus >= 0 and mask[i_prev-minus][j] == 1:
                                edge_mask[i_prev-minus][j] == 1
    return edge_mask

def viz_data(im, seg, color, inner_alpha = 0.3, edge_alpha = 1, edge_width = 5):

    edge = mask_edge_detection(seg, edge_width)

    color_mask = np.zeros((edge.shape[0]*edge.shape[1], 3))
    l_loc = np.where(seg.flatten() == 1)[0]
    color_mask[l_loc, : ] = color
    color_mask = np.reshape(color_mask, im.shape)
    mask = np.concatenate((seg[:,:,np.newaxis],seg[:,:,np.newaxis],seg[:,:,np.newaxis]), axis = -1)

    color_edge = np.zeros((edge.shape[0]*edge.shape[1], 3))
    l_col = np.where(edge.flatten() == 1)[0]
    color_edge[l_col,:] = color
    color_edge = np.reshape(color_edge, im.shape)
    edge = np.concatenate((edge[:,:,np.newaxis],edge[:,:,np.newaxis],edge[:,:,np.newaxis]), axis = -1)

    im_new = im*(1-mask) + im*mask*(1-inner_alpha) + color_mask * inner_alpha
    im_new =  im_new*(1-edge) + im_new*edge*(1-edge_alpha) + color_edge*edge_alpha

    return im_new

def arg_parse():
    parser = argparse.ArgumentParser(description='Tools to visualize semantic segmentation map.')

    # Datasets parameters
    parser.add_argument('--img_path', type=str, default='',
                    help="path to RGB image")
    parser.add_argument('--seg_path', type=str, default='',
                    help="path to seg")

    args = parser.parse_args()

    return args

if __name__ == '__main__':
    cmap = cls_color
    args = arg_parse()

    img_path = args.img_path
    seg_path = args.seg_path

    img = scipy.misc.imread(img_path)
    seg = scipy.misc.imread(seg_path)

    cs = np.unique(seg)

    for c in cs:
        if not c == 0:
            mask = np.zeros((img.shape[0], img.shape[1]))
            ind = np.where(seg==c)
            mask[ind[0], ind[1]] = 1
            img = viz_data(img, mask, color=cmap[c])
    scipy.misc.imsave('./exp.png', img)

reason: import numpy,import scipy


text:
# command line arguments
# [preset names, ...] checkpoint file

import numpy as np

def test_model(env, model, record=False):
    from presets import preset
    enable_vae = preset.vae.enable

    env.set_mode(1)

    s_norm = model["s_norm"]
    actor = model["policy"]

    import torch
    T = lambda x: torch.FloatTensor(x)

    if enable_vae:
        from utils.vae import VAE
        vae = VAE(preset.vae.latent_dim, use_gpu=False)

        vae_path = 'results/vae/models/%s.pth' % preset.vae.model
        vae.load_state_dict(torch.load(vae_path))

        with open(vae_path + '.norm.npy', 'rb') as f:
            vae_mean = np.load(f)
            vae_std = np.load(f)

    while True:
        ob = env.reset()
        while True:
            with torch.no_grad():
                obt = T(ob)
                obt_normed = s_norm(obt)
                ac = actor.act_deterministic(obt_normed)
                if enable_vae:
                    decoded = vae.decode(T(ac[0:preset.vae.latent_dim]*vae_std + vae_mean)).cpu().numpy()
                    ac = np.concatenate((decoded, ac[preset.vae.latent_dim:]))

            ob, rwd, done, info = env.step(ac)

            if done:
                break

if __name__=="__main__":
    import sys

    from presets import preset
    preset.load_default()

    exp_settings = preset.experiment
    env_settings = preset.env
    env_settings.enable_rendering = True

    for i in range(1, len(sys.argv) - 1):
        preset.load_custom(sys.argv[i])

    preset.load_env_override()

    checkpoint = sys.argv[-1]

    from env import get_env
    test_env = get_env(exp_settings.env)(seed=0, evaluate=True)

    from algorithm import algorithm_bundle_dict
    algorithm_name = preset.experiment.algorithm
    model_class = algorithm_bundle_dict[algorithm_name].model_class

    from model import model_dict
    load_model = model_dict[model_class].loader

    model = load_model(checkpoint)
    test_model(test_env, model, exp_settings.record_motion)

max_stars_repo_name: arpspoof/Jump

max_issues_repo_name: squalidux/Jump

max_forks_repo_name: squalidux/Jump


text:
#!/usr/bin/env python
import random
import re, os, sys, joblib
import pandas as pd
import numpy as np
from collections import defaultdict
from functools import reduce
import codefast as cf

max_stars_repo_name: GaoangLiu/pyendlessh

max_issues_repo_name: GaoangLiu/pyendlessh

max_forks_repo_name: GaoangLiu/pyendlessh


text:
import sys
import argparse
import os
import subprocess
import numpy as np
import time
from iwal import utils
import re

def run_command(cmd):
    p = subprocess.Popen(cmd, shell=True, stderr=subprocess.PIPE)
    a = p.stderr.readlines()
    return a

def sed_maker(start, end, file):
    assert start < end
    return "sed -n '{},{}p;{}q' {}".format(start, end, end, file)

parser = argparse.ArgumentParser()
parser.add_argument('--data', help='dataset', type=str, choices=['cifar64', 'mnist_half', 'fashion', 'svhn'])
parser.add_argument('--batch', help='how often to save the model', type=int)
parser.add_argument('--passes', help='passes', type=int)
parser.add_argument('--c0', nargs='+', type=float, default=[1e-5, 1e-4, 1e-3, 1e-2, 1e-1, 1])
parser.add_argument('--ora', default=False, action='store_true')
opts = parser.parse_args()
# TODO alpha and beta hyperparameters, see parsimonious appendix
# TODO learning rates

for c0 in opts.c0:
    run_command('rm -r tmp/*')

    print("=" * 20 + " C0 = {} ".format(c0) + "=" * 20)
    run_command("rm -r tmp/*")

    if opts.data == 'cifar64':
        X_tr, Y_tr, X_te, Y_te = utils.get_cifar64_data()
        tr_idxs = np.random.permutation(X_tr.shape[0])
        utils.vw_format(X_tr[tr_idxs], Y_tr[tr_idxs], 'cifar_train')
        utils.vw_format(X_te, Y_te, 'cifar_test')
        train_set = 'cifar_train.vw'
        test_set = 'cifar_test.vw'
        train_set_passes = 'cifar_train_passes.vw'
    elif opts.data == 'mnist_half':
        X_tr, Y_tr, X_te, Y_te = utils.get_MNISThalf()
        tr_idxs = np.random.permutation(X_tr.shape[0])
        utils.vw_format(X_tr[tr_idxs], Y_tr[tr_idxs], 'mnist_train')
        utils.vw_format(X_te, Y_te, 'mnist_test')
        train_set = 'mnist_train.vw'
        test_set = 'mnist_test.vw'
        train_set_passes = 'mnist_train_passes.vw'
    elif opts.data == 'svhn':
        X_tr, Y_tr, X_te, Y_te = utils.get_svhn_data()
        tr_idxs = np.random.permutation(X_tr.shape[0])
        utils.vw_format(X_tr[tr_idxs], Y_tr[tr_idxs], 'svhn_train')
        utils.vw_format(X_te, Y_te, 'svhn_test')
        train_set = 'svhn_train.vw'
        test_set = 'svhn_test.vw'
        train_set_passes = 'svhn_train_passes.vw'
    elif opts.data == 'fashion':
        X_tr, Y_tr, X_te, Y_te = utils.get_fashion_data()
        tr_idxs = np.random.permutation(X_tr.shape[0])
        utils.vw_format(X_tr[tr_idxs], Y_tr[tr_idxs], 'fashion_train')
        utils.vw_format(X_te, Y_te, 'fashion_test')
        train_set = 'fashion_train.vw'
        test_set = 'fashion_test.vw'
        train_set_passes = 'fashion_train_passes.vw'

    test_data = open(test_set, 'r').readlines()
    test_labels = [float(l.split(' ')[0]) for l in test_data]
    test_len = len(test_labels)

    train_data = open(train_set, 'r').readlines()
    train_labels = np.array([float(l.split(' ')[0]) for l in train_data])
    train_len = len(train_labels)

    batch = opts.batch
    passes = opts.passes

    # duplicate the data
    with open(train_set_passes, 'w') as f:
        for i in range(passes):
            for l in train_data:
                f.write(l)

    vw_initial_train = 'vw -f tmp/{} -P 1 --active_cover --mellowness={} --loss_function=logistic --binary --learning_rate {}' + (
        ' --oracular' if opts.ora else '')
    vw_train = 'vw -f tmp/{} -P 1 --active_cover --mellowness={} --loss_function=logistic --binary --learning_rate {}' + (
        ' --oracular' if opts.ora else '')

    vw_train_predict = 'vw -i tmp/{} -t --loss_function=logistic -p tmp/pred_train_{}.txt --binary {}'
    vw_test_predict = 'vw -i tmp/{} -t --loss_function=logistic -p tmp/pred_test_{}.txt --binary {}'

    checkpoints = []
    num_trained = []
    tr_acc = []
    te_acc = []

    total_rounds = train_labels.shape[0] * passes / batch + 1
    total_queries_taken = 0

    k = 0
    print('total_rounds', total_rounds)
    while k < total_rounds:
        start = 1  # batch*k+1
        end = min(batch * (k + 1), len(train_labels) * passes)
        num_trained.append(end)
        # TRAINING PHASE
        print('running phase', k, start, end)
        cmd = sed_maker(start, end, train_set_passes) + '|' + vw_initial_train.format('model{}.vw'.format(k), c0,
                                                                                .5 if opts.data == 'cifar64' else 1)
        print(cmd)
        a = run_command(cmd)
        # try:
        #     counter = None
        #     for line in a:
        #         if len(line) >= 4 and line.decode()[:4] == 'loss':
        #             counter = []
        #             continue
        #         if counter is not None:
        #             if line.decode() == '\n':
        #                 break
        #             else:
        #                 lst = re.sub(' +', ' ', line.decode()).split(' ')
        #                 counter.append(int(lst[2]) - 1)
        #     queries_taken = len(np.unique(np.array(counter) % train_labels.shape[0]))
        # except:
        #     print('error')
        #     print(a[-1].decode())
        #     print(a)
        #     sys.exit()
        # checkpoints.append(queries_taken)
        # print('took', queries_taken)

        # TEST PHASE
        pred_test_file = 'tmp/pred_test_{}.txt'.format(k)
        pred_train_file = 'tmp/pred_train_{}.txt'.format(k)

        cmd_test = vw_test_predict.format('model{}.vw'.format(k), k, test_set)
        cmd_train = vw_train_predict.format('model{}.vw'.format(k), k, train_set)

        print(cmd_test)
        b_test = run_command(cmd_test)
        while not os.path.exists(pred_test_file):
            print('sleeping')
            time.sleep(1)
        predictions = np.array([int(x) for x in open(pred_test_file, 'r').readlines()])
        te_acc.append(np.sum(test_labels == predictions) / len(predictions))
        print('cmd line test error', b_test[-4], 'accuracy', te_acc[-1])

        b_train = run_command(cmd_train)
        while not os.path.exists(pred_train_file):
            print('sleeping')
            time.sleep(1)
        predictions = np.array([int(x) for x in open(pred_train_file, 'r').readlines()])
        tr_acc.append(np.sum(train_labels == predictions) / len(predictions))
        print('cmd line train error', b_train[-4], 'accuracy', tr_acc[-1])
        print('\n')
        k += 1

    # print(num_trained)
    # print(checkpoints)
    print(te_acc)
    print(tr_acc)

    run_command(cmd)
    time.sleep(1)

    checkpoints = np.zeros(len(tr_acc))
    with open("./queries.txt", "r") as f:
        idxs = []
        for line in f:
            idxs.append(int(line))
        print(len(idxs))
        j = 0
        queried = []
        for i in idxs:
            idx = i % len(train_labels)
            if idx not in queried:
                queried.append(idx)
                checkpoints[i // batch] += 1
    checkpoints = list(np.cumsum(checkpoints).astype(int))
    num_trained = checkpoints
    print(checkpoints)

    with open('vw_{}_{}_{}_{}{}.npy'.format(opts.data, opts.batch, opts.passes, c0, '_ora' if opts.ora else ''),
              'wb') as f:
        np.save(f, [tr_acc, te_acc, num_trained, checkpoints])

max_stars_repo_name: jifanz/ACED

max_issues_repo_name: jifanz/ACED

max_forks_repo_name: jifanz/ACED


text:
# run svd in parallel
import numpy as np
import pandas as pd
import os
import scipy as sp
import math
from scipy import sparse
from numpy import random
import argparse
import time
from joblib import Parallel, delayed
import pyfftw
import ray
import logging

logging.basicConfig(
    format="%(asctime)s - %(levelname)s - %(name)s - %(message)s",
    datefmt="%m/%d/%Y %H:%M:%S",
    level=logging.INFO,
)
logger = logging.getLogger(__name__)

def Parser():
   parser=argparse.ArgumentParser('')
   parser.add_argument('-i', help = 'input file path')
   parser.add_argument('-o', help = 'output file path')
   parser.add_argument('-s', help = 'suffix')
   parser.add_argument('-t', help = 'tolerance of change of gradient',default = 1e-4,type=float)
   parser.add_argument('-max_iter', help = 'max iteration',default = 500,type=int)
   parser.add_argument('-mu', help = 'stepsize for dual variable',default = 1e-4,type=float)
   parser.add_argument('-max_mu', help = 'max stepsize for dual variable',default = 1e10,type=int)
   parser.add_argument('-rho', help = '>=1, increase mu',default = 1.1,type=float)
   parser.add_argument('-n_core',help= 'number of cores for paralization',default=10,type=int)
   return parser.parse_args()

def tensor_fft(tensor):
    t_s = time.time()
    I = len(tensor)
    J = len(tensor[0])
    K = len(tensor[0][0])
    tmp_tensor = pyfftw.empty_aligned((I,J,K), dtype='complex128')
    fft_tensor = pyfftw.empty_aligned((I,J,K), dtype='complex128')
    fft_object = pyfftw.FFTW(tmp_tensor, fft_tensor,axes=(0,),threads=n_cores)
    tmp_tensor[:,:,:] = tensor
    fft_res = fft_object(tmp_tensor)
    t_e = time.time()
    logger.info("fft time:{}".format(t_e-t_s))
    return fft_res

def tensor_ifft(fft_tensor):
    t_s = time.time()
    I = len(fft_tensor)
    J = len(fft_tensor[0])
    K = len(fft_tensor[0][0])
    tmp_tensor = pyfftw.empty_aligned((I,J,K), dtype='complex128')
    ifft_tensor = pyfftw.empty_aligned((I,J,K), dtype='complex128')
    ifft_object = pyfftw.FFTW(tmp_tensor,ifft_tensor,axes=(0,), direction='FFTW_BACKWARD',threads=n_cores)
    tmp_tensor[:,:,:] = fft_tensor
    ifft_res = ifft_object(tmp_tensor)
    t_e = time.time()
    logger.info("ifft time:{}".format(t_e-t_s))
    return ifft_res

def prox_tnn(Y,rho):
    t_s = time.time()
    I,J,K=Y.shape
    X=np.zeros((I,J,K))
    X=X.astype('complex128')
    Y=Y.astype('complex128')
    Y=tensor_fft(Y)
    #Y=np.fft.fft(Y,axis=0)
    tnn=0
    trank=0
    U,S,V = np.linalg.svd(Y[0],hermitian=True)
    r = np.sum(S>rho)
    if r>=1:
        S=S[0:r]-rho
        X[0] = np.dot(np.dot(U[:,0:r],np.diag(S)),V[0:r,])
        tnn += sum(S)
        trank = max(trank,r)
    halfn3 = math.ceil(I/2)
    svd_list = ray_svd(Y[2:(halfn3+1)],n_cores)
    for i in range(2,halfn3+1):
        U,S,V = svd_list[i-2]
        r = np.sum(S>rho)
        if r>=1:
            S=S[0:r]-rho
            X[i] = np.dot(np.dot(U[:,0:r],np.diag(S)),V[0:r,])
            tnn += sum(S)*2
            trank = max(trank,r)
        X[I-i] = X[i].conjugate()

    if I%2 == 0:
        i = halfn3+1
        U,S,V = np.linalg.svd(Y[i],hermitian=True)
        r = np.sum(S>rho)
        if r>=1:
            S=S[0:r]-rho
            X[i] = np.dot(np.dot(U[:,0:r],np.diag(S)),V[0:r,])
            tnn += sum(S)
            trank = max(trank,r)
    tnn /= I
    X=tensor_ifft(X)
    #X=np.fft.ifft(X,axis=0)
    t_e = time.time()
    logger.info("prox_tnn time:{}".format(t_e-t_s))
    return X,tnn

def opt(M,omega,tol,max_iter,mu,max_mu,rho):
    I,J,K=M.shape
    X=np.zeros((I,J,K))
    X=X.astype('complex128')
    for i in omega:
        X[i] = M[i]
    E=np.zeros((I,J,K))
    Y=E
    iter = 0
    log = []
    t_p=time.time()
    for iter in range(max_iter):
        Xk = X
        Ek = E
        X,tnnX=prox_tnn(-E+M+Y/mu,1/mu)
        E=M-X+Y/mu
        for i in omega:
            E[i] = 0
        dY = M-X-E
        D = I*J*K
        chgX = np.max( np.abs(Xk.reshape(D,)-X.reshape(D,)) )
        chgE = np.max( np.abs(Ek.reshape(D,)-E.reshape(D,)) )
        chg = np.max(  [ chgX,chgE,np.max(np.abs(dY.reshape(D,)))  ] )
        if chg < tol:
            break
        Y = Y+mu*dY
        mu = min(rho*mu,max_mu)
        #logger.info(X)
        logger.info("Iter:%0.2f" % iter)
        logger.info("Error: %0.2f" % np.linalg.norm( dY.reshape(D,) ))
        log.append([iter,np.linalg.norm( dY.reshape(D,) ),chg])
        t_e = time.time()
        logger.info(t_e-t_p)
        t_p=t_e
    obj = tnnX
    err = np.linalg.norm( dY.reshape(D,) )
    return X,log

@ray.remote
def mysvd(mat,idx):
    if idx%50 == 0:
        logger.info('Cell {} starts'.format(idx))
    #logger.info('Cell {} starts'.format(idx))
    res = np.linalg.svd(mat,hermitian=True)
    return res

def par_svd(tensor,n_core):
    res = Parallel(n_jobs=n_core)( delayed(mysvd)(i,idx) for idx,i in enumerate(tensor))
    return res

def ray_svd(tensor,n_core):
    ray.put(tensor)
    N=tensor.shape[0]
    res = ray.get([mysvd.remote(i,idx) for idx,i in enumerate(tensor) ])
    return res

args = Parser()
input_path = args.i
output_path = args.o
tol = args.t
max_iter = args.max_iter
mu = args.mu
max_mu = args.max_mu
rho = args.rho
suffix = args.s
n_cores = args.n_core

#initialize ray
ray.init(num_cpus=n_cores)
# read in data and store as a tensor
sc_tensor = []
all_file = os.listdir(input_path)
#KR_file = [i for i in all_file if 'RawCount' in i]
KR_file = all_file
KR_file.sort()
#KR_file=[]
#for i in range(1,41):
#   KR_file.append('slice_'+str(i)+'.txt')

# logger.info raw data file name
logger.info('These files are in use:\n')
[logger.info(i) for i in KR_file]
KR_file_frame = pd.DataFrame(KR_file)
KR_file_frame.to_csv(output_path+'/'+suffix+'_input_file_index.csv')
sc_tensor = pd.read_csv(input_path+"/"+KR_file[0],header=None,sep='\t')
sc_tensor = sc_tensor.values
N = sc_tensor.shape[0]
sc_tensor = np.empty(shape=[len(KR_file),N,N])

for idx,f in enumerate(KR_file):
    tmp_file = pd.read_csv(input_path+"/"+f,header=None,sep='\t')
    tmp_file = tmp_file.values
    # check if the input data has same dimensions
    if idx == 0:
        N = len(tmp_file)
    else:
        if len(tmp_file) != N:
            raise ValueError('scHi-C in this chromosome has different dimensions')
    sc_tensor[idx,:,:]=tmp_file
    logger.info('{idx} outof {n_file} done'.format(idx=idx,n_file=len(KR_file)))

sc_tensor=np.array(sc_tensor) # A[i,j,k] represent (j,k) element of single cell i
# define measurement set: non zero entries in sc_tensor
I = len(sc_tensor)
J = len(sc_tensor[0])
K = J
omega = set()
non_zero_entries=np.transpose(sc_tensor.nonzero())
for i,j,k in non_zero_entries:
    omega.add( (i,k,j) )

logger.info(len(omega))
X,log = opt(sc_tensor,omega,tol,max_iter,mu,max_mu,rho)
np.save(output_path+'/'+suffix, X)
log = pd.DataFrame(log)
log.to_csv(output_path+'/'+suffix+".txt",sep='\t')

reason: import numpy,from numpy,import scipy,from scipy


text:
import pandas as pd
import numpy as np
from datasets import Dataset, load_metric
from sklearn.model_selection import train_test_split
from transformers import TrainingArguments, Trainer
from transformers import BertForSequenceClassification, BertTokenizer
import json
from transformers import TrainerCallback, EarlyStoppingCallback

#get test split from aman
test = pd.read_csv('test')
custom_dataset_train = Dataset.from_pandas(test)

expert1 = BertForSequenceClassification.from_pretrained('', num_labels=2)
expert2 = BertForSequenceClassification.from_pretrained('', num_labels=2)
expert3 = BertForSequenceClassification.from_pretrained('', num_labels=2)

predictions = expert1(**test)
preds = np.argmax(predictions.predictions, axis=-1)

#inputs : 3 predictions (one from each expert)
#MLP
import torch.optim as optim
import torch.nn as nn
import torch.nn.functional as F

criterion = nn.BinaryCrossEntropyLoss()
optimizer = optim.SGD(net.parameters(), lr=0.001, momentum=0.9)

class BinaryClassification(nn.Module):
    def __init__(self):
        super(BinaryClassification, self).__init__()
        self.layer_1 = nn.Linear(3, 16)
        self.layer_2 = nn.Linear(16, 16)
        self.layer_out = nn.Linear(16, 1)
        self.sig = nn.Sigmoid()
        # self.dropout = nn.Dropout(p=0.1)
        # self.batchnorm1 = nn.BatchNorm1d(64)
        # self.batchnorm2 = nn.BatchNorm1d(64)

    def forward(self, inputs):
        x = F.relu(self.layer_1(inputs))
        x = F.relu(self.layer_2(x))
        # x = self.dropout(x)
        x = self.layer_out(x)
        x = self.sig(x)

        return x

#each entry shld be three cached predictions and target
#Heiarcharical

for epoch in range(2):  # loop over the dataset multiple times

    running_loss = 0.0
    for i, data in enumerate(trainloader, 0):
        # get the inputs; data is a list of [inputs, labels]
        inputs, labels = data

        # zero the parameter gradients
        optimizer.zero_grad()

        # forward + backward + optimize
        outputs = net(inputs)
        loss = criterion(outputs, labels)
        loss.backward()
        optimizer.step()

        # print statistics
        running_loss += loss.item()
        if i % 2000 == 1999:    # print every 2000 mini-batches
            print(f'[{epoch + 1}, {i + 1:5d}] loss: {running_loss / 2000:.3f}')
            running_loss = 0.0

print('Finished Training')


text:
"""Swap updator for spin state with value encoding."""
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

from typing import Tuple

import numpy as np
import tensorflow as tf

from updators import updator

class SpinSwapUpdatorValue(updator.Updator):

    def __init__(self, shape: Tuple[int, ...]) -> None:
        """Initializes swap updator for spin state with value encoding.

        Args:
            shape: The shape of the spin state.
        """
        self._shape = shape
        self._n_sites = int(np.prod(self._shape))

    def masks_placeholder(self) -> tf.Tensor:
        return tf.placeholder(dtype=tf.int32, shape=[None, self._n_sites],
                              name='Masks')

    def generate_masks(self, n_sample: int) -> np.ndarray:
        rang = range(n_sample)
        swaps = np.random.randint(0, self._n_sites, (2, n_sample))
        masks = np.arange(self._n_sites)[None, :].repeat(n_sample, axis=0)
        masks[rang, swaps[0]], masks[rang, swaps[1]] = (
            masks[rang, swaps[1]], masks[rang, swaps[0]])
        return masks

    def __call__(self, state: tf.Tensor, mask: tf.Tensor) -> tf.Tensor:
        state = tf.reshape(state, [-1])
        return tf.reshape(tf.gather(state, mask), self._shape)


text:
"""
    Code depeloved by Laura Cabayol
    email: lcabayol@ifae.es

"""

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.preprocessing import normalize, RobustScaler, scale, StandardScaler,MinMaxScaler
import keras
from keras.utils import to_categorical
from keras.models import Sequential,Input,Model
from keras.layers import Dense, Dropout, Flatten, Activation,Conv1D, MaxPooling1D
from keras.layers.advanced_activations import LeakyReLU
from sklearn.metrics import roc_curve, precision_recall_curve, roc_auc_score, average_precision_score,average_precision_score
from sklearn.utils import shuffle
from imblearn.over_sampling import SMOTE
from keras.callbacks import *

def dataset_dealer(dataset):
    """takes the initial dataset, applies a cut on magnitude and creates two differentiated datasets: one for stars and another for galaxies"""
    print('The magnitude cut is at', magCut)
    dataset = dataset[dataset.i_auto < magCut]
    dataset = random_sample(dataset)#shuffle the dataset to avoid having stars and galaxies clustered
    dataset = dataset.dropna()#avoid any object that can cause problems in the training
    dataset = dataset.reset_index(drop = True)
    Stars = dataset[dataset.Label == 1]
    Stars = Stars.reset_index(drop = False)
    Galaxies = dataset[dataset.Label == 0]
    Galaxies = Galaxies.reset_index(drop = False)
    return Stars, Galaxies

def sample_split(Galaxies, Stars, nStarsTRN, nStarsVal, nGalTRN, nGalVal):
    """#takes the star and galaxy datasets and generates a training and validation sample with the number of objects per class defined before"""
    train = pd.concat((Galaxies.loc[1:(nGalTRN+1), : ],Stars.loc[1:(nStarsTRN+1), :]), axis = 0,ignore_index ='True')
    val = pd.concat((Galaxies.loc[(nGalTRN+1): (nGalTRN + 1 + nGalVal) ,:],Stars.loc[(nStarsTRN+1): (nStarsTRN + 1 + nStarsVal) ,: ]), axis = 0, ignore_index ='True' )
    train = random_sample(train)
    val = random_sample(val)
    return train, val

def random_sample(DataFrame):
    DataFrame = shuffle(DataFrame)
    DataFrame = DataFrame.reset_index(drop = True)
    return DataFrame

def Flux_input(Fluxes):
    """ Format that the network needs"""
    Fluxes = Fluxes.reshape(-1, Fluxes.shape[1] ,1)
    Fluxes = Fluxes.astype('float32')
    return Fluxes

def CNN():
    ConvNNet = Sequential()
    ConvNNet.add(Conv1D(filters = 32, kernel_size = 10,activation='linear',input_shape=(FeaturesTRN.shape[1],1),padding='same'))
    ConvNNet.add(LeakyReLU(alpha=0.1))
    ConvNNet.add(MaxPooling1D(pool_size = 4, padding='same'))
    ConvNNet.add(Conv1D(filters =64, kernel_size = 5, activation='linear',padding='same'))
    ConvNNet.add(LeakyReLU(alpha=0.1))
    ConvNNet.add(MaxPooling1D(pool_size= 2, padding='same'))
    ConvNNet.add(Conv1D(filters = 128, kernel_size =3 , activation='linear',padding='same'))
    ConvNNet.add(LeakyReLU(alpha=0.1))
    ConvNNet.add(MaxPooling1D(pool_size= 2 ,padding='same'))
    ConvNNet.add(Flatten())
    ConvNNet.add(Dense(128, activation='linear'))
    ConvNNet.add(LeakyReLU(alpha=0.1))
    #ConvNNet.add(Dropout(0.1))
    ConvNNet.add(Dense(num_classes, activation='softmax'))
    return ConvNNet

#---------------------------------------------MAIN-------------------------------------------

#-------------------parameters--------------------
nStarsTRN = 5000                                #|
nGalTRN = 15000                                 #|
nStarsVal = 10000                               #|
nGalVal = 10000                                 #|
                                                #|
magCut = 22.5                                   #|
                                                #|
Labels = ['Label']                              #|
magnitude = ['i_auto']                          #|
coord = ['ra','dec']                            #|
magerr = ['magerr']                             #|
                                                #|
batch_size = 8                                  #|
epochs = 5                                      #|
num_classes = 2                                 #|
                                                #|
bands = ['NB455','NB465','NB475','NB485','NB495','NB505','NB515','NB525','NB535','NB545','NB555','NB565','NB575','NB585','NB595','NB605','NB615','NB625','NB635','NB645', 'NB655','NB665','NB675','NB685','NB695','NB705','NB715','NB725','NB735','NB745','NB755','NB765','NB775','NB785','NB795','NB805','NB815','NB825','NB835','NB845']
                                                #|
#-----------------------------------------------#|

print('This code classifies stars and galaxies with a ConvNet. The trainig and validation samples must be labeled:')
print(' --> stars = 1')
print(' --> galaxies = 0')
print('The dataset contains the 40NB fluxes, the truth label and i_auto and it consists of PAUs objects on the COSMOS field ')

"""This file should be a csv file containing:
- red_id,
- PAUS fluxes in the 40 NB: [NB455,...,NB845]
- i_auto
- ra
- dec
- Label: From Leauthaud et al. 2007
"""
DataPAU = pd.read_table('YourFile.csv', sep  = ',', comment = '#', header = 0)

Stars, Galaxies = dataset_dealer(DataPAU)#input:original dataset, output: star and galaxy datasets

print('The total number of stars is',Stars.shape[0])
print('The total number of galaxies is',Galaxies.shape[0])

DataTRN, DataVal = sample_split(Galaxies, Stars, nStarsTRN, nStarsVal, nGalTRN, nGalVal)

#Define input features and input labels
FeaturesTRN = DataTRN.loc[:,bands]
LabelsTRN = DataTRN.loc[:,'Label']
FeaturesVal = DataVal.loc[:,bands]
LabelsVal = DataVal.loc[:,'Label']

#input array into the format needed for the cnn
FeaturesTRN = Flux_input(FeaturesTRN.values)
FeaturesVal = Flux_input(FeaturesVal.values)

#input labels must be hot encoded: 0 --> [1,0]; 1 --> [0,1]
LabelsTRN_cat = to_categorical(LabelsTRN)
LabelsVal_cat = to_categorical(LabelsVal)

#----------------------- cnn ---------------------------
ConvNNet = CNN()

ConvNNet.compile(loss=keras.losses.categorical_crossentropy, optimizer=keras.optimizers.Adam(lr =0.001, decay = 0.001),metrics={'output_a': 'accuracy'})#define the network: loss function, optimization method, metrics

print(ConvNNet.summary())

ConvNNet_train = ConvNNet.fit(FeaturesTRN, LabelsTRN_cat,batch_size=batch_size,epochs=epochs,verbose=1,validation_data=(FeaturesVal, LabelsVal_cat), initial_epoch = 0) #train

Stellar_flag_TRN = ConvNNet.predict(FeaturesTRN, batch_size=None, verbose=0) #validate with the network trained
Stellar_flag_Val = ConvNNet.predict(FeaturesVal, batch_size=None, verbose=0) #validate with the network trained


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

from glob import glob
from typing import Dict, List, Any
from collections import defaultdict
from tensorboard.backend.event_processing.event_accumulator import (
    EventAccumulator)
from .utils.parser_utils import ArgparseFormatter
from .utils.logging_utils import stdout_root_logger
from .arg_parser import logging_arg_parser, tensorboard_arg_parser
import numpy as np
import argparse
import os
import csv

def dict2csv(out: Dict[str, List[Any]], dpath: str) -> None:
    """
    Function to write dictionary object as csv file

    Args:
        out: Dictionary containing values/steps to write
        dpath: Path of the directory containing tensorboard logs
    """
    with open(os.path.join(dpath, "%s.csv" % os.path.basename(dpath)),
              "w") as f:
        writer = csv.DictWriter(f, out.keys())
        writer.writeheader()
        for i in range(len(out["steps"])):
            writer.writerow({key: out[key][i] for key in out.keys()})

def tabulate_events(dpath: str) -> Dict[str, List[Any]]:
    """
    Function to tabulate and aggregate event logs into single dictionary
    with post-processing to ensure data sanity

    Args:
        dpath: Path of the directory containing tensorboard logs

    Returns:
        out: Dictionary containing relevant tensorboard data
    """
    summary_iterators = [
        EventAccumulator(os.path.join(dpath, dname)).Reload()
        for dname in os.listdir(dpath) if ".csv" not in dname
    ]
    # find lowest common denominator in tags
    tags = set(summary_iterators[0].Tags()["scalars"])
    for summary_iterator in summary_iterators[1:]:
        tags.intersection_update(summary_iterator.Tags()["scalars"])
    tags = list(tags)
    # create variable dictionary
    out = defaultdict(list)
    # loop over summary iterators
    for summary_iterator in summary_iterators:
        # find lowest common denominator of step elements
        steps_list = [[event.step for event in summary_iterator.Scalars(tag)]
                      for tag in tags]
        steps = set(steps_list[0])
        for inner_list in steps_list[1:]:
            steps.intersection_update(inner_list)
        # sort these in ascending order
        steps = sorted(list(steps))
        out["steps"].extend(steps)
        for tag in tags:
            hold: List[Any] = []
            for event in summary_iterator.Scalars(tag):
                current_steps = [] if hold == [] else list(zip(*hold))[0]
                if event.step in steps and event.step not in current_steps:
                    hold.append([event.step, event.value])
            # sort hold again in case this was disrupted during event access
            hold = [element[1] for element in sorted(hold, key=lambda x: x[0])]
            # ensure all lengths are the same
            assert len(hold) == len(steps)
            # append to out after sorting/cleaning
            out[tag].extend(hold)
    # sort everything based on steps
    sorting_indices = np.argsort(out["steps"])
    for tag in tags + ["steps"]:
        out[tag] = [out[tag][i] for i in sorting_indices]
    return out

def main(args: argparse.Namespace) -> None:
    """ Main function to tabulate tensorboard data and write to disk as csv """
    # parse for tensorboard logs
    tb_event_directories = glob(args.tb_event_directory)
    # loop over log directories
    for tb_event_directory in tb_event_directories:
        if os.path.exists(
                os.path.join(
                    tb_event_directory, "%s.csv" %
                    os.path.basename(tb_event_directory))) and not args.force:
            continue
        else:
            LOGGER.info("Processing: %s", tb_event_directory)
            out = tabulate_events(tb_event_directory)
            LOGGER.info("Writing results to directory: %s", tb_event_directory)
            dict2csv(out, tb_event_directory)

if __name__ == "__main__":
    parser = argparse.ArgumentParser(formatter_class=ArgparseFormatter,
                                     parents=[
                                         tensorboard_arg_parser(),
                                         logging_arg_parser(),
                                     ])
    LOGGER = stdout_root_logger(
        logging_arg_parser().parse_known_args()[0].logging_level)
    main(parser.parse_args())

max_stars_repo_name: atreyasha/spp-explainability

max_issues_repo_name: atreyasha/spp-explainability

max_forks_repo_name: atreyasha/spp-explainability


text:
from functools import partial
from concurrent import futures
import numpy as np
from skimage.transform import resize
from skimage.measure import block_reduce
from tqdm import tqdm

from .util import blocking, grow_bounding_box, open_file

def ds_interpolate(data, scale_factor, out_shape, order):
    dtype = data.dtype
    out = resize(data, out_shape, order=order, mode='constant',
                 anti_aliasing=order > 0, preserve_range=True)
    return out.astype(dtype)

def ds_block_reduce(data, scale_factor, out_shape, function):
    out = block_reduce(data, tuple(scale_factor), function)
    # crop if necessary
    if out.shape != out_shape:
        bb = tuple(slice(0, osh) for osh in out_shape)
        out = out[bb]
    return out

def sample_shape(shape, factor, add_incomplete_blocks=False):
    if add_incomplete_blocks:
        sampled = tuple(sh // scale_factor + int((sh % scale_factor) != 0)
                        for sh, scale_factor in zip(shape, factor))
    else:
        sampled = tuple(sh // scale_factor for sh, scale_factor in zip(shape, factor))
    sampled = tuple(max(1, sh) for sh in sampled)
    return sampled

def get_downsampler(mode):
    if mode == 'nearest':
        downsample_function = partial(ds_interpolate, order=0)
    elif mode == 'mean':
        downsample_function = partial(ds_block_reduce, function=np.mean)
    elif mode == 'max':
        downsample_function = partial(ds_block_reduce, function=np.max)
    elif mode == 'min':
        downsample_function = partial(ds_block_reduce, function=np.min)
    elif mode == 'interpolate':
        downsample_function = partial(ds_interpolate, order=3)
    else:
        raise ValueError("Downsampling mode %s is not supported" % mode)
    return downsample_function

def downsample_in_memory(input_volume,
                         downscale_factors,
                         downscale_mode,
                         block_shape,
                         n_threads):
    downscaled_volumes = []
    downsample_function = get_downsampler(downscale_mode)

    def sample_chunk(bb, in_vol, out_vol, scale_factor, halo):

        # grow the bounding box if we have a halo
        if halo is None:
            bb_grown = bb
            bb_local = np.s_[:]
        else:
            bb_grown, bb_local = grow_bounding_box(bb, halo, shape)

        bb_up = tuple(slice(b.start * scale_factor, b.stop * scale_factor)
                      for b, scale_factor in zip(bb_grown, factor))
        inp = in_vol[bb_up]

        # don't sample empty blocks
        if inp.sum() == 0:
            return

        out_shape = tuple(b.stop - b.start for b in bb_grown)
        outp = downsample_function(inp, factor, out_shape)
        out_vol[bb] = outp[bb_local]

    in_vol = input_volume
    for factor in downscale_factors:
        shape = in_vol.shape

        halo = None
        # TODO for the interpolate mode we should use a halo,
        # otherwise the result is not fully correct at the block boundaries
        # however, the current implementation with halo is not correct
        # halo = 2 * factor if mode == 'interpolate' else None

        ds_shape = sample_shape(shape, factor)
        ds_vol = np.zeros(ds_shape, dtype=input_volume.dtype)

        sampler = partial(sample_chunk,
                          in_vol=in_vol,
                          out_vol=ds_vol,
                          scale_factor=factor,
                          halo=halo)

        with futures.ThreadPoolExecutor(n_threads) as tp:
            list(tp.map(sampler, blocking(ds_shape, block_shape)))

        downscaled_volumes.append(ds_vol)
        in_vol = ds_vol

    return downscaled_volumes

def downsample(path, in_key, out_key, factor, mode, n_threads=1, overwrite=False):
    """ Downsample input volume.
    """

    downsample_function = get_downsampler(mode)
    halo = None
    # TODO for the interpolate mode we should use a halo,
    # otherwise the result is not fully correct at the block boundaries
    # however, the current implementation with halo is not correct
    # halo = 2 * factor if mode == 'interpolate' else None

    with open_file(path, 'a') as f:
        ds_in = f[in_key]
        shape = ds_in.shape
        chunks = ds_in.chunks

        sampled_shape = sample_shape(shape, factor)
        chunks = tuple(min(sh, ch) for sh, ch in zip(sampled_shape, ds_in.chunks))

        if overwrite and out_key in f:
            del f[out_key]

        ds_out = f.create_dataset(out_key, shape=sampled_shape, chunks=chunks,
                                  compression='gzip', dtype=ds_in.dtype)

        def sample_chunk(bb):

            # grow the bounding box if we have a halo
            if halo is None:
                bb_grown = bb
                bb_local = np.s_[:]
            else:
                bb_grown, bb_local = grow_bounding_box(bb, halo, shape)

            bb_up = tuple(slice(b.start * scale_factor, b.stop * scale_factor)
                          for b, scale_factor in zip(bb_grown, factor))
            inp = ds_in[bb_up]

            # don't sample empty blocks
            if inp.sum() == 0:
                return

            out_shape = tuple(b.stop - b.start for b in bb_grown)
            outp = downsample_function(inp, factor, out_shape)
            ds_out[bb] = outp[bb_local]

        blocks = list(blocking(sampled_shape, chunks))
        if n_threads > 1:
            with futures.ThreadPoolExecutor(n_threads) as tp:
                list(tqdm(tp.map(sample_chunk, blocks), total=len(blocks)))
        else:
            for bb in tqdm(blocks, total=len(blocks)):
                sample_chunk(bb)

max_stars_repo_name: jhennies/pybdv

max_issues_repo_name: jhennies/pybdv

max_forks_repo_name: jhennies/pybdv


text:
from __future__ import division
import pandas as pd
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
import numpy.linalg as npla
import scipy.linalg as spla
import scipy.special

# Code may print some warnings due to pandas in Pyton 2.7
# It can be run in Python3 as well. The warnings are completely benign.

z = 20

J_50 = scipy.special.jv(50, z)
J_49 = scipy.special.jv(49, z)

bessel_val_by_rec = [J_50, J_49]

for i in range(48, -1, -1):

    J_i = ( ( ( ( 2 * ( i + 1 ) ) / z ) * bessel_val_by_rec[50 - i - 1] ) - bessel_val_by_rec[50 - i - 2] )

    bessel_val_by_rec.append(J_i)

bessel_val_by_func = [scipy.special.jv(i, z) for i in range(50, -1, -1)]

relative_err = []
for i in range(2, 51):
    relative_err.append((bessel_val_by_rec[i] - bessel_val_by_func[i])/bessel_val_by_func[i])

# Again, I am just printing magnitude of the Relative Error
# Comment this line to print just the ratio
relative_err = [abs(i) for i in relative_err]

df = pd.DataFrame({"n" : [i for i in range(2, 51)], "From Function" : bessel_val_by_func[2:], "From Recurrence" : bessel_val_by_rec[2:], "Relative Error" : relative_err}, columns=["n", "From Function", "From Recurrence", "Relative Error"])
print(df.to_string(index_names=False))

plt.xlabel("Value of n")
plt.ylabel("Absolute Value of Relative Error")
plt.title("Problem 3(d)")
plt.plot( np.arange(2, 51), relative_err)
plt.savefig("problem_3d.png")
# Uncomment to show plot
# plt.show()

max_stars_repo_name: justachetan/scientific-computing

max_issues_repo_name: justachetan/scientific-computing

max_forks_repo_name: justachetan/scientific-computing


text:
from ._version import get_versions
from .aggregate_purepy import aggregate as aggregate_py

def dummy_no_impl(*args, **kwargs):
    raise NotImplementedError("You may need to install another package (numpy, "
                              "weave, or numba) to access a working implementation.")

aggregate = aggregate_py

try:
    import numpy as np
except ImportError:
    aggregate_np = aggregate_ufunc = dummy_no_impl
    multi_arange = multi_cumsum = label_contiguous_1d = dummy_no_impl
else:
    from .aggregate_numpy import aggregate
    aggregate_np = aggregate
    from .aggregate_numpy_ufunc import aggregate as aggregate_ufunc
    from .utils_numpy import (label_contiguous_1d, multi_arange, relabel_groups_masked,
                              relabel_groups_unique, unpack)

try:
    try:
        import weave
    except ImportError:
        from scipy import weave
except ImportError:
    aggregate_wv = None
else:
    from .aggregate_weave import aggregate as aggregate_wv, step_count, step_indices
    aggregate = aggregate_wv

try:
    import numba
except ImportError:
    aggregate_nb = None
else:
    from .aggregate_numba import aggregate as aggregate_nb, step_count, step_indices
    aggregate = aggregate_nb

def uaggregate(group_idx, a, **kwargs):
    return unpack(group_idx, aggregate(group_idx, a, **kwargs))

__version__ = get_versions()['version']
del get_versions

reason: import numpy,from scipy,import numba


text:
# Author: Mark McDonnell, mark.mcdonnell@unisa.edu.au
import numpy as np
import config

if config.tf:
    from tensorflow.keras import backend as K
    from tensorflow.keras.layers import InputSpec, Conv2D
else:
    from keras import backend as K
    from keras.layers import InputSpec, Conv2D

# from tensorflow.keras import initializers

class BinaryConv2D(Conv2D):
    '''Binarized Convolution2D layer
    References:
    "BinaryNet: Training Deep Neural Networks with Weights and Activations Constrained to +1 or -1" [http://arxiv.org/abs/1602.02830]

    adapated by Mark McDonnell from https://github.com/DingKe/nn_playground/blob/master/binarynet/binary_layers.py
    '''

    def __init__(self, filters, **kwargs):
        super(BinaryConv2D, self).__init__(filters, **kwargs)

    def build(self, input_shape):
        if self.data_format == 'channels_first':
            channel_axis = 1
        else:
            channel_axis = -1

        input_dim = int(input_shape[channel_axis])

        if input_dim is None:
            raise ValueError('The channel dimension of the inputs '
                             'should be defined. Found `None`.')

        self.multiplier = np.sqrt(
            2.0 / np.float(self.kernel_size[0]) / np.float(self.kernel_size[1]) / float(input_dim))

        self.kernel = self.add_weight(shape=self.kernel_size + (input_dim, self.filters),
                                      initializer=self.kernel_initializer,
                                      name='kernel',
                                      regularizer=self.kernel_regularizer,
                                      constraint=self.kernel_constraint)

        # Set input spec.
        self.input_spec = InputSpec(ndim=4, axes={channel_axis: input_dim})
        self.built = True

    def call(self, inputs):

        binary_kernel = self.kernel + K.stop_gradient(K.sign(self.kernel) - self.kernel)
        binary_kernel = binary_kernel + K.stop_gradient(binary_kernel * self.multiplier - binary_kernel)

        outputs = K.conv2d(inputs,
                           binary_kernel,
                           strides=self.strides,
                           padding=self.padding,
                           data_format=self.data_format,
                           dilation_rate=self.dilation_rate)

        return outputs

    def get_config(self):
        config = {'multiplier': self.multiplier}
        base_config = super(BinaryConv2D, self).get_config()
        return dict(list(base_config.items()) + list(config.items()))


text:
import random
import numpy as np
from scipy.integrate import odeint

import file_utils

def lorenz(state, t):
    x = state[0]
    y = state[1]
    z = state[2]

    sigma = 10.0
    rho = 28.0
    beta = 8.0 / 3.0

    # compute state derivatives
    xd = sigma * (y - x)
    yd = (rho - z) * x - y
    zd = x * y - beta * z

    # return the state derivatives
    return [xd, yd, zd]

def get_lorenz(step, ts, tf, x0):
    td = int((tf - ts) / step)
    t = np.linspace(ts, tf, td)
    x_t = np.asarray([odeint(lorenz, x0i, t) for x0i in x0])
    if len(x_t[0]) == td:
        print('New system calculated')
    return x_t

if __name__ == "__main__":
    step = 0.01
    ts = 0
    tf = 100
    num_traj = 2
    x0 = [(random.randint(-25, 25), random.randint(-35, 35), random.randint(5, 55)) for traj in range(num_traj)]
    x_t = get_lorenz(step, ts, tf, x0)
    file_utils.save_to_file(x_t)
    array = file_utils.open_file()
    print(array == x_t)

max_stars_repo_name: nikolisan/Diploma-Thesis

max_issues_repo_name: nikolisan/Diploma-Thesis

max_forks_repo_name: nikolisan/Diploma-Thesis


text:
import argparse
import os
import numpy as np

parser = argparse.ArgumentParser(description=__doc__)
def get_config():
    config, unparsed = parser.parse_known_args()
    return config, parser

# Feature extraction procedure
rawdatapath = '/storage/ge/voxceleb2'
data_path = '/storageNVME/ge/voxceleb2'
data_arg = parser.add_argument_group('Data')
data_arg.add_argument("-data_i", type=str, default=os.path.join(rawdatapath, 'dev', 'aac'), help="input folder path with raw files")
data_arg.add_argument("-data_o", type=str, default=os.path.join(data_path, 'train'), help="output folder with feature data")
data_arg.add_argument("-sr", type=int, default=16000, help="sampling rate")
data_arg.add_argument("-seg_len", type=float, default=3.9, help="minimum duration of segments(s)")
data_arg.add_argument("-eer_length", type=float, default=3.9, help="duration of eer computing segments(s)")

# Training Parameters
train_arg = parser.add_argument_group('Training')
train_path = '/home/ge/SV-wav-pytorch/experiments/dummy/'
train_arg.add_argument("-train_source", type=str, action='append', help="data source for training")
train_arg.add_argument("-wavencoder", type=str, default='cnn', help="wavform encoder")
train_arg.add_argument("-backbone", type=str, default='rawcnn', help="architecture for training")
train_arg.add_argument("-train_path", type=str, default=train_path, help="train path")
train_arg.add_argument('-epoch', type=int, default=30, help="number of epoch for training")
train_arg.add_argument('-batch_size', type=int, default=128, help="batch size")
train_arg.add_argument('-embed_dim', type=int, default=256, help="embedding vector dimension")
train_arg.add_argument('-loss_type', type=str, default='softmaxloss', help="loss type for training")
train_arg.add_argument("-vd", type=bool, help="whether use variational dropout")
train_arg.add_argument('-final_layer', type=str, default='None', help="loss type for training")
train_arg.add_argument('-optim', type=str, default='SGD', help="optimizer")
train_arg.add_argument('-lr', type=float, default=0.01, help="learning rate")
train_arg.add_argument('-gpu_idx', type=str, default='0', help="gpu index")

# Evaluation Parameters
eva_arg = parser.add_argument_group('Evaluation')
testpath = '/storage/ge/voxceleb/test'
# testpath = '/data/ge/voxceleb/test'
eva_arg.add_argument("-layer", type=str, default='fc2', help="layer name")
eva_arg.add_argument("-eval", type=str, default=os.path.join(testpath, 'eval'), help="eval folder path with raw files")
eva_arg.add_argument("-extract_dataset", type=str, default='vox2', help="dataset name for vector extraction")
eva_arg.add_argument("-eval_njob", type=int, default=3, help="number of eval multipreocessing jobs")
eva_arg.add_argument("-list", type=str, default=os.path.join(testpath, 'wav'), help="eval list, it should also contain eval files")
eva_arg.add_argument("-test_i", type=str, default=os.path.join(testpath, 'test'), help="eval folder path with pickle files")
eva_arg.add_argument("-test_o", type=str, default='/storageNVME/ge/voxceleb', help="eval folder path with pickle files")
config = get_config()


text:
import numpy as np
from skimage import io

def get_tensor_as_image(tensor):
    tensor = tensor[0]
    np_img = tensor.cpu().data.numpy()
    np_img = np.reshape(np_img, (3,150,-1))
    np_img = np.moveaxis(np_img,0,2)
    return np_img

def get_tensor_as_image_grayscale(tensor):
    tensor = tensor[0]
    np_img = tensor.cpu().data.numpy()
    np_img = np.reshape(np_img, (150,-1))
    #np_img = np.moveaxis(np_img,0,2)
    return np.dstack((np_img, np_img, np_img))

def gt_to_cont():
    img = io.imread("/home/ola/projects/masteroppgave/msth-datasets/laser-scan/corner-scan-v1/dataset/ground-truth/img000.png")
    height, width = img.shape
    gt_cont = np.argmax(img, axis=1)
    print(gt_cont.shape)
    print(type(np.argmax(img, axis=1)))

if __name__ == '__main__':
    gt_to_cont()