text:
"""
See LICENSE for details on how to use this code
"""
# Libraries
from functools import lru_cache
import argparse
from math import sin, cos, tan, atan, exp, pi, sqrt
from mpmath import csch, coth
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate
import scipy.signal

class MechanicalProperties:
    """
    A class to represent the mechanical properties of the cell membrane and the particle.

    Attributes:
        ----------
        gamma_bar: float
            adimensional linear adhesion energy between membrane and the particle
        sigma_bar: float
            adimensional membrane tension

    Methods:
        -------
        None

    """

    def __init__(self, gamma_bar, sigma_bar):
        """
        Constructs all the necessary attributes for the mechanical properties object.

        Parameters:
            ----------
            gamma_bar: float
                adimensional linear adhesion energy between membrane and the particle
            sigma_bar: float
                adimensional membrane tension

        Returns:
            -------
            None
        """
        self.gamma_bar = gamma_bar
        self.sigma_bar = sigma_bar

class ParticleGeometry:
    """
    A class to represent the particle.

    Attributes:
        ----------
        semi_minor_axis: float
            semi-minor axis of the elliptic particle
        semi_major_axis: float
            semi-major axis of the elliptic particle
        sampling_points_circle: int
            number of points to describe a circular particle with radius 1
        f: float
            wrapping degree
        theta: float
            angle used as polar coordinate to define the ellipse, see figure * of readme

    Methods:
        -------
        define_particle_geometry_variables(self, f):
            Returns useful geometrical parameters to perform further calculations
        compute_r_coordinate(self, f, theta):
            Returns the value of the coordinate r given the position on the ellipse
        compute_z_coordinate(self, f, theta):
            Returns the value of the coordinate z given the position on the ellipse
        get_alpha_angle(self, f):
            Returns the value of the alpha angle (see figure *)
        compute_psi1_psi3_angles(self, f):
            Returns the curvature angles in the particle (see figure *)
        get_squared_dpsi_region3(self, f):
            Returns the values of dpsi3**2

    """
    def __init__(self,
                 r_bar,
                 particle_perimeter,
                 sampling_points_circle):
        """
        Constructs all the necessary attributes for the particle object.

        Parameters:
            ----------
            r_bar: float
                particle's aspect ratio
            particle_perimeter: float
                particle's perimeter
            sampling_points_circle: int
                number of points to describe a circular particle with radius semi_major axis.

        Returns:
            -------
            None
        """
        self.sampling_points_circle = sampling_points_circle
        self.r_bar = r_bar
        self.particle_perimeter = particle_perimeter

        # computes the semi-minor and semi-major axes lengths of the elliptic particle
        # using Ramanujan's formula [2]

        h = ((self.r_bar-1)/(self.r_bar+1))**2
        self.semi_minor_axis = self.particle_perimeter / \
                               (pi*(1 + self.r_bar)*(1+ 3*h/(10 + sqrt(4-3*h))))
        self.semi_major_axis = self.semi_minor_axis * self.r_bar
        self.effective_radius = self.particle_perimeter / (2 * pi)

        # amount of points to sample the particle.
        self.sampling_points_ellipse = self.sampling_points_circle * self.particle_perimeter / \
                                      (2*pi*self.semi_major_axis)

    @lru_cache(maxsize=10)
    def define_particle_geometry_variables(self, f):
        """
        Defines all the necessary variables to describe the elliptic particle
        for a given wrapping degree f.

        Parameters:
            ----------
            f: float
                wrapping degree (between 0 and 1)

        Returns:
            -------
            beta: float
                trigonometric angle at intersection between regions 1, 2r and 3 (see figure *)
            beta_left: float
                trigonometric angle at intersection between regions 1, 2l and 3 (see figure *)
            theta_list_region1: array
                trigonomic angle theta into the region 1 (see figure *)
            theta_list_region3: array
                trigonomic angle theta into the region 3 (see figure *)
            l1: float
                arclength of the region 1 (see figure *)
            l3: float
                arclength of the region 3 (see figure *)
            s_list_region1: array
                sampling of the arclength of the region 1 (see figure *)
            s_list_region3: array
                sampling of the arclength of the region 3 (see figure *)

        """
        beta = pi*f + 1.5 * pi
        beta_left = 3*pi - beta
        n3 = int(f * self.sampling_points_ellipse)
        n1 = int((1-f) * self.sampling_points_ellipse)
        theta_list_region1 = np.linspace(beta, 2*pi + beta_left, n1)
        theta_list_region3 = np.linspace(beta_left, beta, n3)
        l1 = self.particle_perimeter * (1 - f)
        l3 = self.particle_perimeter * f
        s_list_region1 = np.linspace(0, l1, n1)
        s_list_region3 = np.linspace(0, l3, n3)
        return beta, beta_left, theta_list_region1, theta_list_region3, s_list_region3, l3, s_list_region1, l1

    @lru_cache(maxsize=128)
    def compute_r_coordinate(self, f, theta):
        """
        Computes the value of the coordinate r given the position on the ellipse,
            depicted by the angle theta, for a given wrapping degree f

        Parameters:
            ----------
            f: float
                wrapping degree
            theta: float
                trigonomic angle in the ellipse

        Returns:
            -------
            r_coordinate: float
                r coordinate (see figure *)
        """
        compute_x_coordinate = lambda t : self.semi_major_axis * cos(t)
        _, beta_left, _, _, _, _, _, _ = self.define_particle_geometry_variables(f)
        r_coordinate = compute_x_coordinate(theta) - compute_x_coordinate(beta_left)
        return r_coordinate

    @lru_cache(maxsize=128)
    def compute_z_coordinate(self, f, theta):
        """
        Computes the value of the coordinate z given the position on the ellipse,
            depicted by the angle theta, for a given wrapping degree f

        Parameters:
            ----------
            f: float
                wrapping degree
            theta: float
                trigonomic angle in the ellipse

        Returns:
            -------
            z_coordinate: float
                z coordinate (see figure *)
        """
        compute_y_coordinate = lambda t : self.semi_minor_axis * sin(t)
        _, beta_left, _, _, _, _, _, _ = self.define_particle_geometry_variables(f)
        z_coordinate = compute_y_coordinate(theta) - compute_y_coordinate(beta_left)
        return z_coordinate

    @lru_cache(maxsize=10)
    def get_alpha_angle(self, f):
        """
        Computes the value of the alpha angle (see figure *), for a given wrapping degree f

        Parameters:
            ----------
            f: float
                wrapping degree

        Returns:
            -------
            alpha: float
                curvature angle at intersection between regions 1, 2l and 3 (see figure *)

        """
        psi_list_region1, _ = self.compute_psi1_psi3_angles(f)
        alpha = psi_list_region1[0]
        return alpha

    @lru_cache(maxsize=10)
    def compute_psi1_psi3_angles(self, f):
        """
        Computes the curvature angles in the particle (see figure *),
            for a given wrapping degree f

        Parameters:
            ----------
            f: float
                wrapping degree

        Returns:
            -------
            psi_list_region1: list
                psi angle in region 1 (see figure *)
            psi_list_region3: list
                psi angle in region 3 (see figure *)

        """
        beta, beta_left, theta_list_region1, theta_list_region3, s_list_region3, _, s_list_region1, _ \
            = self.define_particle_geometry_variables(f)
        x_bl = self.semi_major_axis * cos(beta_left)

        def compute_psi_from_r_z(theta):
            """
            Computes the curvature angle given the position in the ellipse,
                depicted by theta (see figure *), using r and z coordinates

            Parameters:
                ----------
                theta: float
                    trigonomic angle in the ellipse

            Returns:
                -------
                delta: float
                    angle between the tangent to the particle and horizontal (see figure *)

            """
            r_elli = self.compute_r_coordinate(f, theta)
            z_elli = self.compute_z_coordinate(f, theta)

            def compute_tangent_to_ellipse_at_rtan_and_theta(r_tan):
                """
                Computes the position of the tangent to the particle (z with respect to r)

                Parameters:
                    ----------
                    r_tan: float
                        r coordinates where the tangent equation is evaluated

                Returns:
                    -------
                    z_tan: float
                        z coordinate of the tangent at r = r_tan

                """
                r_elli = self.compute_r_coordinate(f, theta)
                z_elli = self.compute_z_coordinate(f, theta)
                # managing possible singularities
                if r_elli == (-x_bl - self.semi_major_axis):
                    r_elli = r_elli + 0.01 * self.semi_major_axis
                elif r_elli == (-x_bl + self.semi_major_axis):
                    r_elli = r_elli - 0.01 * particle.semi_major_axis
                slope1 = 1
                """
                Returns:
                    Eq of the tangent to the ellipse at theta = beta
                """
                # depending on the side of the particle, the tangent's slope is positive or negative
                if theta > pi:
                    slope1 = -1
                dz = -self.semi_minor_axis * (r_elli + x_bl) / (self.semi_major_axis **2) / \
                    sqrt(1 - ((r_elli + x_bl)/self.semi_major_axis)**2)
                z_tan = z_elli  + slope1 * dz * (r_tan - r_elli)
                return z_tan

            r_tan_list = np.linspace(0.5 * r_elli - 0.5, 1.5 * r_elli + 0.5, 10)
            r1 = max(r_tan_list)
            z1 = compute_tangent_to_ellipse_at_rtan_and_theta(r1)
            delta = atan(abs((z1 - z_elli)/(r1 - r_elli)))
            return delta

        delta_list_region1 = [compute_psi_from_r_z(t) for t in theta_list_region1]
        delta_list_region3 = [compute_psi_from_r_z(t) for t in theta_list_region3]
        psi_list_region1 = np.zeros_like(s_list_region1)
        psi_list_region3 = np.zeros_like(s_list_region3)

        if f < 0.5:  # psi angle is defined differently depending on the position on the particle
            for i in range(len(s_list_region3)):
                theta = theta_list_region3[i]
                delta = delta_list_region3[i]
                if theta < 1.5*pi:
                    psi = 2*pi - delta
                elif theta == 1.5*pi:
                    psi = 2*pi
                elif theta <= beta:
                    psi = 2*pi + delta
                psi_list_region3[i] = psi
            for i in range(len(s_list_region1)):
                theta = theta_list_region1[i]
                delta = delta_list_region1[i]
                if theta <= 2*pi:
                    psi = delta
                elif theta <= 2*pi + pi/2:
                    psi = pi - delta
                elif theta <= 3*pi:
                    psi = pi + delta
                else:
                    psi = 2*pi - delta
                psi_list_region1[i] = psi
        else:
            for i in range(len(s_list_region3)):
                theta = theta_list_region3[i]
                delta = delta_list_region3[i]
                if theta <pi:
                    psi = pi + delta
                elif theta == pi:
                    psi = 1.5*pi
                elif theta < 1.5*pi:
                    psi = 2*pi - delta
                elif theta == 1.5*pi:
                    psi = 2*pi
                elif theta < 2*pi:
                    psi = 2*pi + delta
                elif theta == 2*pi:
                    psi = 2.5*pi
                elif theta <=beta:
                    psi = 3*pi - delta
                psi_list_region3[i] = psi

            for i in range(len(s_list_region1)):
                theta = theta_list_region1[i]
                delta = delta_list_region1[i]
                if theta < 2.5*pi:
                    psi = pi - delta
                elif theta == 2.5 * pi:
                    theta = pi
                elif theta <= 2*pi + beta_left:
                    psi = pi + delta
                psi_list_region1[i] = psi

        return psi_list_region1, psi_list_region3

    def get_squared_dpsi_region3(self, f):
        """
        Computes the values of dpsi3**2, necessary to evaluate the bending energy
            of the region 3, for a given wrapping degree f

        Parameters:
            ----------
            f: float
                wrapping degree

        Returns:
            -------
            squared_dpsi_list_region3: list
                dpsi angle power 2 in region 3 (see figure *)

        """
        _, _, _, _, s_list_region3, _, _, _ = self.define_particle_geometry_variables(f)
        _, psi_list_region3 = self.compute_psi1_psi3_angles(f)
        ds = s_list_region3[1] - s_list_region3[0]
        # computes the derivative of psi in region 3 using finite differences method. The value of
        # ds was set after a convergence study.
        dpsi3_list_region3 = [(psi_list_region3[i+1] - psi_list_region3[i])/ds
                               for i in range(0, len(psi_list_region3)-1)]
        squared_dpsi_list_region3 = [p**2 for p in dpsi3_list_region3]
        return squared_dpsi_list_region3

class MembraneGeometry:
    """
    A class to represent the membrane object.

    Attributes:
        ----------
        particle: class
            ParticleGeometry class
        mechanics: class
            MechanicalProperties class
        sampling_points_membrane: int
            number of points to sample the regions 2r and 2l

    Methods:
        -------
        compute_r2r_r2l_z2r_z2l_from_analytic_expression(self, f, particle, mechanics):
            Returns r and z coordinates in the regions 2r and 2l

    """
    def __init__(self, particle, sampling_points_membrane):
        """
        Constructs all the necessary attributes for the membrane object.

        Parameters:
            ----------
            particle: class
                ParticleGeometry class
            sampling_points_membrane: int
                number of points to sample the regions 2r and 2l

        Returns:
            -------
            None
        """
        self.sampling_points_membrane = sampling_points_membrane
        self.l2 = 20 * particle.effective_radius
        S2a = np.linspace(0, (self.l2 / 2), int((0.8*self.sampling_points_membrane) + 1))
        S2b = np.linspace(1.2 * self.l2 / 2, self.l2, int((0.2 * self.sampling_points_membrane)))
        self.S2 = np.concatenate((S2a, S2b), axis=None)

    @lru_cache(maxsize=10)
    def compute_r2r_r2l_z2r_z2l_from_analytic_expression(self, f, particle, mechanics):
        """
        Computes the r and z coordinates to describe the regions 2r and 2l,
            for a given wrapping degree f

        Parameters:
            ----------
            f: float
                wrapping degree
            particle: class
                ParticleGeometry class
            mechanics: class
                MechanicalProperties class
        Returns:
            -------
            r2r: list
                r coordinate in the region 2r
            z2r: list
                z coordinate in the region 2l
            r2l: list
                r coordinate in the region 2r
            z2l: list
                z coordinate in the region 2l
        """
        _, _, _, theta_list_region3, _, _, _, _ = particle.define_particle_geometry_variables(f)
        alpha = particle.get_alpha_angle(f)
        r2r_0 = particle.compute_r_coordinate(f, theta_list_region3[-1])
        z2r_0 = particle.compute_z_coordinate(f, theta_list_region3[-1])
        r2r = np.zeros_like(self.S2)
        z2r = np.zeros_like(self.S2)
        sigma = mechanics.sigma_bar

        for i in range(1, len(self.S2)):
            s = self.S2[i]
            r = r2r_0 + s - sqrt(2 / sigma) * (1 - cos(alpha)) / \
                (coth(s * sqrt(0.5*sigma)) + cos(alpha * 0.5))
            z = z2r_0 + sqrt(8 / sigma) * sin(0.5*alpha) * \
                (1 - (csch(s *sqrt(0.5*sigma))) / (coth(s* sqrt(0.5*sigma)) + cos(0.5*alpha)))
            r2r[i] = r
            z2r[i] = z

        r2r[0] = r2r_0
        z2r[0] = z2r_0
        r2l = np.array([r2r[0] - r2r[s] for s in range(len(self.S2))])
        z2l = z2r
        return r2r, z2r, r2l, z2l

class Wrapping:
    """
    A class to represent the wrapping of the particle.

    Attributes:
        ----------
        wrapíng_list: list
            list of wrapping degrees at which the system is evaluated

    """
    def __init__(self, wrapping_list):
        """
        Constructs all the necessary attributes for the membrane object.

        Parameters:
            ----------
            wrapíng_list: list
                list of wrapping degrees at which the system is evaluated

        Returns:
            -------
            None
        """
        self.wrapping_list = wrapping_list

def compute_adimensioned_energy_variation(f, particle, mechanics, membrane):
    """
    Computes the adimensional energy variation between a given wrapping degree f
        and the initial state at wero wrapping, as expressed in [1]

    Parameters:
        ----------
        f: float
            wrapping degree
        particle: class
            ParticleGeometry class
        mechanics: class
            MechanicalProperties class
        membrane: class
            MembraneGeometry class

    Returns:
        -------
        adimensional_total_energy_variation: float
            computed value of the adimensional total energy variation
    """
    _, _, _, _, s_list_region3, l3, _, _ = particle.define_particle_geometry_variables(f)
    r2r, _, r2l, _ = membrane.compute_r2r_r2l_z2r_z2l_from_analytic_expression(f,
                                                                               particle,
                                                                               mechanics)

    def bending_energy3():
        """
        Args:
            f: wrapping degree
            particle: Parameters class
        Returns:
            float - the variation of the bending energy in the region 3 (bending of the membrane)
            between the state at wrapping f and the initial state: E(f) - E(0)
        """
        dpsi3_list_region3 = particle.get_squared_dpsi_region3(f)
        energy = scipy.integrate.simps(dpsi3_list_region3, s_list_region3[0:-1])
        return energy

    def bending_energy2d():
        """
        Args:
            f: wrapping degree
            particle: Parameters class
        Returns:
            float - the variation of the bending energy in the region 2d (free membrane)
            between the state at wrapping f and the initial state: E(f) - E(0)
        """
        a = particle.get_alpha_angle(f)
        t = tan(0.25 * a)
        t2 = t**2
        b = sqrt(mechanics.sigma_bar / 2)
        E2 = - 8* b / particle.effective_radius * t2 * \
             ((1/(t2 + exp(2* b *membrane.l2 / particle.effective_radius))) - (1/(t2 + 1)))
        return E2

    Eb3 =  0.25 * particle.effective_radius * bending_energy3()
    Eb2d = 0.25 * particle.effective_radius * bending_energy2d()
    Eb2g = Eb2d
    Eb = Eb2d + Eb2g + Eb3
    Eadh = - mechanics.gamma_bar * l3 * 0.25 / particle.effective_radius
    Etens =  mechanics.sigma_bar * ( l3 + 2*membrane.l2 - (r2r[-1] - r2l[-1])) * \
             0.25 / particle.effective_radius
    adimensional_total_energy_variation = Eb + Eadh + Etens
    return adimensional_total_energy_variation

def plot_energy(particle, mechanics, membrane, wrapping):
    """
    Plots the evolution of the adimensional variation of energy during wrapping

    Parameters:
        ----------
        particle: class
            ParticleGeometry class
        mechanics: class
            MechanicalProperties class
        membrane: class
            MembraneGeometry class

    Returns:
        -------
        None
    """
    energy_list = np.array([compute_adimensioned_energy_variation(f,
                                                                  particle,
                                                                  mechanics,
                                                                  membrane) for f in wrapping.wrapping_list])
    plt.figure()
    plt.plot(wrapping.wrapping_list,
             energy_list,
             '-k',
             label = '$\overline{r} = $' + str(np.round(particle.r_bar, 2)) +
                     ' ; $\overline{\gamma} = $' + str(mechanics.gamma_bar) +
                     ' ; $\overline{\sigma} = $' + str(mechanics.sigma_bar))
    plt.xlabel('wrapping degree f [-]')
    plt.ylabel('$\Delta E [-]$')
    plt.title('$\Delta E(f)$')
    plt.legend()
    plt.xlim((0, 1))

def identify_wrapping_phase(particle, mechanics, membrane, wrapping):
    """
    Identifies the wrapping phase following the process introduced in [1]

    Parameters:
        ----------
        particle: class
            ParticleGeometry object
        mechanics: class
            MechanicalProperties object
        membrane: class
            MembraneGeometry object

    Returns:
        -------
        wrapping_phase_number: float
            phase number (1, 2 or 3)
        wrapping_phase: str
            the wrapping phase as an intelligible string
    """
    def get_local_energy_minima():
        """
        Identifies the local minima of energy variation with their location

        Parameters:
            ----------
            None

        Returns:
            -------
            energy_list:list
                list of adimensional energy variation during wrapping
            min_energy_list: list
                list of local minimna of adimensional energy with respect to wrapping degree f
            f_min_energy_list: list
                list of wrapping degree values where local minima of adimensional energy are reached
        """
        energy_list = np.array([compute_adimensioned_energy_variation(f,
                                                                      particle,
                                                                      mechanics,
                                                                      membrane) for f in wrapping.wrapping_list])
        min_energy_index_list = scipy.signal.argrelextrema(energy_list, np.less)

        # check if the minimum is reached for wrapping.wrapping_list[-1]
        if energy_list[-1] <  energy_list[-2]:
            min_energy_index_list = np.concatenate((min_energy_index_list,
                                                    np.array([-1])), axis=None)

        # check if the minimum is reached for wrapping.wrapping_list[0]
        if energy_list[0] < energy_list[1]:
            min_energy_index_list = np.concatenate((np.array(wrapping.wrapping_list[0]),
                                                    min_energy_index_list), axis=None)

        min_energy_list = [energy_list[int(k)] for k in min_energy_index_list]
        f_min_energy_list = [wrapping.wrapping_list[int(k)] for k in min_energy_index_list]

        return energy_list, min_energy_list, f_min_energy_list

    _, _, f_min_energy_list = get_local_energy_minima()
    f_eq = f_min_energy_list[0]
    wrapping_phase_number = 0
    wrapping_phase = '0'

    if f_eq < 0.2:  # check if wrapping phase is phase 1, according to [1]
        wrapping_phase_number = 1
        wrapping_phase = 'no wrapping'
    else:
        r2r, _, r2l, _ = membrane.compute_r2r_r2l_z2r_z2l_from_analytic_expression(f_eq,
                                                                                particle,
                                                                                mechanics)
        intersection_membrane = min(r2r) - max(r2l)
        wrapping_phase_number = 3 if intersection_membrane < 0 else 2
        wrapping_phase = 'full wrapping' if intersection_membrane < 0 else "partial wrapping"

    def plot_geometry(f, particle, mechanics, membrane):
        """
        Plots the morphology of the interaction between the cell and the particle,
            for a given wrapping degree f

        Parameters:
            ----------
            f: float
                wrapping degree
            particle: class
                ParticleGeometry class
            mechanics: class
                MechanicalProperties class
            membrane: class
                MembraneGeometry class

        Returns:
            -------
            None
        """
        _, _, theta_list_region1, theta_list_region3, _, _, _, _ = \
        particle.define_particle_geometry_variables(f)
        r1 = [particle.compute_r_coordinate(f, theta) for theta in theta_list_region1]
        r3 = [particle.compute_r_coordinate(f, theta) for theta in theta_list_region3]
        z1 = [particle.compute_z_coordinate(f, theta) for theta in theta_list_region1]
        z3 = [particle.compute_z_coordinate(f, theta) for theta in theta_list_region3]
        r2r, z2r, r2l, z2l = \
        membrane.compute_r2r_r2l_z2r_z2l_from_analytic_expression(f, particle, mechanics)
        plt.figure()
        plt.plot(r1, z1, '-b', label = 'region 1')
        plt.plot(r2r, z2r, '-k', label = 'region 2')
        plt.plot(r2l, z2l, '-k')
        plt.plot(r3, z3, '-r', label = 'region 3')
        plt.xlabel('r (x100) (nm)')
        plt.ylabel('z (x100) (nm)')
        plt.ylim((-4, 4))
        plt.title('wrapping morphology')
        plt.gca().set_aspect('equal' , adjustable='box')
        plt.legend()

    plot_geometry(f_eq, particle, mechanics, membrane)
    plt.title('wrapping morphology at equilibrium: ' + wrapping_phase)

    return wrapping_phase_number, wrapping_phase

def parse_arguments():
    """
    Parses arguments to run the code in terminal

    Parameters:
        ----------
        None

    Returns:
        -------
        args: class
            #TODO complete here
    """
    parser = argparse.ArgumentParser()
    parser.add_argument('-g', '--gamma_bar', required=False, default=10.0, type=float,
        help='adimensional lineic adhesion between the membrane and the particle. Default value = 10.')
    parser.add_argument('-s', '--sigma_bar', required=False, default=2.0, type=float,
        help='adimensional membrane tension. Default value = 2.')
    parser.add_argument('-r', '--r_bar', required=False, default=1.0, type=float,
        help='particle aspect ratio. Default value = 1.')
    parser.add_argument('-p', '--particle_perimeter', required=False, default=2*pi, type=float,
        help='particle perimeter. Default value = 2pi.')
    args = parser.parse_args()
    return args

if __name__ == "__main__":
    args = parse_arguments()

    particle = ParticleGeometry(r_bar = args.r_bar,
                                particle_perimeter = args.particle_perimeter,
                                sampling_points_circle=300)

    mechanics = MechanicalProperties(gamma_bar=args.gamma_bar,
                                     sigma_bar=args.sigma_bar)

    membrane = MembraneGeometry(particle,
                                sampling_points_membrane=100)

    wrapping = Wrapping(wrapping_list=np.arange(0.03, 0.97, 0.003125))

    plot_energy(particle, mechanics, membrane, wrapping)
    wrapping_phase_number, wrapping_phase = identify_wrapping_phase(particle, mechanics, membrane, wrapping)
    print('wrapping phase at equilibrium: ', wrapping_phase)
    plt.show()

reason: import numpy,import scipy,from mpmath


text:
import matplotlib
matplotlib.use('Agg')  # noqa: E402
import matplotlib.pyplot as plt
import numpy as np
import pytest
import skimage.data

from correlateim import transform

@pytest.fixture(scope="session")
def source_coords():
    src = np.array(
        [[83.09970486, 73.68900001],
         [70.22116071, 490.09526094],
         [518.39449721, 525.29661496],
         [356.9834105, 202.47444154],
         [212.74371599, 351.0069841]])
    return src

@pytest.fixture(scope="session")
def destination_coords():
    dst = np.array(
        [[79.66542642, 69.48044789],
         [137.1895903, 440.38251948],
         [533.84875019, 381.14121638],
         [340.67058791, 128.72175099],
         [241.07651313, 290.99140731]])
    return dst

@pytest.fixture(scope="session")
def example_affine_matrix():
    r = 0.12
    c, s = np.cos(r), np.sin(r)
    matrix_transform = np.array([[c, -s, 0],
                                 [s, c, 50],
                                 [0, 0, 1]])
    return matrix_transform

@pytest.mark.mpl_image_compare
def test_apply_transform(example_affine_matrix):
    image = skimage.data.astronaut()
    output = transform.apply_transform(image, example_affine_matrix)
    fig = plt.figure()
    ax = fig.add_subplot(1, 1, 1)
    ax.axis('off')
    ax.imshow(output)
    return fig

def test_apply_transform_dimension_mismatch(example_affine_matrix):
    image = np.random.random((10, 10, 3))
    with pytest.raises(ValueError):
        transform.apply_transform(image, example_affine_matrix,
                                  multichannel=False)

@pytest.mark.mpl_image_compare
def test_apply_transform_multichannel_false(example_affine_matrix):
    image = skimage.img_as_float(skimage.data.camera())
    output = transform.apply_transform(image, example_affine_matrix,
                                       multichannel=False)
    fig = plt.figure()
    ax = fig.add_subplot(1, 1, 1)
    ax.axis('off')
    ax.imshow(output)
    return fig

@pytest.mark.mpl_image_compare
def test_apply_transform_inverse(example_affine_matrix):
    image = skimage.data.astronaut()
    output = transform.apply_transform(image, example_affine_matrix,
                                       inverse=False)
    fig = plt.figure()
    ax = fig.add_subplot(1, 1, 1)
    ax.axis('off')
    ax.imshow(output)
    return fig

def test_calculate_transform(source_coords, destination_coords):
    output = transform.calculate_transform(source_coords, destination_coords)
    expected_output = np.array(
        [[0.87231142, 0.16466797, -4.19019986],
         [-0.20180383, 0.88444066, 21.69926004],
         [0., 0., 1.]])
    assert np.allclose(output, expected_output)

def test_point_coords(source_coords, destination_coords):
    example_matched_points_dict = [
        {'point_id': 1,
         'img1_x': 73.68900001398778,
         'img1_y': 83.09970485795407,
         'img2_x': 69.48044788854543,
         'img2_y': 79.66542641731121},
        {'point_id': 2,
         'img1_x': 490.0952609419367,
         'img1_y': 70.2211607055433,
         'img2_x': 440.3825194779763,
         'img2_y': 137.18959029807945},
        {'point_id': 3,
         'img1_x': 525.2966149585263,
         'img1_y': 518.3944972094389,
         'img2_x': 381.14121637688663,
         'img2_y': 533.8487501923319},
        {'point_id': 4,
         'img1_x': 202.4744415380957,
         'img1_y': 356.9834104992237,
         'img2_x': 128.7217509896351,
         'img2_y': 340.67058790617},
        {'point_id': 5,
         'img1_x': 351.0069840959002,
         'img1_y': 212.74371599222275,
         'img2_x': 290.991407310011,
         'img2_y': 241.07651312752648}]
    src, dst = transform.point_coords(example_matched_points_dict)
    assert np.allclose(src, source_coords)
    assert np.allclose(dst, destination_coords)

ext: py

max_stars_repo_name: DeMarcoLab/correlateim

max_issues_repo_name: DeMarcoLab/correlateim

max_forks_repo_name: DeMarcoLab/correlateim

reason: import numpy


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed May 22 11:21:00 2019

@author: shlomi
"""

def concat_to_years():
    from pathlib import Path
    import numpy as np
    import xarray as xr
    import pandas as pd
    data11 = Path('/data11/ziskin/MERRA2/')
    years = np.arange(2001, 2019)
    for year in years:
        ds_year = []
        print('proccessing year {}:'.format(year))
        part = '100'
        if year > 1991:
            part = '200'
        if year > 2000:
            part = '300'
        if year > 2010:
            part = '400'
        year_date = pd.date_range(
            str(year) + '-01-01',
            str(year) + '-12-31',
            freq='D')
        for date in year_date:
            filename = 'MERRA2_' + part + '.inst6_3d_ana_Nv.' + \
                date.strftime('%Y%m%d') + '.SUB.nc'
            try:
                merra2 = xr.open_dataset(data11.as_posix() + '/' + filename)
            except FileNotFoundError:
                print('{} not found...skipping'.format(filename))
                continue
            ds_year.append(merra2)
        ds = xr.concat(ds_year, 'time')
        comp = dict(zlib=True, complevel=9)  # best compression
        encoding = {var: comp for var in ds.data_vars}
        file_to_save = 'MERRA2_T_cold_' + str(year) + '.nc'
        ds.to_netcdf(data11 / file_to_save, 'w', encoding=encoding)
    return

def concat_to_one_big_file():
    from pathlib import Path
    import xarray as xr
    data11 = Path('/data11/ziskin/MERRA2/')
    from dask.diagnostics import ProgressBar
    merra2 = xr.open_mfdataset(data11.as_posix() + '/' 'MERRA2_T_cold_*.nc')
    comp = dict(zlib=True, complevel=9)  # best compression
    encoding = {var: comp for var in merra2.data_vars}
    file_to_save = 'MERRA2_Tcold_BIG.nc'
    with ProgressBar():
        merra2.to_netcdf(data11 / file_to_save, 'w', encoding=encoding)
    return

concat_to_one_big_file()


text:
# Authors: Jose C. Garcia Alanis <alanis.jcg@gmail.com>
#
# License: BSD-3-Clause

import os
import os.path as op
import shutil
import zipfile
from sys import stdout

import numpy as np

from ...channels import make_standard_montage
from ...epochs import EpochsArray
from ...io.meas_info import create_info
from ...utils import _fetch_file, _check_pandas_installed, verbose
from ..utils import _get_path, _do_path_update

# root url for LIMO files
root_url = 'https://files.de-1.osf.io/v1/resources/52rea/providers/osfstorage/'

# subject identifier
subject_ids = {'S1': '5cde823c8d6e050018595862',
               'S2': '5cde825e23fec40017e0561a',
               'S3': '5cf7eedee650a2001ad560f2',
               'S4': '5cf7eee7d4c7d700193defcb',
               'S5': '5cf7eeece650a20017d5b153',

               'S6': '5cf8300fe650a20018d59cef',
               'S7': '5cf83018a542b8001bc7c75f',
               'S8': '5cf8301ea542b8001ac7cc47',
               'S9': '5cf830243a4d9500178a692b',
               'S10': '5cf83029e650a20017d600b1',

               'S11': '5cf834bfa542b8001bc7cae0',
               'S12': '5cf834c53a4d9500188a6311',
               'S13': '5cf834caa542b8001cc8149b',
               'S14': '5cf834cf3a4d9500178a6c6c',
               'S15': '5cf834d63a4d9500168ae5d6',

               'S16': '5cf834dbe650a20018d5a123',
               'S17': '5cf834e23a4d9500198a911f',
               'S18': '5cf834e73a4d9500198a9122'}

@verbose
def data_path(subject, path=None, force_update=False, update_path=None,
              verbose=None):
    """Get path to local copy of LIMO dataset URL.

    This is a low-level function useful for getting a local copy of the
    remote LIMO dataset :footcite:`Rousselet2016`. The complete dataset is
    available at datashare.is.ed.ac.uk/.

    Parameters
    ----------
    subject : int
        Subject to download. Must be of class ìnt in the range from 1 to 18.
    path : None | str
        Location of where to look for the LIMO data storing directory.
        If None, the environment variable or config parameter
        ``MNE_DATASETS_LIMO_PATH`` is used. If it doesn't exist, the
        "~/mne_data" directory is used. If the LIMO dataset
        is not found under the given path, the data
        will be automatically downloaded to the specified folder.
    force_update : bool
        Force update of the dataset even if a local copy exists.
    update_path : bool | None
        If True, set the MNE_DATASETS_LIMO_PATH in mne-python
        config to the given path. If None, the user is prompted.
    %(verbose)s

    Returns
    -------
    path : str
        Local path to the given data file.

    Notes
    -----
    For example, one could do:

        >>> from mne.datasets import limo
        >>> limo.data_path(subject=1, path=os.getenv('HOME') + '/datasets') # doctest:+SKIP

    This would download the LIMO data file to the 'datasets' folder,
    and prompt the user to save the 'datasets' path to the mne-python config,
    if it isn't there already.

    References
    ----------
    .. footbibliography::
    """  # noqa: E501
    # set destination path for download
    key = 'MNE_DATASETS_LIMO_PATH'
    name = 'LIMO'
    path = _get_path(path, key, name)
    limo_dir = op.join(path, 'MNE-limo-data')
    subject_id = 'S%s' % subject
    destination = op.join(limo_dir, '%s.zip') % subject_id

    # url for subject in question
    url = root_url + subject_ids[subject_id] + '/?zip='

    # check if LIMO directory exists; update if desired
    if not op.isdir(limo_dir) or force_update:
        if op.isdir(limo_dir):
            shutil.rmtree(limo_dir)
        if not op.isdir(limo_dir):
            os.makedirs(limo_dir)

    # check if subject in question exists
    if not op.isdir(op.join(limo_dir, subject_id)):
        os.makedirs(op.join(limo_dir, subject_id))
        _fetch_file(url, destination, print_destination=False)

        # check if download is a zip-folder
        if any(group.endswith(".zip") for group in op.splitext(destination)):
            if not op.isdir(op.join(limo_dir, subject_id)):
                os.makedirs(op.join(limo_dir, subject_id))
            with zipfile.ZipFile(destination) as z1:
                files = [op.join(limo_dir, file) for file in z1.namelist()]
                stdout.write('Decompressing %g files from\n'
                             '"%s" ...' % (len(files), destination))
                z1.extractall(op.join(limo_dir, subject_id))
                stdout.write(' [done]\n')
                z1.close()
                os.remove(destination)

    # update path if desired
    _do_path_update(path, update_path, key, name)

    return limo_dir

@verbose
def load_data(subject, path=None, force_update=False, update_path=None,
              verbose=None):
    """Fetch subjects epochs data for the LIMO data set.

    Parameters
    ----------
    subject : int
        Subject to use. Must be of class ìnt in the range from 1 to 18.
    path : str
        Location of where to look for the LIMO data.
        If None, the environment variable or config parameter
        ``MNE_DATASETS_LIMO_PATH`` is used. If it doesn't exist, the
        "~/mne_data" directory is used.
    force_update : bool
        Force update of the dataset even if a local copy exists.
    update_path : bool | None
        If True, set the MNE_DATASETS_LIMO_PATH in mne-python
        config to the given path. If None, the user is prompted.
    %(verbose)s

    Returns
    -------
    epochs : instance of Epochs
        The epochs.
    """  # noqa: E501
    pd = _check_pandas_installed()
    from scipy.io import loadmat

    # subject in question
    if isinstance(subject, int) and 1 <= subject <= 18:
        subj = 'S%i' % subject
    else:
        raise ValueError('subject must be an int in the range from 1 to 18')

    # set limo path, download and decompress files if not found
    limo_path = data_path(subject, path, force_update, update_path)

    # -- 1) import .mat files
    # epochs info
    fname_info = op.join(limo_path, subj, 'LIMO.mat')
    data_info = loadmat(fname_info)
    # number of epochs per condition
    design = data_info['LIMO']['design'][0][0]['X'][0][0]
    data_info = data_info['LIMO']['data'][0][0][0][0]
    # epochs data
    fname_eeg = op.join(limo_path, subj, 'Yr.mat')
    data = loadmat(fname_eeg)

    # -- 2) get epochs information from structure
    # sampling rate
    sfreq = data_info['sampling_rate'][0][0]
    # tmin and tmax
    tmin = data_info['start'][0][0]
    # create events matrix
    sample = np.arange(len(design))
    prev_id = np.zeros(len(design))
    ev_id = design[:, 1]
    events = np.array([sample, prev_id, ev_id]).astype(int).T
    # event ids, such that Face B == 1
    event_id = {'Face/A': 0, 'Face/B': 1}

    # -- 3) extract channel labels from LIMO structure
    # get individual labels
    labels = data_info['chanlocs']['labels']
    labels = [label for label, *_ in labels[0]]
    # get montage
    montage = make_standard_montage('biosemi128')
    # add external electrodes (e.g., eogs)
    ch_names = montage.ch_names + ['EXG1', 'EXG2', 'EXG3', 'EXG4']
    # match individual labels to labels in montage
    found_inds = [ind for ind, name in enumerate(ch_names) if name in labels]
    missing_chans = [name for name in ch_names if name not in labels]
    assert labels == [ch_names[ind] for ind in found_inds]

    # -- 4) extract data from subjects Yr structure
    # data is stored as channels x time points x epochs
    # data['Yr'].shape  # <-- see here
    # transpose to epochs x channels time points
    data = np.transpose(data['Yr'], (2, 0, 1))
    # initialize data in expected order
    temp_data = np.empty((data.shape[0], len(ch_names), data.shape[2]))
    # copy over the non-missing data
    for source, target in enumerate(found_inds):
        # avoid copy when fancy indexing
        temp_data[:, target, :] = data[:, source, :]
    # data to V (to match MNE's format)
    data = temp_data / 1e6
    # create list containing channel types
    types = ["eog" if ch.startswith("EXG") else "eeg" for ch in ch_names]

    # -- 5) Create custom info for mne epochs structure
    # create info
    info = create_info(ch_names, sfreq, types).set_montage(montage)
    # get faces and noise variables from design matrix
    event_list = list(events[:, 2])
    faces = ['B' if event else 'A' for event in event_list]
    noise = list(design[:, 2])
    # create epochs metadata
    metadata = {'face': faces, 'phase-coherence': noise}
    metadata = pd.DataFrame(metadata)

    # -- 6) Create custom epochs array
    epochs = EpochsArray(data, info, events, tmin, event_id, metadata=metadata)
    epochs.info['bads'] = missing_chans  # missing channels are marked as bad.

    return epochs

max_stars_repo_name: vagechirkov/mne-python

max_issues_repo_name: jaromilfrossard/mne-python

max_forks_repo_name: jaromilfrossard/mne-python

reason: import numpy,from scipy


text:
#!/usr/bin/env python3

import math
import matplotlib.pyplot as plt
import numpy as np

plt.rcParams["figure.dpi"] = 140
x = np.linspace(2, 10, 1000000)

def f(input):
    return math.sqrt(input - 2) + 3

vfunc = np.vectorize(f)
y = vfunc(x)

fig, ax = plt.subplots()

ax.spines["left"].set_position("zero")
ax.spines["bottom"].set_position("zero")
ax.spines["right"].set_color("none")
ax.spines["top"].set_color("none")
ax.xaxis.set_ticks_position("bottom")
ax.yaxis.set_ticks_position("left")

ax.set_aspect("equal")
plt.xlim(-3, 10)
plt.ylim(-1, 10)

ax.plot(x, y, 'r-', linewidth=2, label='$\sqrt{x-2} + 3$', alpha=0.6)

ax.set_title('Horizental Line text')
ax.legend(loc='best')

fig.savefig('../out/inverse.png')

max_stars_repo_name: GUC-Fractals/math-curriculum

max_issues_repo_name: GUC-Fractals/math-curriculum

max_forks_repo_name: GUC-Fractals/math-curriculum


text:
def example(Simulator):
    from csdl import Model
    import csdl
    import numpy as np

    class ExampleTensorTensor(Model):
        def define(self):

            m = 3
            n = 4
            p = 5

            # Shape of the tensors
            ten_shape = (m, n, p)

            # Number of elements in the tensors
            num_ten_elements = np.prod(ten_shape)

            # Values for the two tensors
            ten1 = np.arange(num_ten_elements).reshape(ten_shape)
            ten2 = np.arange(num_ten_elements,
                             2 * num_ten_elements).reshape(ten_shape)

            # Adding tensors to csdl
            ten1 = self.declare_variable('ten1', val=ten1)
            ten2 = self.declare_variable('ten2', val=ten2)

            # Tensor-Tensor Inner Product specifying the first and last axes
            self.register_output(
                'TenTenInner',
                csdl.inner(ten1, ten2, axes=([0, 2], [0, 2])),
            )

    sim = Simulator(ExampleTensorTensor())
    sim.run()

    print('ten1', sim['ten1'].shape)
    print(sim['ten1'])
    print('ten2', sim['ten2'].shape)
    print(sim['ten2'])
    print('TenTenInner', sim['TenTenInner'].shape)
    print(sim['TenTenInner'])

    return sim

max_stars_repo_name: LSDOlab/csdl

max_issues_repo_name: LSDOlab/csdl

max_forks_repo_name: LSDOlab/csdl


text:
# -*- coding: utf-8 -*-
# File generated according to Generator/ClassesRef/Simulation/SliceModel.csv
# WARNING! All changes made in this file will be lost!
"""Method code available at https://github.com/Eomys/pyleecan/tree/master/pyleecan/Methods/Simulation/SliceModel
"""

from os import linesep
from sys import getsizeof
from logging import getLogger
from ._check import set_array, check_var, raise_
from ..Functions.get_logger import get_logger
from ..Functions.save import save
from ..Functions.copy import copy
from ..Functions.load import load_init_dict
from ..Functions.Load.import_class import import_class
from ._frozen import FrozenClass

# Import all class method
# Try/catch to remove unnecessary dependencies in unused method
try:
    from ..Methods.Simulation.SliceModel.get_distribution import get_distribution
except ImportError as error:
    get_distribution = error

try:
    from ..Methods.Simulation.SliceModel.get_data import get_data
except ImportError as error:
    get_data = error

try:
    from ..Methods.Simulation.SliceModel.plot import plot
except ImportError as error:
    plot = error

from numpy import array, array_equal
from ._check import InitUnKnowClassError

class SliceModel(FrozenClass):
    """Class to hande slices in magnetics model"""

    VERSION = 1

    # Check ImportError to remove unnecessary dependencies in unused method
    # cf Methods.Simulation.SliceModel.get_distribution
    if isinstance(get_distribution, ImportError):
        get_distribution = property(
            fget=lambda x: raise_(
                ImportError(
                    "Can't use SliceModel method get_distribution: "
                    + str(get_distribution)
                )
            )
        )
    else:
        get_distribution = get_distribution
    # cf Methods.Simulation.SliceModel.get_data
    if isinstance(get_data, ImportError):
        get_data = property(
            fget=lambda x: raise_(
                ImportError("Can't use SliceModel method get_data: " + str(get_data))
            )
        )
    else:
        get_data = get_data
    # cf Methods.Simulation.SliceModel.plot
    if isinstance(plot, ImportError):
        plot = property(
            fget=lambda x: raise_(
                ImportError("Can't use SliceModel method plot: " + str(plot))
            )
        )
    else:
        plot = plot
    # save and copy methods are available in all object
    save = save
    copy = copy
    # get_logger method is available in all object
    get_logger = get_logger

    def __init__(
        self,
        type_distribution=None,
        Nslices=5,
        z_list=None,
        angle_rotor=None,
        angle_stator=None,
        L=None,
        is_step=None,
        is_skew=None,
        init_dict=None,
        init_str=None,
    ):
        """Constructor of the class. Can be use in three ways :
        - __init__ (arg1 = 1, arg3 = 5) every parameters have name and default values
            for pyleecan type, -1 will call the default constructor
        - __init__ (init_dict = d) d must be a dictionary with property names as keys
        - __init__ (init_str = s) s must be a string
        s is the file path to load

        ndarray or list can be given for Vector and Matrix
        object or dict can be given for pyleecan Object"""

        if init_str is not None:  # Load from a file
            init_dict = load_init_dict(init_str)[1]
        if init_dict is not None:  # Initialisation by dict
            assert type(init_dict) is dict
            # Overwrite default value with init_dict content
            if "type_distribution" in list(init_dict.keys()):
                type_distribution = init_dict["type_distribution"]
            if "Nslices" in list(init_dict.keys()):
                Nslices = init_dict["Nslices"]
            if "z_list" in list(init_dict.keys()):
                z_list = init_dict["z_list"]
            if "angle_rotor" in list(init_dict.keys()):
                angle_rotor = init_dict["angle_rotor"]
            if "angle_stator" in list(init_dict.keys()):
                angle_stator = init_dict["angle_stator"]
            if "L" in list(init_dict.keys()):
                L = init_dict["L"]
            if "is_step" in list(init_dict.keys()):
                is_step = init_dict["is_step"]
            if "is_skew" in list(init_dict.keys()):
                is_skew = init_dict["is_skew"]
        # Set the properties (value check and convertion are done in setter)
        self.parent = None
        self.type_distribution = type_distribution
        self.Nslices = Nslices
        self.z_list = z_list
        self.angle_rotor = angle_rotor
        self.angle_stator = angle_stator
        self.L = L
        self.is_step = is_step
        self.is_skew = is_skew

        # The class is frozen, for now it's impossible to add new properties
        self._freeze()

    def __str__(self):
        """Convert this object in a readeable string (for print)"""

        SliceModel_str = ""
        if self.parent is None:
            SliceModel_str += "parent = None " + linesep
        else:
            SliceModel_str += "parent = " + str(type(self.parent)) + " object" + linesep
        SliceModel_str += (
            'type_distribution = "' + str(self.type_distribution) + '"' + linesep
        )
        SliceModel_str += "Nslices = " + str(self.Nslices) + linesep
        SliceModel_str += (
            "z_list = "
            + linesep
            + str(self.z_list).replace(linesep, linesep + "\t")
            + linesep
        )
        SliceModel_str += (
            "angle_rotor = "
            + linesep
            + str(self.angle_rotor).replace(linesep, linesep + "\t")
            + linesep
            + linesep
        )
        SliceModel_str += (
            "angle_stator = "
            + linesep
            + str(self.angle_stator).replace(linesep, linesep + "\t")
            + linesep
            + linesep
        )
        SliceModel_str += "L = " + str(self.L) + linesep
        SliceModel_str += "is_step = " + str(self.is_step) + linesep
        SliceModel_str += "is_skew = " + str(self.is_skew) + linesep
        return SliceModel_str

    def __eq__(self, other):
        """Compare two objects (skip parent)"""

        if type(other) != type(self):
            return False
        if other.type_distribution != self.type_distribution:
            return False
        if other.Nslices != self.Nslices:
            return False
        if other.z_list != self.z_list:
            return False
        if not array_equal(other.angle_rotor, self.angle_rotor):
            return False
        if not array_equal(other.angle_stator, self.angle_stator):
            return False
        if other.L != self.L:
            return False
        if other.is_step != self.is_step:
            return False
        if other.is_skew != self.is_skew:
            return False
        return True

    def compare(self, other, name="self", ignore_list=None):
        """Compare two objects and return list of differences"""

        if ignore_list is None:
            ignore_list = list()
        if type(other) != type(self):
            return ["type(" + name + ")"]
        diff_list = list()
        if other._type_distribution != self._type_distribution:
            diff_list.append(name + ".type_distribution")
        if other._Nslices != self._Nslices:
            diff_list.append(name + ".Nslices")
        if other._z_list != self._z_list:
            diff_list.append(name + ".z_list")
        if not array_equal(other.angle_rotor, self.angle_rotor):
            diff_list.append(name + ".angle_rotor")
        if not array_equal(other.angle_stator, self.angle_stator):
            diff_list.append(name + ".angle_stator")
        if other._L != self._L:
            diff_list.append(name + ".L")
        if other._is_step != self._is_step:
            diff_list.append(name + ".is_step")
        if other._is_skew != self._is_skew:
            diff_list.append(name + ".is_skew")
        # Filter ignore differences
        diff_list = list(filter(lambda x: x not in ignore_list, diff_list))
        return diff_list

    def __sizeof__(self):
        """Return the size in memory of the object (including all subobject)"""

        S = 0  # Full size of the object
        S += getsizeof(self.type_distribution)
        S += getsizeof(self.Nslices)
        if self.z_list is not None:
            for value in self.z_list:
                S += getsizeof(value)
        S += getsizeof(self.angle_rotor)
        S += getsizeof(self.angle_stator)
        S += getsizeof(self.L)
        S += getsizeof(self.is_step)
        S += getsizeof(self.is_skew)
        return S

    def as_dict(self, type_handle_ndarray=0, keep_function=False, **kwargs):
        """
        Convert this object in a json serializable dict (can be use in __init__).
        type_handle_ndarray: int
            How to handle ndarray (0: tolist, 1: copy, 2: nothing)
        keep_function : bool
            True to keep the function object, else return str
        Optional keyword input parameter is for internal use only
        and may prevent json serializability.
        """

        SliceModel_dict = dict()
        SliceModel_dict["type_distribution"] = self.type_distribution
        SliceModel_dict["Nslices"] = self.Nslices
        SliceModel_dict["z_list"] = (
            self.z_list.copy() if self.z_list is not None else None
        )
        if self.angle_rotor is None:
            SliceModel_dict["angle_rotor"] = None
        else:
            if type_handle_ndarray == 0:
                SliceModel_dict["angle_rotor"] = self.angle_rotor.tolist()
            elif type_handle_ndarray == 1:
                SliceModel_dict["angle_rotor"] = self.angle_rotor.copy()
            elif type_handle_ndarray == 2:
                SliceModel_dict["angle_rotor"] = self.angle_rotor
            else:
                raise Exception(
                    "Unknown type_handle_ndarray: " + str(type_handle_ndarray)
                )
        if self.angle_stator is None:
            SliceModel_dict["angle_stator"] = None
        else:
            if type_handle_ndarray == 0:
                SliceModel_dict["angle_stator"] = self.angle_stator.tolist()
            elif type_handle_ndarray == 1:
                SliceModel_dict["angle_stator"] = self.angle_stator.copy()
            elif type_handle_ndarray == 2:
                SliceModel_dict["angle_stator"] = self.angle_stator
            else:
                raise Exception(
                    "Unknown type_handle_ndarray: " + str(type_handle_ndarray)
                )
        SliceModel_dict["L"] = self.L
        SliceModel_dict["is_step"] = self.is_step
        SliceModel_dict["is_skew"] = self.is_skew
        # The class name is added to the dict for deserialisation purpose
        SliceModel_dict["__class__"] = "SliceModel"
        return SliceModel_dict

    def _set_None(self):
        """Set all the properties to None (except pyleecan object)"""

        self.type_distribution = None
        self.Nslices = None
        self.z_list = None
        self.angle_rotor = None
        self.angle_stator = None
        self.L = None
        self.is_step = None
        self.is_skew = None

    def _get_type_distribution(self):
        """getter of type_distribution"""
        return self._type_distribution

    def _set_type_distribution(self, value):
        """setter of type_distribution"""
        check_var("type_distribution", value, "str")
        self._type_distribution = value

    type_distribution = property(
        fget=_get_type_distribution,
        fset=_set_type_distribution,
        doc=u"""Type of slice distribution to use for rotor skew if linear and continuous ("uniform", "gauss", "user-defined")

        :Type: str
        """,
    )

    def _get_Nslices(self):
        """getter of Nslices"""
        return self._Nslices

    def _set_Nslices(self, value):
        """setter of Nslices"""
        check_var("Nslices", value, "int")
        self._Nslices = value

    Nslices = property(
        fget=_get_Nslices,
        fset=_set_Nslices,
        doc=u"""Number of slices

        :Type: int
        """,
    )

    def _get_z_list(self):
        """getter of z_list"""
        return self._z_list

    def _set_z_list(self, value):
        """setter of z_list"""
        if type(value) is int and value == -1:
            value = list()
        check_var("z_list", value, "list")
        self._z_list = value

    z_list = property(
        fget=_get_z_list,
        fset=_set_z_list,
        doc=u"""List of slice positions

        :Type: list
        """,
    )

    def _get_angle_rotor(self):
        """getter of angle_rotor"""
        return self._angle_rotor

    def _set_angle_rotor(self, value):
        """setter of angle_rotor"""
        if type(value) is int and value == -1:
            value = array([])
        elif type(value) is list:
            try:
                value = array(value)
            except:
                pass
        check_var("angle_rotor", value, "ndarray")
        self._angle_rotor = value

    angle_rotor = property(
        fget=_get_angle_rotor,
        fset=_set_angle_rotor,
        doc=u"""Array of rotor skew angles in case of skew

        :Type: ndarray
        """,
    )

    def _get_angle_stator(self):
        """getter of angle_stator"""
        return self._angle_stator

    def _set_angle_stator(self, value):
        """setter of angle_stator"""
        if type(value) is int and value == -1:
            value = array([])
        elif type(value) is list:
            try:
                value = array(value)
            except:
                pass
        check_var("angle_stator", value, "ndarray")
        self._angle_stator = value

    angle_stator = property(
        fget=_get_angle_stator,
        fset=_set_angle_stator,
        doc=u"""Array of stator skew angles in case of skew

        :Type: ndarray
        """,
    )

    def _get_L(self):
        """getter of L"""
        return self._L

    def _set_L(self, value):
        """setter of L"""
        check_var("L", value, "float")
        self._L = value

    L = property(
        fget=_get_L,
        fset=_set_L,
        doc=u"""Machine length (mean of rotor/stator lengths)

        :Type: float
        """,
    )

    def _get_is_step(self):
        """getter of is_step"""
        return self._is_step

    def _set_is_step(self, value):
        """setter of is_step"""
        check_var("is_step", value, "bool")
        self._is_step = value

    is_step = property(
        fget=_get_is_step,
        fset=_set_is_step,
        doc=u"""True to define slices as steps

        :Type: bool
        """,
    )

    def _get_is_skew(self):
        """getter of is_skew"""
        return self._is_skew

    def _set_is_skew(self, value):
        """setter of is_skew"""
        check_var("is_skew", value, "bool")
        self._is_skew = value

    is_skew = property(
        fget=_get_is_skew,
        fset=_set_is_skew,
        doc=u"""True if slices account for skew

        :Type: bool
        """,
    )

reason: from numpy


text:
import pandas as pd
import numpy as np
import xgboost as xgb

df = pd.read_csv("dset2.csv")

dtrain = xgb.DMatrix(df[["x1","x2","x3","x4"]], label=df["y"])

def gamma_grad(preds, acts):
 return (acts-preds)/preds

def gamma_hess(preds, acts):
 return -acts/preds

def gamma_obj(preds, dtrain):
 acts = dtrain.get_label()
 preds = np.exp(preds)
 grad = -gamma_grad(preds, acts)
 hess = -gamma_hess(preds, acts)
 return grad,hess

bst = xgb.train({"max_depth":1}, dtrain, 1000, obj = gamma_obj)

preds = np.exp(np.array(bst.predict(dtrain)))
acts = np.array(df["y"])

print(preds[:12])
print(acts[:12])

Lump = np.array([0]*len(df))

def gamma_grad_around_Lump(preds, acts): #Haskell forgive me for this parody of functional programming
 return (preds/((preds+Lump)**2))*(acts-(preds+Lump))

def gamma_hess_around_Lump(preds, acts):
 return (preds/((preds+Lump)**3))*(acts*Lump-Lump*preds-Lump*Lump-preds*acts)

def gamma_obj_around_Lump(preds, dtrain):
 acts = dtrain.get_label()
 preds = np.exp(preds)
 grad = -gamma_grad_around_Lump(preds,acts)
 hess = -gamma_hess_around_Lump(preds,acts)
 return grad, hess

bst = xgb.train({"max_depth":1}, dtrain, 1000, obj = gamma_obj_around_Lump)

preds = np.exp(np.array(bst.predict(dtrain)))
acts = np.array(df["y"])

print(preds[:12])
print(acts[:12])

BV = np.log(0.5*sum(df["y"])/len(df))

bstA = xgb.train({"max_depth":1, "base_score":BV}, dtrain, 0, obj = gamma_obj_around_Lump)
bstB = xgb.train({"max_depth":1, "base_score":BV}, dtrain, 0, obj = gamma_obj_around_Lump)

for i in range(30):
 Lump = np.array(np.exp(bstB.predict(dtrain)))
 bstA = xgb.train({"max_depth":1, "base_score":BV}, dtrain, i, obj = gamma_obj_around_Lump, xgb_model = bstA)
 Lump = np.array(np.exp(bstA.predict(dtrain)))
 bstB = xgb.train({"max_depth":1, "base_score":BV}, dtrain, i, obj = gamma_obj_around_Lump, xgb_model = bstB)

predsA = np.array(np.exp(bstA.predict(dtrain)))
predsB = np.array(np.exp(bstB.predict(dtrain)))
preds = predsA + predsB
acts = np.array(df["y"])

print(predsA[:12])
print(predsB[:12])
print(preds[:12])
print(acts[:12])

print(df[:12])

max_stars_repo_name: H-B-P/parallelized-gamma-xgboost

max_issues_repo_name: H-B-P/parallelized-gamma-xgboost

max_forks_repo_name: H-B-P/parallelized-gamma-xgboost


text:
import numpy as np
from matplotlib import pyplot as plt
from sklearn.model_selection import KFold
from sklearn.tree import DecisionTreeRegressor

def get_features_targets(data):
  features = np.zeros((data.shape[0], 4)) #n lines, 4 columns
  features[:,0] = data['u'] - data['g']
  features[:,1] = data['g'] - data['r']
  features[:,2] = data['r'] - data['i']
  features[:,3] = data['i'] - data['z']
  targets = data['redshift']
  return (features, targets)

def median_diff(predicted, actual):
  diff = np.median(np.absolute(predicted - actual))
  return diff

def cross_validate_predictions(model, features, targets, k):
  kf = KFold(n_splits=k, shuffle=True)

  # declare an array for predicted redshifts from each iteration
  all_predictions = np.zeros_like(targets)

  for train_indices, test_indices in kf.split(features):
    print (train_indices, len(train_indices), test_indices, len(test_indices))
    # split the data into training and testing
    train_features, test_features = features[train_indices], features[test_indices]
    train_targets, test_targets = targets[train_indices], targets[test_indices]

    # fit the model for the current set
    model.fit(train_features, train_targets)

    # predict using the model
    predictions = model.predict(test_features)

    # put the predicted values in the all_predictions array defined above
    all_predictions[test_indices] = predictions

    # note: the line above is filling multiple disjoint elements of the
    #array at the same time. for example, if:
    #test_indexes is 8 11 23
    #with test_features z8, z11, z23
    #the code does: allpredictions[8,11,23] = z8, z11, z23
    #and hence, once per galaxy

  # return the predictions
  return all_predictions

if __name__ == "__main__":
  data = np.load('./sdss_galaxy_colors.npy')
  features, targets = get_features_targets(data)

  # initialize model
  dtr = DecisionTreeRegressor(max_depth=19)

  # call your cross validation function
  predictions = cross_validate_predictions(dtr, features, targets, 10)

  # calculate and print the rmsd as a sanity check
  diffs = median_diff(predictions, targets)
  print('Median difference: {:.3f}'.format(diffs))

  # plot the results to see how well our model looks
  plt.scatter(targets, predictions, s=0.4)
  plt.xlim((0, targets.max()))
  plt.ylim((0, predictions.max()))
  plt.xlabel('Measured Redshift')
  plt.ylabel('Predicted Redshift')
  plt.show()

max_stars_repo_name: lokijota/datadrivenastronomymooc

max_issues_repo_name: lokijota/datadrivenastronomymooc

max_forks_repo_name: lokijota/datadrivenastronomymooc


text:
from itertools import combinations, product

import numpy as np
from pyDOE import lhs
from scipy.special import comb

def normalize(vectors):
    """
    Normalize a set of vectors.

    The length of the returned vectors will be unity.

    Parameters
    ----------
    vectors : np.ndarray
        Set of vectors of any length, except zero.

    """
    if len(np.asarray(vectors).shape) == 1:
        return vectors / np.linalg.norm(vectors)
    norm = np.linalg.norm(vectors, axis=1)
    return vectors / norm[:, np.newaxis]

def shear(vectors, degrees: float = 5):
    """
    Shear a set of vectors lying on the plane z=0 towards the z-axis, such that the
    resulting vectors 'degrees' angle away from the z axis.

    z is the last element of the vector, and has to be equal to zero.

    Parameters
    ----------
    vectors : numpy.ndarray
        The final element of each vector should be zero.
    degrees : float, optional
        The angle that the resultant vectors make with the z axis. Unit is radians.
        (the default is 5)

    """
    angle = degrees * np.pi / 180
    m = 1 / np.tan(angle)
    norm = np.linalg.norm(vectors, axis=1)
    vectors[:, -1] += norm * m
    return normalize(vectors)

def rotate(initial_vector, rotated_vector, other_vectors):
    """Calculate the rotation matrix that rotates the initial_vector to the
    rotated_vector. Apply that rotation on other_vectors and return.
    Uses Householder reflections twice to achieve this."""

    init_vec_norm = normalize(initial_vector)
    rot_vec_norm = normalize(np.asarray(rotated_vector))
    middle_vec_norm = normalize(init_vec_norm + rot_vec_norm)
    first_reflector = init_vec_norm - middle_vec_norm
    second_reflector = middle_vec_norm - rot_vec_norm
    Q1 = householder(first_reflector)
    Q2 = householder(second_reflector)
    reflection_matrix = np.matmul(Q2, Q1)
    rotated_vectors = np.matmul(other_vectors, np.transpose(reflection_matrix))
    return rotated_vectors

def householder(vector):
    """Return reflection matrix via householder transformation."""
    identity_mat = np.eye(len(vector))
    v = vector[np.newaxis]
    denominator = np.matmul(v, v.T)
    numerator = np.matmul(v.T, v)
    rot_mat = identity_mat - (2 * numerator / denominator)
    return rot_mat

def rotate_toward(initial_vector, final_vector, other_vectors, degrees: float = 5):
    """
    Rotate other_vectors (with the centre at initial_vector) towards final_vector
    by an angle degrees.

    Parameters
    ----------
    initial_vector : np.ndarray
        Centre of the vectors to be rotated.
    final_vector : np.ndarray
        The final position of the center of other_vectors.
    other_vectors : np.ndarray
        The array of vectors to be rotated
    degrees : float, optional
        The amount of rotation (the default is 5)

    Returns
    -------
    rotated_vectors : np.ndarray
        The rotated vectors
    reached: bool
        True if final_vector has been reached
    """
    final_vector = normalize(final_vector)
    initial_vector = normalize(initial_vector)
    cos_phi = np.dot(initial_vector, final_vector)
    theta = degrees * np.pi / 180
    cos_theta = np.cos(theta)
    phi = np.arccos(cos_phi)
    if phi < theta:
        return (rotate(initial_vector, final_vector, other_vectors), True)
    cos_phi_theta = np.cos(phi - theta)
    A = np.asarray([[cos_phi, 1], [1, cos_phi]])
    B = np.asarray([cos_phi_theta, cos_theta])
    x = np.linalg.solve(A, B)
    rotated_vector = x[0] * initial_vector + x[1] * final_vector
    return (rotate(initial_vector, rotated_vector, other_vectors), False)

class ReferenceVectors:
    """Class object for reference vectors."""

    def __init__(
        self,
        lattice_resolution: int = None,
        number_of_objectives: int = None,
        creation_type: str = "Uniform",
        vector_type: str = "Spherical",
        ref_point: list = None,
    ):
        """Create a Reference vectors object.

        A simplex lattice is formed

        Parameters
        ----------
        lattice_resolution : int
            Number of divisions along an axis when creating the simplex lattice.
        number_of_objectives : int
            Number of objectives.
        creation_type : str, optional
            'Uniform' creates the reference vectors uniformly using simplex lattice
            design. 'Focused' creates reference vectors symmetrically around a central
            reference vector. 'Reversed' coming soon.By default 'Uniform'.
        vector_type : str, optional
            'Spherical' normalizes the vectors to a hypersphere, i.e. the second norm
            is equal to 1. 'Planar' normalizes vectors to a plane, i.e. the first norm
            is equal to 1. By default 'Spherical'.
        ref_point : list, optional
            User preference information for a priori methods.
        """

        self.number_of_objectives = number_of_objectives
        self.lattice_resolution = lattice_resolution
        self.number_of_vectors = 0
        self.creation_type = creation_type
        self.vector_type = vector_type
        self.values = []
        self.values_planar = []
        self.ref_point = [1] * number_of_objectives if ref_point is None else ref_point
        self._create(creation_type)
        self.initial_values = np.copy(self.values)
        self.initial_values_planar = np.copy(self.values_planar)
        self.neighbouring_angles()
        # self.iteractive_adapt_1() Can use this for a priori preferences!

    def _create(self, creation_type: str = "Uniform"):
        """Create the reference vectors.

        Parameters
        ----------
        creation_type : str, optional
            'Uniform' creates the reference vectors uniformly using simplex lattice
            design. 'Focused' creates reference vectors symmetrically around a central
            reference vector. By default 'Uniform'.
        """
        if creation_type == "Uniform":
            number_of_vectors = comb(
                self.lattice_resolution + self.number_of_objectives - 1,
                self.number_of_objectives - 1,
                exact=True,
            )
            self.number_of_vectors = number_of_vectors
            temp1 = range(1, self.number_of_objectives + self.lattice_resolution)
            temp1 = np.array(list(combinations(temp1, self.number_of_objectives - 1)))
            temp2 = np.array(
                [range(self.number_of_objectives - 1)] * self.number_of_vectors
            )
            temp = temp1 - temp2 - 1
            weight = np.zeros(
                (self.number_of_vectors, self.number_of_objectives), dtype=int
            )
            weight[:, 0] = temp[:, 0]
            for i in range(1, self.number_of_objectives - 1):
                weight[:, i] = temp[:, i] - temp[:, i - 1]
            weight[:, -1] = self.lattice_resolution - temp[:, -1]
            self.values = weight / self.lattice_resolution
            self.values_planar = np.copy(self.values)
            self.normalize()
            return
        elif creation_type == "Focused":
            point_set = [[0, 1, -1]] * (self.number_of_objectives - 1)
            # The cartesian product of point_set.
            initial = np.array(list(product(*point_set)))[1:]
            # First element was removed because of the error during normalization.
            initial = normalize(initial)
            initial = np.hstack((initial, np.zeros((initial.shape[0], 1))))
            final = shear(initial, degrees=5)
            # Adding the first element back
            final = np.vstack(([0] * (self.number_of_objectives - 1) + [1], final))
            self.number_of_vectors = final.shape[0]
            self.values = rotate(final[0], self.ref_point, final)
            self.values_planar = np.copy(self.values)
            self.normalize()
            self.add_edge_vectors()
        elif creation_type == "Sparse_Focused":
            initial = np.eye(self.number_of_objectives - 1)
            initial = np.vstack((initial, -initial))
            initial = normalize(initial)
            initial = np.hstack((initial, np.zeros((initial.shape[0], 1))))
            final = shear(initial, degrees=5)
            # Adding the first element back
            final = np.vstack(([0] * (self.number_of_objectives - 1) + [1], final))
            self.number_of_vectors = final.shape[0]
            self.values = rotate(final[0], self.ref_point, final)
            self.values_planar = np.copy(self.values)
            self.normalize()
            self.add_edge_vectors()

    def normalize(self):
        """Normalize the reference vectors to a unit hypersphere."""
        self.number_of_vectors = self.values.shape[0]
        norm_2 = np.linalg.norm(self.values, axis=1)
        norm_1 = np.sum(self.values_planar, axis=1)
        norm_2 = np.repeat(norm_2, self.number_of_objectives).reshape(
            self.number_of_vectors, self.number_of_objectives
        )
        norm_1 = np.repeat(norm_1, self.number_of_objectives).reshape(
            self.number_of_vectors, self.number_of_objectives
        )
        norm_2[norm_2 == 0] = np.finfo(float).eps
        self.values = np.divide(self.values, norm_2)
        self.values_planar = np.divide(self.values_planar, norm_1)

    def neighbouring_angles(self) -> np.ndarray:
        """Calculate neighbouring angles for normalization."""
        cosvv = np.dot(self.values, self.values.transpose())
        cosvv.sort(axis=1)
        cosvv = np.flip(cosvv, 1)
        cosvv[cosvv > 1] = 1
        acosvv = np.arccos(cosvv[:, 1])
        self.neighbouring_angles_current = acosvv
        return acosvv

    def adapt(self, fitness: np.ndarray):
        """Adapt reference vectors. Then normalize.

        Parameters
        ----------
        fitness : np.ndarray
        """
        max_val = np.amax(fitness, axis=0)
        min_val = np.amin(fitness, axis=0)
        self.values = np.multiply(
            self.initial_values,
            np.tile(np.subtract(max_val, min_val), (self.number_of_vectors, 1)),
        )
        self.normalize()

    def interactive_adapt_1(self, z: np.ndarray, n_solutions: int, translation_param: float = 0.2) -> None:
        """
        Adapt reference vectors using the information about prefererred solution(s) selected by the Decision maker.

        Args:
            z (np.ndarray): Preferred solution(s).
            n_solutions (int): Number of solutions in total.
            translation_param (float): Parameter determining how close the reference vectors are to the central vector
            **v** defined by using the selected solution(s) z.

        Returns:

        """

        if z.shape[0] == n_solutions:
            # if dm specifies all solutions as preferred, reinitialize reference vectors
            self.values = self.initial_values
            self.values_planar = self.initial_values_planar

        else:
            if z.shape[0] == 1:
                # single preferred solution
                # calculate new reference vectors
                self.values = translation_param * self.initial_values + ((1 - translation_param) * z)
                self.values_planar = translation_param * self.initial_values_planar + ((1 - translation_param) * z)

            else:
                # multiple preferred solutions
                # calculate new reference vectors for each preferred solution
                values = [translation_param * self.initial_values + ((1 - translation_param) * z_i) for z_i in z]
                values_planar = [translation_param * self.initial_values_planar + ((1 - translation_param) * z_i)
                                 for z_i in z]

                # combine arrays of reference vectors into a single array and update reference vectors
                self.values = np.concatenate(values)
                self.values_planar = np.concatenate(values_planar)

        self.normalize()

    def interactive_adapt_2(self, z: np.ndarray, n_solutions: int, predefined_distance: float = 0.2) -> None:
        """
        Adapt reference vectors by using the information about non-preferred solution(s) selected by the Decision maker.
        After the Decision maker has specified non-preferred solution(s), Euclidian distance between normalized solution
        vector(s) and each of the reference vectors are calculated. Those reference vectors that are **closer** than a
        predefined distance are either **removed** or **re-positioned** somewhere else.

        Note:
            At the moment, only the **removal** of reference vectors is supported. Repositioning of the reference
            vectors is **not** supported.

        Note:
            In case the Decision maker specifies multiple non-preferred solutions, the reference vector(s) for which the
            distance to **any** of the non-preferred solutions is less than predefined distance are removed.

        Note:
            Future developer should implement a way for a user to say: "Remove some percentage of
            objecive space/reference vectors" rather than giving a predefined distance value.

        Args:
            z (np.ndarray): Non-preferred solution(s).
            n_solutions (int): Number of solutions in total.
            predefined_distance (float): The reference vectors that are closer than this distance are either removed or
            re-positioned somewhere else.
            Default value: 0.2

        Returns:

        """

        if z.shape[0] == n_solutions:
            # if dm specifies all solutions as non-preferred ones, reinitialize reference vectors
            self.values = self.initial_values
            self.values_planar = self.initial_values_planar
            self.normalize()

        else:
            # calculate L1 norm of non-preferred solution(s)
            z = np.atleast_2d(z)
            norm = np.linalg.norm(z, ord=1, axis=1).reshape(np.shape(z)[0], 1)

            # non-preferred solutions normalized
            v_c = np.divide(z, norm)

            # distances from non-preferred solution(s) to each reference vector
            distances = np.array([list(map(lambda solution: np.linalg.norm(solution - value, ord=2), v_c))
                                  for value in self.values_planar])

            # find out reference vectors that are not closer than threshold value to any non-preferred solution
            mask = [all(d >= predefined_distance) for d in distances]

            # set those reference vectors that met previous condition as new reference vectors, drop others
            self.values = self.values[mask]
            self.values_planar = self.values_planar[mask]

    def iteractive_adapt_3(self, ref_point, translation_param=0.2):
        """Adapt reference vectors linearly towards a reference point. Then normalize.

        The details can be found in the following paper: Hakanen, Jussi &
        Chugh, Tinkle & Sindhya, Karthik & Jin, Yaochu & Miettinen, Kaisa.
        (2016). Connections of Reference Vectors and Different Types of
        Preference Information in Interactive Multiobjective Evolutionary
        Algorithms.

        Parameters
        ----------
        ref_point :

        translation_param :
            (Default value = 0.2)

        """
        self.values = self.initial_values * translation_param + (
            (1 - translation_param) * ref_point
        )
        self.values_planar = self.initial_values_planar * translation_param + (
            (1 - translation_param) * ref_point
        )
        self.normalize()

    def interactive_adapt_4(self, preferred_ranges: np.ndarray) -> None:
        """
        Adapt reference vectors by using the information about the Decision maker's preferred range for each of the
        objective. Using these ranges, Latin hypercube sampling is applied to generate m number of samples between
        within these ranges, where m is the number of reference vectors. Normalized vectors constructed of these samples
        are then set as new reference vectors.

        Args:
            preferred_ranges (np.ndarray): Preferred lower and upper bound for each of the objective function values.

        Returns:

        """

        # bounds
        lower_limits = np.array([ranges[0] for ranges in preferred_ranges])
        upper_limits = np.array([ranges[1] for ranges in preferred_ranges])

        # generate samples using Latin hypercube sampling
        w = lhs(self.number_of_objectives, samples=self.number_of_vectors)

        # scale between bounds
        w = w * (upper_limits - lower_limits) + lower_limits

        # set new reference vectors and normalize them
        self.values = w
        self.values_planar = w
        self.normalize()

    def slow_interactive_adapt(self, ref_point):
        """Basically a wrapper around rotate_toward. Slowly rotate ref vectors toward
        ref_point. Return a boolean value to tell if the ref_point has been reached.

        Parameters
        ----------
        ref_point : list or np.ndarray
            The reference vectors will slowly move towards the ref_point.

        Returns
        -------
        boolean
            True if ref_point has been reached. False otherwise.
        """
        assert self.creation_type == "Focused" or self.creation_type == "Sparse_Focused"
        if np.array_equal(self.values[0], ref_point):
            return
        self.values, reached = rotate_toward(
            self.values[0], ref_point, self.values[0 : -self.number_of_objectives]
        )
        self.values_planar = self.values
        self.add_edge_vectors()
        self.normalize()
        return reached

    def add_edge_vectors(self):
        """Add edge vectors to the list of reference vectors.

        Used to cover the entire orthant when preference information is
        provided.

        """
        edge_vectors = np.eye(self.values.shape[1])
        self.values = np.vstack([self.values, edge_vectors])
        self.values_planar = np.vstack([self.values_planar, edge_vectors])
        self.number_of_vectors = self.values.shape[0]
        self.normalize()

max_stars_repo_name: gialmisi/desdeo-emo

max_issues_repo_name: gialmisi/desdeo-emo

max_forks_repo_name: gialmisi/desdeo-emo


text:
# -*- coding: utf-8 -*-
import numpy as np
import json
import re

"""
1 = train
2 = test
3 = dev
"""
TEST = "./data/test.json"
DEV = "./data/dev.json"
TRAIN = "./data/train.json"
DATA_FILE = [None, TRAIN, TEST, DEV]
DATA_JSON = [None, [], [], []]
LABEL = 5

def clean_str(string):
    string = re.sub(r"[^A-Za-z0-9(),.!?\-\'\`]", " ", string)
    string = re.sub(r"\'s", " \'s", string)
    string = re.sub(r"\'ve", " \'ve", string)
    string = re.sub(r"n\'t", " n\'t", string)
    string = re.sub(r"\'re", " \'re", string)
    string = re.sub(r"\'d", " \'d", string)
    string = re.sub(r"\'ll", " \'ll", string)
    string = re.sub(r",", " , ", string)
    string = re.sub(r"!", " ! ", string)
    string = re.sub(r"\(", "-LRB-", string)
    string = re.sub(r"\)", "-RRB-", string)
    string = re.sub(r"\?", " \? ", string)
    string = re.sub(r"\s{2,}", " ", string)
    return string.strip().lower()

def gat_label(sent):
    # [0, 0.2], (0.2, 0.4], (0.4, 0.6], (0.6, 0.8], (0.8, 1.0]
    if sent >= 0 and sent <= .2:
        return 0
    elif sent <= .4:
        return 1
    elif sent <= .6:
        return 2
    elif sent <= .8:
        return 3
    elif sent <= 1.0:
        return 4
    else:
        return 2

def load_data_and_labels(data_file):
    # Load data from json file
    with open(data_file, "r", encoding='utf8') as f:
        items = json.load(f)
    text = [item[0] for item in items]
    label = []
    for item in items:
        la = [0] * LABEL
        la[item[1]] = 1
        label.append(la)
    # label = [item[1] for item in items]
    return (text, label)

def batch_iter(data, batch_size, num_epochs, shuffle=True):
    data = np.array(data)
    data_size = len(data)
    num_batches_per_epoch = int((data_size-1)/batch_size) + 1
    for epoch in range(num_epochs):
        # Shuffle the data at each epoch
        if shuffle:
            shuffle_indices = np.random.permutation(np.arange(data_size))
            shuffled_data = data[shuffle_indices]
        else:
            shuffled_data = data
        for batch_num in range(num_batches_per_epoch):
            start_index = batch_num * batch_size
            end_index = min((batch_num + 1) * batch_size, data_size)
            print ("start: %d, end: %d" % (start_index, end_index))
            yield shuffled_data[start_index:end_index]

def main():
    dictionary = {}
    with open("./stanfordSentimentTreebank/dictionary.txt", "r", encoding="utf8") as f:
        for line in f.readlines():
            key, value = line.strip().split("|")
            dictionary[clean_str(key)] = value

    split = {}
    with open("./stanfordSentimentTreebank/datasetSplit.txt", "r", encoding="utf8") as f:
        for i, line in enumerate(f.readlines()):
            if i > 0:
                key, value = line.strip().split(",")
                split[key] = value

    sentiment = {}
    with open("./stanfordSentimentTreebank/sentiment_labels.txt", "r", encoding="utf8") as f:
        for i, line in enumerate(f.readlines()):
            if i > 0:
                key, value = line.strip().split("|")
                sentiment[key] = value

    with open("./stanfordSentimentTreebank/datasetSentences.txt", "r", encoding="utf8") as f:
        for i, line in enumerate(f.readlines()):
            if i > 0:
                index, sentence = line.strip().split("\t")
                # print (": ".join((index, sentence)))
                sentence = clean_str(sentence)
                DATA_JSON[int(split[index])].append((sentence, gat_label(float(sentiment[dictionary[sentence]]))))

    for i in range(1, 4):
        print ("json.dump %s" % i)
        with open(DATA_FILE[i], "w", encoding="utf8") as f:
            json.dump(DATA_JSON[i], f)

if __name__ == '__main__':
    main()

max_stars_repo_name: DesperadoAdil/SentimentAnalysis

max_issues_repo_name: DesperadoAdil/SentimentAnalysis

max_forks_repo_name: DesperadoAdil/SentimentAnalysis


text:
#-*- coding:utf-8 -*-

from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import os
import sys
import pickle
import numpy as np

from keras.models import Model
from keras.preprocessing import image
from keras.layers import Flatten, Input
from keras.applications.resnet50 import ResNet50
from keras.applications.resnet50 import preprocess_input, decode_predictions
from keras.utils.vis_utils import plot_model
import gc

res50_base_model = ResNet50(weights='imagenet', pooling=max, include_top=False)
image_input = Input(shape=(224, 224, 3), name="image_input")
x = res50_base_model(image_input)
x = Flatten()(x)
res50_model = Model(inputs=image_input, outputs=x)
res50_model.summary()
# plot_model(res50_base_model, to_file='model.png', show_shapes=True)

batch_size = 1000
# image file abs_path, label
image_dir = "/home/ai-i-liuguiyang/ImageRetireval/dataset/OxBuild/src/"
image_list_file_path = "/home/ai-i-liuguiyang/ImageRetireval/dataset/OxBuild/src/index_file.csv"
nn_feature_save_dir = "/home/ai-i-liuguiyang/ImageRetireval/dataset/OxBuild/src/"

with open(image_list_file_path, "r") as fl_reader:

    image_nn_feature_dict = dict()

    def __fetch_nn_feature(batch_image, batch_file_name):
        batch_image = np.concatenate(batch_image, axis=0)
        x = preprocess_input(batch_image)

        features = res50_model.predict(x)
        features_reduce = features.squeeze()
        for idx in range(len(batch_file_name)):
            image_nn_feature_dict[batch_file_name[idx]] = features_reduce[idx]
        # print(features_reduce)
        print(features_reduce.shape)

    batch_image, batch_file_name = [], []
    for line in fl_reader.readlines():
        image_name, image_label = line.strip().split(",")
        image_path = image_dir + image_name
        if not os.path.exists(image_path):
            print("{} not found !".format(image_path))
            continue
        image_data = image.load_img(image_path, target_size=(224, 224))
        x = image.img_to_array(image_data)
        x = np.expand_dims(x, axis=0)
        batch_image.append(x)
        batch_file_name.append(image_name)
        if len(batch_image) == batch_size:
            __fetch_nn_feature(batch_image, batch_file_name)
            batch_image = list()
            batch_file_name = list()

    if len(batch_image):
        __fetch_nn_feature(batch_image, batch_file_name)

    print("Before dump the data !")
    pickle.dump(image_nn_feature_dict, open(nn_feature_save_dir+"nn_features.pkl", "wb"), True)
    print("After dump the data !")

gc.collect()

max_stars_repo_name: liuguiyangnwpu/MassImageRetrieval


text:
import pandas as pd
import numpy as np
from collections import OrderedDict

def calculate_profits(
        bids,
        transactions,
        reservation_prices=None,
        fees=None,
        **kwargs):
    """Extras from the transactions and the bids the profit
    of each player and the market maker

    Parameters
    ----------
    bids : pd.DataFrame
        Collections of bids to be used

    transactions : pd.DataFrame
        Collection of transactions to be taken into account
    reservation_prices : dict, (Default value = None)
        Maping between users and their reservation prices. If None,
        it is assumed that each user bided truthfully and the
        information is extracted from the bid.
    fees : np.ndarray, (Default value = None)
        List of fees that each user has to pay to trade in the market.

    Returns
    -------
    profit: dict
        A dictionary with three values:
        * player_bid: A list with the profits of each user
        using their bids as reservation prices
        * player_reservation: Same as above but using
        their reservation prices, if none are provided,
        it is the same as `player_bid`
        * market: profit of the market maker

    Examples
    ---------
    >>> tm = pm.TransactionManager()
    >>> bm = pm.BidManager()
    >>> bm.add_bid(1, 3, 0)
    0
    >>> bm.add_bid(1, 2, 1)
    1
    >>> bm.add_bid(1.5, 1, 2, False)
    2
    >>> tm.add_transaction(0, 1, 2, 2, False)
    0
    >>> tm.add_transaction(2, 1, 2, 0, False)
    1
    >>> rp = {2: 0}
    >>> profits = calculate_profits(bm.get_df(), tm.get_df(),
    ...        reservation_prices=rp)
    >>> profits['player_bid']
    array([1., 0., 1.])
    >>> profits['player_reservation']
    array([1., 0., 2.])
    >>> profits['market']
    0.0
    """
    users = sorted(bids.user.unique())
    buyers = bids.loc[bids['buying']].index.values
    sellers = bids.loc[~bids['buying']].index.values

    if reservation_prices is None:
        reservation_prices = OrderedDict()
    for i, x in bids.iterrows():
        if i not in reservation_prices:
            reservation_prices[i] = x.price

    if fees is None:
        fees = np.zeros(bids.user.unique().shape[0])

    profit = OrderedDict()
    for case in ['bid', 'reservation']:
        tmp = bids.reset_index().rename(columns={'index': 'bid'}).copy()
        tmp = tmp[['bid', 'price', 'buying', 'user']]
        if case == 'reservation':
            tmp.price = tmp.apply(
                lambda x: reservation_prices.get(x.bid, x.price), axis=1)
        merged = transactions.merge(tmp, on='bid').copy()
        merged['gain'] = merged.apply(lambda x: get_gain(x), axis=1)
        profit_player = merged.groupby('user')['gain'].sum()
        # print(profit_player)
        profit_player = np.array([profit_player.get(x, 0) for x in users])
        profit['player_{}'.format(case)] = profit_player.astype('float64')

        if case == 'bid':
            # print(merged)
            mb = merged.loc[merged['buying']]
            ms = merged.loc[~merged['buying']]
            # print(ms)
            # print(ms.quantity.sum(), mb.quantity.sum())
            # print(ms.price_x * ms.quantity)
            profit_market = (mb.price_x * mb.quantity).values.sum()
            profit_market -= (ms.price_x * ms.quantity).values.sum()
            profit_market += fees.sum()
            profit['market'] = profit_market.astype('float64')

    return profit

def get_gain(row):
    """Finds the gain of the row

    Parameters
    ----------
    row : pandas row
       Row obtained by merging a transaction with a
       bid dataframe

    Returns
    -------
    gain
        The gain obtained by the row

    """
    gap = row.price_y - row.price_x
    if not row.buying:
        gap = - gap
    return gap * row.quantity

max_stars_repo_name: kiedanski/pymarket

max_issues_repo_name: kiedanski/pymarket


text:
# -*- coding: utf-8 -*-

from ....Classes.OutStruct import OutStruct
from ....Methods.Simulation.Input import InputError
from numpy import ndarray

def gen_input(self):
    """Generate the input for the structural module (skip force computation)

    Parameters
    ----------
    self : InputForce
        An InputForce object
    """

    output = OutStruct()
    # Load and check time
    if self.time is None:
        raise InputError("ERROR: InputForce.time missing")
    output.time = self.time.get_data()

    if not isinstance(output.time, ndarray) or len(output.time.shape) != 1:
        # time should be a vector
        raise InputError(
            "ERROR: InputForce.time should be a vector, "
            + str(output.time.shape)
            + " shape found"
        )
    Nt_tot = len(output.time)

    # Load and check angle
    if self.angle is None:
        raise InputError("ERROR: InputForce.angle missing")
    output.angle = self.angle.get_data()
    if not isinstance(output.angle, ndarray) or len(output.angle.shape) != 1:
        # angle should be a vector
        raise InputError(
            "ERROR: InputForce.angle should be a vector, "
            + str(output.angle.shape)
            + " shape found"
        )
    Na_tot = len(output.angle)

    if self.Prad is None:
        raise InputError("ERROR: InputForce.Prad missing")
    output.Prad = self.Prad.get_data()
    if not isinstance(output.Prad, ndarray) or output.Prad.shape != (Nt_tot, Na_tot):
        raise InputError(
            "ERROR: InputForce.Prad must be a matrix with the shape "
            + str((Nt_tot, Na_tot))
            + " (len(time), stator phase number), "
            + str(output.Prad.shape)
            + " returned"
        )

    if self.Ptan is not None:
        output.Ptan = self.Ptan.get_data()
        if not isinstance(output.Ptan, ndarray) or output.Ptan.shape != (
            Nt_tot,
            Na_tot,
        ):
            raise InputError(
                "ERROR: InputForce.Ptan must be a matrix with the shape "
                + str((Nt_tot, Na_tot))
                + " (len(time), rotor phase number), "
                + str(output.Ptan.shape)
                + " returned"
            )
    else:
        output.Ptan = None

    if self.parent.parent is None:
        raise InputError(
            "ERROR: The Simulation object must be in an Output object to run"
        )
    # Save the Output in the correct place
    self.parent.parent.struct = output

max_stars_repo_name: helene-t/pyleecan

max_issues_repo_name: magnetron/pyleecan

max_forks_repo_name: magnetron/pyleecan