text:
"""Plot tools: Reusable utilities for working with matplotlib plots.

Contents:
params_(reports|poster|...)
                    Standard parameter sets for final-quality plots
colors              Easily distinguishable color sets

Functions:
update_color_cycle  Use the distinguishable colors by default
plot_smoothed       Plot a time series smoothed with a cosine kernel
thin_points         Thin data points that are too close to each other
scatter_thin_points Plot a set of thinned data points with weight info
scatter_mark_outliers
                    Mark outliers at the edges of a scatterplot
scatter_outliers_size
                    Make a scatterplot with sizes and outliers

"""

import matplotlib as mpl
from matplotlib import pyplot
import numpy as np
from scipy import signal
from scipy.spatial import KDTree

# Matplotlib params for report-quality plots
params_orig = mpl.rcParams.copy()
params_reports = {
    'font.family': 'TeX Gyre Schola',
    'font.size': 10,
    'legend.fontsize': 'small',
    'axes.labelsize': 'small',
    #'font.serif': ['New Century Schoolbook', 'CMU Serif'],
    'font.sans-serif': [],
    'mathtext.default': 'regular',
    'pgf.texsystem': 'pdflatex',
    'pgf.preamble': ['\\usepackage{siunitx}'],
    'figure.figsize': (6.0, 3.7)
    #'savefig.dpi': 150,
}
params_poster_a0 = dict(params_reports)
params_poster_a0.update({
    'font.size': 25,
    'figure.figsize': (12.0, 8.0),
    'font.family': 'serif',
    'font.serif': ['Linux Libertine', 'CMU Serif'],
    'font.sans-serif': [],
    'legend.fontsize': 'large',
    'axes.labelsize': 'large',
})

# Some colours for easily distinguishable and colourblind-friendly coding
# Use the 4th colour if necessary; may not be colourblind-friendly though
# from colorbrewer.org
colors = {}
colors['green'] = '#1b9e77'
colors['orange'] = '#d95f02'
colors['purple'] = '#7570b3'
colors['pink'] = '#e7298a'
def update_color_cycle():
    mpl.rcParams.update({'axes.color_cycle': [colors['purple'],
                                              colors['orange'],
                                              colors['green']]})

def plot_smoothed(time, data, avg_period, ax=None, **plot_args):
    """Plot a time series smoothed with a cosine kernel

    Parameters:
        time        Array of equally spaced time points
        data        Array of data to plot at those times
        avg_period  Number of time points to average over
                    (width of the averaging kernel)
        ax          Axes to plot into (default: current/new Axes)
        plot_args   Any additional arguments to pass to pyplot.plot()
    """
    print("Averaging period: {:g} fs".format(avg_period * (time[1] - time[0])))
    #krnl = np.ones(avg_period) / avg_period
    krnl = signal.cosine(avg_period)
    krnl = krnl / np.sum(krnl)
    data_ma = np.convolve(data, krnl, mode='valid')
    if ax is None:
        ax = pyplot.gca()
    ax.plot(time[avg_period/2:-(avg_period/2 - 1)], data_ma, **plot_args)

def thin_points(data, r=None, nmax=1, density=None, len_scale=0.01):
    """Thin a set of data points down to some target density.

    Uses "poisonous seedling" thinning algorithm: Pick a point at
    random, and if it has more than a certain number of neighbours
    within some predefined radius, randomly remove neighbours until the
    number is reached.  Repeat in turn for all the remaining points.

    Parameters:
    data        Data matrix, size (N,d) - N is num points, d is problem
                dimensionality
    Specify either:
    density     Target density, maximum density in thinned set
        with optional:
    len_scale   Length scale of density averaging, as a fraction of the
                range in the data (maximum range along any dimension).
                Default 0.01.
    or:
    r           Absolute length scale of distance averaging (radius of
                sphere within which neighbours are counted)
        with optional:
    nmax        Maximum number of neighbours for a point (default 1)

    Note that the default of nmax=1 ensures that no two points are
    closer than r, thinning points to a maximum (approximately uniform)
    density of 1/r^dim

    NB density only implemented in two dimensions (i.e. d==2).

    """
    dim = data.shape[1]
    if density is None and r is None:
        raise ValueError("Must specify r if not using density")
    elif density is not None and dim != 2:
        raise NotImplementedError("Density not implemented for dimension != 2")
    elif density is not None:
        # Choose nmax and r such that r is about 1/100 of the range of
        # the data, but make nmax an integer.
        data_range = np.max(np.max(data, axis=0) - np.min(data, axis=0))
        nmax = np.ceil(density * np.pi * (data_range * len_scale)**2)
        r = np.sqrt(nmax / (density * np.pi))
    # Technically the leafsize should be equal to the ratio of volumes of the
    # d-cube to the d-sphere times nmax, but I think this is good enough.
    neightree = KDTree(data, leafsize=max(20, int(nmax) * 2))
    deleted = np.zeros((data.shape[0], ), dtype=bool)
    point_weight = np.ones((data.shape[0], ), dtype=np.float32)
    idces = np.arange(data.shape[0])
    np.random.shuffle(idces)
    for data_idx in idces:
        if deleted[data_idx]:
            continue
        # TODO adjust eps for efficiency?
        neighbours = neightree.query_ball_point(data[data_idx,:], r, eps=0.0)
        neighbours.remove(data_idx)
        if nmax == 1:
            del_idces = neighbours
        else:
            del_idces = np.random.choice(neighbours,
                                         len(neighbours) - nmax + 1)
        deleted[del_idces] = True
        point_weight[data_idx] += np.sum(point_weight[del_idces])
        point_weight[del_idces] = 0
    return data[~deleted,:], point_weight[~deleted]

def scatter_thin_points(data_thin, weights, method='size', ax=None,
                        **plot_args):
    """Make a scatterplot with a thinned set of points

    Paramters:
    data_thin   Thinned data, e.g. from thin_points()
    weights     Weights of the thinned data, e.g. from thin_points()
                (represents how many real points each thin point
                represents
    method      How to represent the weights; options are 'size',
                'color', and 'alpha'
    ax          Existing axis to plot into, if any
    Any remaining keyword arguments are passed to pyplot.scatter().

    """
    if ax is None:
        ax = pyplot.gca()
    if method == 'size':
        # TODO An automatic determination of the base size based on the
        #      actual plot size (in pixels) might be nice...
        s_scale = plot_args.pop('s', 2)
        # Retain the border (linewidth) to make it easier to see points' sizes
        pl = ax.scatter(data_thin[:,0], data_thin[:,1],
                        s=(weights * s_scale), **plot_args)
    elif method == 'color':
        plot_args.pop('c', None)
        pl = ax.scatter(data_thin[:,0], data_thin[:,1], linewidth=0, c=weights,
                        **plot_args)
    elif method == 'alpha':
        alpha_min = 0.0
        c_base = plot_args.pop('c', 'k')
        c_base = matplotlib.colors.colorConverter.to_rgba(c_base)
        c_array = np.zeros((len(weights), 4))
        c_array[:,:] = c_base
        c_array[:,3] = 1.0 + (alpha_min - 1.0) * np.exp(-1.0 * weights)
        pl = ax.scatter(data_thin[:,0], data_thin[:,1], linewidth=0, c=c_array,
                        **plot_args)
    else:
        raise ValueError("Unknown representation method " + method)
    return pl

def scatter_mark_outliers(xs, ys, xlim, ylim, ax=None, **plot_args):
    """Mark outliers at the edges of a scatterplot

    Paramters:
    xs      List of x-coordinates of points, passed to pyplot.scatter()
    ys      List of y-coordinates of points, passed to pyplot.scatter()
    xlim    Limits of x-axis, as a tuple (xmin, xmax)
    ylim    Limits of y-axis, as a tuple (ymin, ymax)
    ax      Existing axis to plot into, if any

    Any other keyword arguments are passed to pyplot.scatter().

    The outliers are marked at the edge of the graph as large red dots,
    by default size 40, color 'r'.  If a label was present, the outliers
    will be labelled as the same with " (outliers)" appended.

    Returns the collection of scatterplots - (top, bottom, left, right).

    """

    if ax is None:
        ax = pyplot.gca()
    outl_top = xs[(ys > ylim[1]) & (xs >= xlim[0]) & (xs <= xlim[1])]
    outl_bot = xs[(ys < ylim[0]) & (xs >= xlim[0]) & (xs <= xlim[1])]
    # Let out_left and outl_right handle the (literal) corner cases
    outl_left = np.array(ys[(xs < xlim[0])])
    outl_right = np.array(ys[(xs > xlim[1])])
    outl_left[outl_left < ylim[0]] = ylim[0]
    outl_left[outl_left > ylim[1]] = ylim[1]
    outl_right[outl_right < ylim[0]] = ylim[0]
    outl_right[outl_right > ylim[1]] = ylim[1]
    if 'linewidth' not in plot_args:
        plot_args['linewidth'] = 0
    if 'c' not in plot_args:
        plot_args['c'] = 'r'
    if 's' not in plot_args:
        plot_args['s'] = 40
    plot_args['clip_on'] = False
    plot_args['label'] = (plot_args['label'] + ' (outliers)'
                          if 'label' in plot_args else '_')
    op_t = ax.scatter(outl_top, ylim[1] * np.ones_like(outl_top), **plot_args)
    plot_args['label'] = '_'
    op_b = ax.scatter(outl_bot, ylim[0] * np.ones_like(outl_bot), **plot_args)
    op_l = ax.scatter(
        xlim[0] * np.ones_like(outl_left), outl_left, **plot_args)
    op_r = ax.scatter(
        xlim[1] * np.ones_like(outl_right), outl_right, **plot_args)
    # Not typically necessary, but just ensure the outlier markers end
    # up at the edge of the plot
    ax.set_xlim(*xlim)
    ax.set_ylim(*ylim)
    return op_t, op_b, op_l, op_r

def scatter_outliers_size(data_thin, weights, xlim, ylim, ax=None, c_out='r',
                          s_out=40, **plot_args):
    """Do scatterplot with outliers, encoding the given weights as size.

    NB This doesn't change the size of the outlier markers - it might be
    possible by modifying the original escatter_outliers() function, but
    does it even make sense to encode the weights of outliers?

    Parameters:
    data_thin   Thinned data, e.g. from thin_points()
    weights     Weights of the thinned data, e.g. from thin_points()
                (represents how many real points each thin point
                represents
    xlim        Limits of x-axis, as a tuple (xmin, xmax)
    ylim        Limits of y-axis, as a tuple (ymin, ymax)
    ax          Existing axis to plot into, if any
    c_out       Colour to use to mark outliers, default 'r'
    s_out       Size for outlier markers, default 40 (px)
    Any remaining keyword arguments will only modify the properties of
    the central (non-outlier) scatterplot.

    Returns the original plot as well as the four outlier scatterplots -
    (orig, top, bottom, left, right)

    """
    pl = scatter_thin_points(data_thin, weights, ax=ax, **plot_args)
    outps = scatter_mark_outliers(
        *data_thin.T, xlim=xlim, ylim=ylim, ax=ax, c=c_out, s=s_out)
    return (pl,) + outps

ext: py

reason: import numpy,from scipy


text:
# Copyright 2021 NREL

# Licensed under the Apache License, Version 2.0 (the "License"); you may not
# use this file except in compliance with the License. You may obtain a copy of
# the License at http://www.apache.org/licenses/LICENSE-2.0

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

# See https://floris.readthedocs.io for documentation

# NOTE: This example requires an additional module, `streamlit`.
#
# TO USE THIS FROM COMMAND TYPE:
#    streamlit run interactive_visual.py

import matplotlib.pyplot as plt
import numpy as np

import streamlit as st
import floris.tools as wfct

# Fixed parameters
minSpeed = 6.0
maxSpeed = 8.0
D  = 126.
HH = 90.

# Options
ws = st.sidebar.slider("Wind Speed", 5., 10., 8., step=0.1)
wd = st.sidebar.slider("Wind Direction", 250, 290, 270, step=2)
yaw_1 = st.sidebar.slider("Yaw angle T1", -30, 30, 0, step=1)
x_loc = st.sidebar.slider("x normal plane intercept", 0., 35., 7., step=.25)
y_loc = st.sidebar.slider("y normal plane intercept", -2., 2., 0., step=0.1)
second_turbine  = st.sidebar.checkbox('Second Turbine?')
GCH = st.sidebar.checkbox('GCH')

# Initialize the FLORIS interface fi (legacy)
fi = wfct.floris_interface.FlorisInterface("../example_input.json")
fi.set_gch(GCH)
if second_turbine:
    fi.reinitialize_flow_field(
        wind_speed = ws, wind_direction=wd, layout_array=((0, 126 * 7), (0, 0))
        # wind_speed = ws, wind_direction=wd, layout_array=((0, 126 * 7, 126 * 14, 126 * 21, 126 * 28), (0, 0, 0, 0, 0))
    )
    fi.calculate_wake(yaw_angles=[yaw_1, 0])
    # fi.calculate_wake(yaw_angles=[yaw_1, 0, 0, 0, 0])
else:
    fi.reinitialize_flow_field(
        wind_speed = ws,wind_direction=wd, layout_array=([0], [0])
    )
    fi.calculate_wake(yaw_angles=[yaw_1])

# Initialize the FLORIS interface fi (legacy)
fi_gauss = wfct.floris_interface.FlorisInterface("gauss_version.json")
fi_gauss.set_gch(GCH)
if second_turbine:
    fi_gauss.reinitialize_flow_field(
        wind_speed = ws, wind_direction=wd, layout_array=((0, 126 * 7), (0, 0))
        # wind_speed = ws, wind_direction=wd, layout_array=((0, 126 * 7, 126 * 14, 126 * 21, 126 * 28), (0, 0, 0, 0, 0))
    )
    fi_gauss.calculate_wake(yaw_angles=[yaw_1, 0])
    # fi_gauss.calculate_wake(yaw_angles=[yaw_1, 0, 0, 0, 0])
else:
    fi_gauss.reinitialize_flow_field(
        wind_speed = ws,wind_direction=wd, layout_array=([0], [0])
    )
    fi_gauss.calculate_wake(yaw_angles=[yaw_1, 0])

st.write("# Horizontal planes")

# Horizontal plane

fig, axarr = plt.subplots(2,1)

hor_plane = fi.get_hor_plane()
ax = axarr[0]
wfct.visualization.visualize_cut_plane(
    hor_plane, ax=ax, minSpeed=minSpeed, maxSpeed=maxSpeed
)
ax.axhline(y_loc * D, color="w", ls="--", lw=1)
ax.axvline(x_loc*D, color="w", ls="--", lw=1)
ax.set_title("Legacy")

hor_plane_gauss = fi_gauss.get_hor_plane()
ax = axarr[1]
wfct.visualization.visualize_cut_plane(
    hor_plane_gauss, ax=ax, minSpeed=minSpeed, maxSpeed=maxSpeed
)
ax.axhline(y_loc * D, color="w", ls="--", lw=1)
ax.axvline(x_loc*D, color="w", ls="--", lw=1)
ax.set_title("Gauss")

st.write(fig)

# print(hor_plane_gauss.df.x2.unique())

# Look at flow center line
x_points = np.arange(-100, D*16, 1)
y_points = np.ones_like(x_points) * y_loc * D
z_points = np.ones_like(x_points) * HH

# Get the values
flow_points = fi.get_set_of_points(x_points, y_points, z_points).sort_values('x')
flow_points_gauss = fi_gauss.get_set_of_points(x_points, y_points, z_points).sort_values('x')

# Compare wind speeds at centerline
fig, ax = plt.subplots()
ax.plot(flow_points.x/D, flow_points.u,color='b',label='Legacy')
ax.plot(flow_points_gauss.x/D, flow_points_gauss.u,color='g',label='Gauss')
ax.grid(True)
ax.legend()
ax.set_title('Center line wind speed')
st.write(fig)

# Compare the wake profile at x
# Show the velocity profile using a rotor sweep
cross_plane = fi.get_cross_plane(x_loc * D)
x1_locs, v_array = wfct.cut_plane.wind_speed_profile(
    cross_plane, D / 2, HH, resolution=100, x1_locs=None
)
cross_plane_gauss = fi_gauss.get_cross_plane(x_loc * D)
x1_locs, v_array_gauss = wfct.cut_plane.wind_speed_profile(
    cross_plane_gauss, D / 2, HH, resolution=100, x1_locs=None
)

fig, ax = plt.subplots()
ax.plot(x1_locs, v_array, color="b",label="Legacy")
ax.plot(x1_locs, v_array_gauss, color="r",label="Gauss")
ax.grid(True)
ax.set_title("Rotor-averaged velocity profile")
ax.set_ylabel("Rotor-averaged wind speed (m/s)")
ax.set_xlabel("Center of rotor disk (m)")
ax.legend()
st.write(fig)

# # Legacy

# Powers downstream
x_array = np.arange(D,D*12,D/2)
leg_res = []
gauss_res = []

for x in x_array:

    fi.reinitialize_flow_field(
        wind_speed = ws, wind_direction=wd, layout_array=([0, x], [0, y_loc * D])
    )
    fi.calculate_wake(yaw_angles=[yaw_1, 0])
    leg_res.append(fi.get_turbine_power()[1])

    fi_gauss.reinitialize_flow_field(
        wind_speed = ws, wind_direction=wd, layout_array=([0, x], [0, y_loc * D])
    )
    fi_gauss.calculate_wake(yaw_angles=[yaw_1, 0])
    gauss_res.append(fi_gauss.get_turbine_power()[1])

fig, ax = plt.subplots()
ax.plot(x_array/D, leg_res, color="b",label="Legacy")
ax.plot(x_array/D, gauss_res, color="r",label="Gauss")
ax.grid(True)
ax.set_title("Power behind turbine")
ax.set_ylabel("Power")
ax.set_xlabel("Distance downstream")
ax.legend()
st.write(fig)

reason: import numpy


text:
from PyAstronomy.pyasl import KeplerEllipse
from PyAstronomy import funcFit as fuf
from numpy import pi
import numpy as np
from PyAstronomy.pyaC import pyaErrors as PE

class KeplerEllipseModel(fuf.OneDFit):
    """
    A model of a Keplerian orbit.

    This class uses the *KeplerEllipse* from the PyA's pyasl
    to calculate a Keplerian orbit. It may be used to fit
    complete 3d position or velocity information on the orbit;
    any individual axes may also be selected.

    The constructor allows to specify *relevant axes*, which
    are those axes considered in the calculation. The actual
    (technical) model is, however, only one-dimensional. The values
    returned by evaluate have the order
    a1, b1, c1, a2, b2, c3, ... . Where a, b, and c represent
    the first, second, and third axis and the number specifies
    the data point. Note that in this case, the resulting
    model has not the same number of points as the time
    axis.

    *Fit parameters*
      -  `a`    - The semi-major axis (same units as the data)
      - `per`   - The period (same time units as data)
      - `e`     - The eccentricity
      - `tau`   - Time of periapsis passage (same time units as data)
      - `Omega` - Longitude of the ascending node [deg]
      - `w`     - Argument of periapsis [deg]
      - `i`     - Inclination angle [deg]

    Parameters
    ----------
    relevantAxes : string
        A string containing any combination of x, y, and z.
        The string specifies the axes (and their order) to
        be considered in the calculations.
    mode : string, {"pos", "vel"}
        Determines whether the output is positions or
        velocities. In this case, the units are determined
        by the units of major axis and time (e.g., AU per day).
    """

    def __init__(self, relevantAxes="xyz", mode="pos"):
        self.ke = KeplerEllipse(1.0, 1.0)
        fuf.OneDFit.__init__(self, ["a", "per", "e", "tau", "Omega", "w", "i"])
        self["a"] = 1.0
        self["per"] = 1.0
        # Which axes to consider?
        # x=0, y=1, z=2
        self.axes = ()
        for i, axis in enumerate("xyz"):
            if axis in relevantAxes:
                self.axes += (i,)
        # Save the mode
        if not mode in ["pos", "vel"]:
            raise(PE.PyAValError("Unknown mode: " + str(mode), \
                                 where="KeplerEllipseModel", \
                                 solution="Choose either 'pos' or 'vel'."))
        self._mode = mode

    def evaluate(self, t):
        """
        Calculates and returns model according to the
        current parameter values.

        Although more than one axis may be relevant
        the output will be one dimensional. If, e.g.,
        the relevant axes are x and y, the order of
        the output will be x0, y0, x1, y1, ... .

        Parameters
        ----------
        t : array
            Times at which to evaluate the model.
        """
        self.ke.i = self["i"]
        self.ke.w = self["w"]
        self.ke.Omega = self["Omega"]
        self.ke.e = abs(self["e"])
        self.ke.a = self["a"]
        self.ke.per = self["per"]
        self.ke.tau = self["tau"]
        self.ke.n = 2.0 * pi / self["per"]
        # Get the data pertaining to the relevant axes
        if self._mode == "pos":
            result = self.ke.xyzPos(t)[::, self.axes]
        else:
            result = self.ke.xyzVel(t)[::, self.axes]
        # Reshape to 1d form
        # If the relevant axes are x and y, the order
        # of the output will be x0, y0, x1, y1, x2, y2, ...
        result = result.reshape(result.size)
        return result

class KeplerRVModel(fuf.OneDFit):
    """
    A model of a Keplerian orbit.

    This class uses the *KeplerEllipse* from the PyA's pyasl
    to calculate radial velocities for a Keplerian orbit. All calculations are
    based on the negligible companion mass hypothesis (i.e., this model has been
    implemented with planetary systems in mind).

    .. note:: Any planet with zero period will be ignored.

    .. note:: ? is a placeholder for an integer larger zero, indicating the number of
              the planet (the total number is controlled by the `mp` keyword).

    *Fit parameters*

      - `per?`   - The period (same time units as data)
      - `e?`     - The eccentricity
      - `tau?`   - Time of periapsis passage (same time units as data)
      - `w?`     - Argument of periapsis [deg]
      - `K?`     - Semi-amplitude of radial velocity (same units as data)
      - `mstar` - Stellar mass in solar masses. This parameter is usually not fitted.
                  It may be used to take into account the uncertainty on stellar mass
                  in Bayesian (MCMC) analysis.

    *Derived parameters (not to be fitted)*
      -  `a?`        - The semi-major axis in AU
      -  `MA?`       - Mean anomaly corresponding to time of first data point [deg]
      -  `msini?`    - Minimum mass (msini) in Jupiter masses

    Parameters
    ----------
    mp : int, optional
        The number of planets considered in the model. Default is one. Note that
        signals are added, i.e., no interaction is taken into account.
    deg : int, optional
        Default is zero (i.e., a constant). The degree of a polynomial used to represent
        a systematic (non-periodic) evolution in the data.
    msun : float, optional
        Solar mass [kg]
    mJ : float, optional
        Jupiter mass [kg]
    au : float, optional
        Astronomical unit [m]
    """

    def __init__(self, mp=1, deg=0, msun=1.988547e30, mJ=1898.6e24, au=1.49597870700e11):
        self._msun = msun
        self._mJ = mJ
        self._au = au
        self._deg = deg
        self._mp = mp

        self.poly = fuf.PolyFit1d(deg)
        # Kepler models for individual planets
        self.kems = [KeplerEllipseModel(relevantAxes="z", mode="vel") for _ in range(1,mp+1)]

        pars = []
        for m in range(1, mp+1):
            # Independent parameters
            pars.extend([pn+str(m) for pn in ["K", "per", "e", "tau", "w"]])
            # Dependent parameters
            pars.extend([pn+str(m) for pn in ["msini", "a", "MA"]])
        # There is only one stellar mass
        pars.append("mstar")
        # Use parameters from polynomial
        pars.extend(self.poly.availableParameters())

        fuf.OneDFit.__init__(self, pars)
        self.setRootName("KeplerRVModel")
        # Dummy SMA and inclination
        for m in range(mp):
            self.kems[m]["a"] = 1.0
            self.kems[m]["i"] = 90.0
        # Default stellar mass = 0 solar mass (force thoughtful input)
        self["mstar"] = 0.0
        self.setRestriction({"mstar":[0,None]})

        for m in range(1, 1+mp):
            self.setRestriction({"K"+str(m):[0.,None], "per"+str(m):[0.,None], "e"+str(m):[0,1]})

    def _dtr(self, d):
        """ Convert degrees into rad """
        return d/180.*np.pi

    def evaluate(self, t):
        if self["mstar"] == 0.0:
            raise(PE.PyAValError("Stellar mass has to be specified before evaluation.", \
                                 where="KeplerRVModel"))

        # Take into account polynomial
        for i in range(self._deg+1):
            p = "c" + str(i)
            self.poly[p] = self[p]
        rvnew = self.poly.evaluate(t)

        for m in range(1, self._mp+1):
            # Name add
            nadd = str(m)
            if self["per"+nadd] == 0.0:
                continue
            for p in ["per", "e", "tau", "w"]:
                # Collect parameters from current model
                self.kems[m-1][p] = self[p+nadd]

                rvmodel = self.kems[m-1].evaluate(t)
            # Prevent dividing by zero
            imax = np.argmax(np.abs(rvmodel))
            cee = np.cos(self.kems[m-1].ke.trueAnomaly(t[imax]) + self._dtr(self.kems[m-1]["w"])) \
                  + self["e"+nadd]*np.cos(self._dtr(self["w"+nadd]))
            rvnew += rvmodel / rvmodel[imax] * cee * self["K"+nadd]

            self["msini"+nadd] = self._getmsini(nadd)
            self["MA"+nadd] = self._MA(m)
            self["a"+nadd] = self._geta(nadd)

        return rvnew

    def _getmsini(self, nadd):
        """
        Get msini
        """
        msini = self["K"+nadd] * ((self["per"+nadd]*86400.) * (self["mstar"]*self._msun)**2 / (2.*np.pi*6.67408e-11))**(1./3.) \
                * np.sqrt(1. - self["e"+nadd]**2)
        return msini/self._mJ

    def _MA(self, m):
        """
        Get mean anomaly corresponding to first data point [deg]
        """
        return self.kems[m-1].ke.meanAnomaly(0.0)/np.pi*180. % 360.0

    def _geta(self, nadd):
        """
        Get SMA in AU
        """
        return ((self["per"+nadd]*86400.)**2 * 6.67408e-11 * (self["mstar"] * self._msun) / (4.*np.pi**2) )**(1./3.) / self._au

max_stars_repo_name: mirofedurco/PyAstronomy

max_issues_repo_name: mirofedurco/PyAstronomy

max_forks_repo_name: mirofedurco/PyAstronomy

reason: import numpy,from numpy


text:
from flask_ngrok import run_with_ngrok
from flask import Flask, render_template, request
import convlab2
from convlab2.nlu.jointBERT.multiwoz import BERTNLU
from convlab2.nlu.milu.multiwoz import MILU
from convlab2.dst.rule.multiwoz import RuleDST
from convlab2.policy.rule.multiwoz import RulePolicy
from convlab2.nlg.template.multiwoz import TemplateNLG
from convlab2.dialog_agent import PipelineAgent, BiSession
from convlab2.evaluator.multiwoz_eval import MultiWozEvaluator
from pprint import pprint
import random
import numpy as np
import torch

app = Flask(__name__)
run_with_ngrok(app)   #starts ngrok when the app is run
@app.route("/")
def home():
    return "
Running Flask on Google Colab!
"
#define app routes
@app.route("/chat")
def index():
    return render_template("index.html")

@app.route("/get")
#function for the bot response
def get_bot_response():
    userText = request.args.get('msg')
    sys_nlu = BERTNLU()
    sys_dst = RuleDST()
    sys_policy = RulePolicy()
    sys_nlg = TemplateNLG(is_user=False)
    sys_agent = PipelineAgent(sys_nlu, sys_dst, sys_policy, sys_nlg, name='sys')
    return str(sys_agent.response(userText))

@app.route('/process',methods=['POST'])
def process():
    user_input = request.form['user_input']
    sys_nlu = BERTNLU()
    sys_dst = RuleDST()
    sys_policy = RulePolicy()
    sys_nlg = TemplateNLG(is_user=False)
    sys_agent = PipelineAgent(sys_nlu, sys_dst, sys_policy, sys_nlg, name='sys')
    bot_response = str(sys_agent.response(user_input))
    print("Cherry: "+bot_response)
    return render_template('index.html',user_input=user_input,bot_response=bot_response)

if __name__ == "__main__":
    app.run()


text:
# from skimage import io
# import numpy as np
from collections import defaultdict
from glob import glob
from os import path

from config import get_config

# import matplotlib.pyplot as plt
# import pickle

# plt.switch_backend('agg')

def main():
    save_fp = 'info.txt'
    cfg = get_config()
    patch_path = cfg.patch_path

    patch_txt_fps = glob(path.join(patch_path, '*.txt'))

    slide_cnt = defaultdict(lambda: defaultdict(int))
    cnt = defaultdict(int)

    for each in patch_txt_fps:
        base_name = path.basename(each)
        img_id = int(base_name.split('.')[0])
        print(f'Reading => {each}')
        with open(each, 'rt', encoding='utf-8') as in_file:
            for line in in_file:
                line_sp = line.strip().split(',')
                label = int(line_sp[-1])
                slide_cnt[img_id][label] += 1
                cnt[label] += 1

    with open(save_fp, 'wt', encoding='utf-8') as out_file:
        out_file.write(f'Total:{cnt[0]+cnt[1]}, Pos:{cnt[1]}, Neg:{cnt[0]}\n')
        out_file.write('=' * 10 + '\n')

        slide_cnt_kv = sorted(slide_cnt.items(), key=lambda key: key[0], reverse=False)

        pos_count, neg_cnt = 0, 0
        for k, v in slide_cnt_kv:
            out_file.write(f'Slide@{k}:{v[1]+v[0]}, Pos:{v[1]}, Neg:{v[0]}\n')
            pos_count += v[1]
            neg_cnt += v[0]
            if k == 100:
                out_file.write(f'Pos_cnt:{pos_count}, Neg_cnt:{neg_cnt}\n')
                out_file.write('=' * 10 + '\n')
                pos_count = 0
                neg_cnt = 0
        out_file.write(f'Pos_cnt:{pos_count}, Neg_cnt:{neg_cnt}\n')
        out_file.write('=' * 10 + '\n')

    print(f'Done!\n, All:{cnt[0]+cnt[1]}, Pos:{cnt[1]}, Neg:{cnt[0]}\n save in {path.abspath(save_fp)}')

if __name__ == '__main__':
    main()

max_stars_repo_name: steermomo/ACDC

max_issues_repo_name: steermomo/ACDC

max_forks_repo_name: steermomo/ACDC


text:
from flask import Flask, render_template, request
# from models import MobileNet
from model import Net
import os
from math import floor

import torch
# import torch.nn as nn
# import torch.nn.functional as F
import torch.optim as optim
from torch.autograd import Variable
# from torchvision import datasets, transforms
from torchvision import transforms
import numpy as np

from PIL import Image

app = Flask(__name__)

app.config['UPLOAD_FOLDER'] = 'uploads'

########################### START - model stuff ###########################

use_cuda = False; #not args.no_cuda and torch.cuda.is_available()
torch.manual_seed(1); #args.seed
device = torch.device("cuda" if use_cuda else "cpu");
kwargs = {'num_workers': 1, 'pin_memory': True} if use_cuda else {};

lr = 0.01;
momentum = 0.5;

model = Net().to(device);
optimizer = optim.SGD(model.parameters(), lr=lr, momentum=momentum);

model_path = "model.pth"
model.load_state_dict(torch.load(model_path))
model.eval()

########################### END - model stuff ###########################

@app.route('/')
def index():
    return render_template('index.html')

@app.route('/about')
def about():
    return render_template('about.html')

def pil_loader(path):
    """Load images from /eval/ subfolder, convert to greyscale and resized it as squared"""
    with open(path, 'rb') as f:
        with Image.open(f) as img:
            sqrWidth = np.ceil(np.sqrt(img.size[0]*img.size[1])).astype(int)
            return img.convert('L').resize((sqrWidth, sqrWidth))

@app.route('/predict', methods=['POST'])
def predict():
    if request.method == 'POST':
        f = request.files['file']
        saveLocation = f.filename
        f.save(saveLocation)

        image_path = saveLocation

        # input_image = Image.open(image_path)
        input_image = pil_loader(image_path)
        preprocess = transforms.Compose([
            # transforms.Resize(256),
            # transforms.CenterCrop(224),
            # transforms.ToTensor(),
            # transforms.Normalize(mean=[0.485, 0.456, 0.406], std=[0.229, 0.224, 0.225]),
            transforms.Resize(28),
            transforms.CenterCrop(28),
            transforms.ToTensor(),
            transforms.Normalize((0.1307,), (0.3081,))
        ])

        input_tensor = preprocess(input_image)
        # create a mini-batch as expected by the model
        input_batch = input_tensor.unsqueeze(0)
        #output = self.model(input_batch)

        # input = Variable(input_batch)
        # input = input.to(device)

        #inference, confidence = model.infer(saveLocation)

        # single_loaded_img = test_loader.dataset.data[test_data_loader_index]
        # #print(single_loaded_img)
        # single_loaded_img = single_loaded_img.to(device)
        # single_loaded_img = single_loaded_img[None, None]
        # #single_loaded_img = single_loaded_img.type('torch.DoubleTensor')
        # single_loaded_img = single_loaded_img.type('torch.FloatTensor')

        #raw_prediction = model(input)
        raw_prediction = model(input_batch)

        #maxed_prediction = raw_prediction.max(1, keepdim=True)[1]
        maxed_prediction = raw_prediction.argmax(dim=1, keepdim=True).item()
        # maxed_prediction = maxed_prediction.item()
        print("maxed_prediction:")
        print(maxed_prediction)
        print("")

        output = "hello world"
        return str(maxed_prediction)

# from torchvision import transforms
# class ConvNet
    # def infer(self, image_path):
    #     input_image = Image.open(image_path)
    #     preprocess = transforms.Compose([
    #         transforms.Resize(256),
    #         transforms.CenterCrop(224),
    #         transforms.ToTensor(),
    #         transforms.Normalize(mean=[0.485, 0.456, 0.406], std=[0.229, 0.224, 0.225]),
    #     ])
    #     input_tensor = preprocess(input_image)

    #     # create a mini-batch as expected by the model
    #     input_batch = input_tensor.unsqueeze(0)

    #     # move the input and model to GPU for speed if available
    #     if torch.cuda.is_available():
    #         input_batch = input_batch.to('cuda')
    #         self.model.to('cuda')

    #     with torch.no_grad():
    #         output = self.model(input_batch)

    #     # The output has unnormalized scores. To get probabilities, you can run a softmax on it.
    #     output = torch.nn.functional.softmax(output[0], dim=0)
    #     confidence, index = torch.max(output, 0)

    #     return (self.classes[index.item()], confidence.item())

# @app.route('/infer', methods=['POST'])
# def success():
#     if request.method == 'POST':
#         f = request.files['file']
#         saveLocation = f.filename
#         f.save(saveLocation)
#         inference, confidence = model.infer(saveLocation)
#         # make a percentage with 2 decimal points
#         confidence = floor(confidence * 10000) / 100
#         # delete file after making an inference
#         os.remove(saveLocation)
#         # respond with the inference
#         return render_template('inference.html', name=inference, confidence=confidence)

if __name__ == '__main__':
    app.debug = True
    port = int(os.environ.get("PORT", 80))
    app.run(host='0.0.0.0', port=port, debug=True)

max_stars_repo_name: pandorasNox/deepo

max_issues_repo_name: pandorasNox/deepo

max_forks_repo_name: pandorasNox/deepo


text:
import os

import numpy as np
import pytest
import xarray as xr

from climpred.constants import CLIMPRED_ENSEMBLE_DIMS
from climpred.preprocessing.shared import (
    load_hindcast,
    rename_SLM_to_climpred_dims,
    rename_to_climpred_dims,
    set_integer_time_axis,
)

on_mistral = False
try:
    host = os.environ["HOSTNAME"]
    for node in ["mlogin", "mistralpp"]:
        if node in host:
            on_mistral = True
            from climpred.preprocessing.mpi import get_path
except KeyError:
    pass

# check for intake_esm to be installed
try:
    import intake
    import intake_esm

    print(intake_esm.__version__)
    intake_esm_loaded = True
except ImportError:
    intake_esm_loaded = False

def preprocess_1var(ds, v="global_primary_production"):
    return ds[v].to_dataset(name=v).squeeze()

@pytest.mark.mistral
@pytest.mark.skipif(not on_mistral, reason="requires to be on mistral.dkrz.de")
@pytest.mark.parametrize(
    "inits,members",
    [(range(1961, 1964), range(3, 6)), (range(1970, 1972), range(1, 3))],
)
def test_load_hindcast(inits, members):
    """Test that `load_hindcast` loads the appropriate files."""
    actual = load_hindcast(
        inits=inits, members=members, preprocess=preprocess_1var, get_path=get_path,
    )
    assert isinstance(actual, xr.Dataset)
    assert (actual["init"].values == inits).all()
    assert (actual["member"].values == members).all()
    assert "global_primary_production" in actual.data_vars
    assert len(actual.data_vars) == 1

@pytest.mark.mistral
@pytest.mark.skipif(not on_mistral, reason="requires to be on mistral.dkrz.de")
@pytest.mark.skipif(not intake_esm_loaded, reason="requires intake_esm to be installed")
def test_climpred_pre_with_intake_esm():
    """Test that `preprocess` including `set_integer_time_axis` enables concatination
    of all hindcast into one xr.object."""
    col_url = "/home/mpim/m300524/intake-esm-datastore/catalogs/mistral-cmip6.json"
    col = intake.open_esm_datastore(col_url)
    # load 2 members for 2 inits from one model
    query = dict(
        experiment_id=["dcppA-hindcast"],
        table_id="Amon",
        member_id=["r1i1p1f1", "r2i1p1f1"],
        dcpp_init_year=[1970, 1971],
        variable_id="tas",
        source_id="MPI-ESM1-2-HR",
    )
    cat = col.search(**query)
    cdf_kwargs = {"chunks": {"time": 12}, "decode_times": False}

    def preprocess(ds):
        # extract tiny spatial and temporal subset
        ds = ds.isel(lon=[50, 51, 52], lat=[50, 51, 52], time=np.arange(12 * 2))
        # make time dim identical
        ds = set_integer_time_axis(ds)
        return ds

    dset_dict = cat.to_dataset_dict(cdf_kwargs=cdf_kwargs, preprocess=preprocess)
    # get first dict value
    ds = dset_dict[list(dset_dict.keys())[0]]
    assert isinstance(ds, xr.Dataset)
    ds = rename_to_climpred_dims(ds)
    # check for all CLIMPRED_DIMS
    for c in CLIMPRED_ENSEMBLE_DIMS:
        assert c in ds.coords
    # check for requested dimsizes
    assert ds.member.size == 2
    assert ds.init.size == 2

def test_rename_SLM(da_SLM):
    """Check that dimensions in input are renamed by rename_SLM_to_climpred_dims."""
    da_renamed = rename_SLM_to_climpred_dims(da_SLM)
    for dim in CLIMPRED_ENSEMBLE_DIMS:
        assert dim in da_renamed.dims
    for dim in da_SLM.dims:
        assert dim not in da_renamed.dims

def test_rename_climpred_dims(da_dcpp):
    """Check that dimensions in input are renamed by rename_to_climpred_dims."""
    da_renamed = rename_to_climpred_dims(da_dcpp)
    for dim in CLIMPRED_ENSEMBLE_DIMS:
        assert dim in da_renamed.dims
    for dim in da_dcpp.dims:
        assert dim not in da_renamed.dims

max_stars_repo_name: bradyrx/climpred

max_issues_repo_name: bradyrx/climpred

max_forks_repo_name: bradyrx/climpred


text:
import os
import logging
import numpy as np
import pandas as pd
from sklearn.model_selection import train_test_split
import opticl
from pyomo import environ
from pyomo.environ import *
import doframework as dof

@dof.resolve
def predict_optimize_opticl(data: np.array, constraints: np.array, **kwargs):
    '''
    OptiCL predict-optimize algorithm.

            Parameters:
                    data (np.array): data for regression [possibly including unfeasible data points].
                    constraints (np.array): array of coefficients of equations that define the feasibility region.

            Returns:
                    optimum: the location of a predicted optimal solution.
                    predicted_value: the predicted optimal value.
                    model: the regression model used for prediction.

    '''

    logger_name = kwargs['logger_name'] if 'logger_name' in kwargs else None
    is_raised = kwargs['is_raised'] if 'is_raised' in kwargs else False
    alg_list = kwargs['alg_list'] if 'alg_list' in kwargs else ['linear','rf','svm','cart','gbm','mlp']
    task_type = kwargs['task_type'] if 'task_type' in kwargs else 'continuous'
    is_trust_region = kwargs['is_trust_region'] if 'is_trust_region' in kwargs else True
    solver = kwargs['solver'] if 'solver' in kwargs else None

    tolerance = 1e-8

    try:

        data_feasible = data[np.all(np.pad(data[:,:-1],((0,0),(0,1)),constant_values=1) @ constraints.T <= -tolerance, axis=1)]

        X = pd.DataFrame(data_feasible[:,:-1]).add_prefix('x')
        y = pd.DataFrame(data_feasible[:,-1:]).add_prefix('y')

        X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=1)

        mfiles = {}
        models = {}
        perfs = pd.DataFrame()
        response_vars = y_train.columns

        for var in response_vars:

            for alg in alg_list:

                alg_run = 'rf_shallow' if alg == 'rf' else alg

                model_path = os.path.join('results','{}'.format(alg))
                model_file = os.path.join(model_path,'{}_model'.format(var))

                m, perf = opticl.run_model(X_train, y_train[var], X_test, y_test[var], alg_run, var,
                                           task = task_type, cv_folds = 5, parameter_grid = None,
                                           save_path = model_file, save = False)

                constraintL = opticl.ConstraintLearning(X_train, y_train, m, alg)
                constraint_add = constraintL.constraint_extrapolation(task_type)
                mfiles[model_file] = constraint_add
                models[alg] = m

                perf['outcome'] = var
                perf['alg'] = alg
                perfs = pd.concat([perfs,perf],ignore_index=True)

        model_master = opticl.model_selection(perfs, objectives_embed = {var:1})
        model_master['lb'] = None
        model_master['ub'] = None
        model_chosen = model_master.loc[model_master.outcome == response_vars[0],'model_type'][0]

        model_pyo = ConcreteModel()

        N = X_train.columns
        model_pyo.x = Var(N, domain=Reals)

        for i,eqn in enumerate(constraints):
            model_pyo.add_component('constr_known{}'.format(i), Constraint(expr=sum(model_pyo.x[var]*eqn[j] for j,var in enumerate(N)) <= -eqn[-1]+tolerance))

        model_pyo.OBJ = Objective(expr=0, sense=minimize)

        final_model_pyo = opticl.optimization_MIP(model_pyo,
                                                  model_pyo.x,
                                                  model_master,
                                                  mfiles,
                                                  X_train,
                                                  tr=is_trust_region)

        opt = SolverFactory(solver)
        _ = opt.solve(final_model_pyo)
        x_solution = getattr(final_model_pyo, 'x')
        optimum = np.array([x_solution[i].value for i in N])
        predicted_value = final_model_pyo.OBJ()

    except ValueError as e:
        if logger_name:
            log = logging.getLogger(logger_name)
            log.error('ValueError: Prediction optimization failed.')
            log.error(e)
        else:
            print('ValueError: Prediction optimization failed.')
            print(e)
        if is_raised:
            raise
        else:
            return None, None, None

    except TypeError as e:
        if logger_name:
            log = logging.getLogger(logger_name)
            log.error('TypeError: Prediction optimization failed.')
            log.error(e)
        else:
            print('TypeError: Prediction optimization failed.')
            print(e)
        if is_raised:
            raise
        else:
            return None, None, None

    except Exception as e:
        if logger_name:
            log = logging.getLogger(logger_name)
            log.error('Exception: Prediction optimization failed.')
            log.error(e)
        else:
            print('Exception: Prediction optimization failed.')
            print(e)
        if is_raised:
            raise
        else:
            return None, None, None

    else:

        return optimum, predicted_value, models[model_chosen]

if __name__ == '__main__':

    configs_path = os.getcwd()
    configs_file = 'configs.yaml'

    dof.run(predict_optimize_opticl,
            os.path.join(configs_path,configs_file),
            alg_list=['linear','rf','svm','cart','gbm'],
            solver='cplex',
            objectives=1, datasets=2, after_idle_for=400, distribute=True, alg_num_cpus=4)

max_stars_repo_name: IBM/doframework

max_issues_repo_name: IBM/doframework

max_forks_repo_name: IBM/doframework

reason: import numpy,from pyomo


text:
# -*- coding: utf-8 -*-
from numpy import abs, cos, exp, pi, prod, sin, sqrt, sum
from .go_benchmark import Benchmark

class Ackley01(Benchmark):

    r"""
    Ackley01 objective function.

    The Ackley01 [1]_ global optimization problem is a multimodal minimization
    problem defined as follows:

    .. math::

        f_{\text{Ackley01}}(x) = -20 e^{-0.2 \sqrt{\frac{1}{n} \sum_{i=1}^n
         x_i^2}} - e^{\frac{1}{n} \sum_{i=1}^n \cos(2 \pi x_i)} + 20 + e

    Here, :math:`n` represents the number of dimensions and :math:`x_i \in
    [-35, 35]` for :math:`i = 1, ..., n`.

    *Global optimum*: :math:`f(x) = 0` for :math:`x_i = 0` for
    :math:`i = 1, ..., n`

    .. [1] Adorio, E. MVF - "Multivariate Test Functions Library in C for
    Unconstrained Global Optimization", 2005

    TODO: the -0.2 factor in the exponent of the first term is given as
    -0.02 in Jamil et al.
    """

    def __init__(self, dimensions=2):
        Benchmark.__init__(self, dimensions)

        self._bounds = list(zip([-35.0] * self.N, [35.0] * self.N))
        self.global_optimum = [[0 for _ in range(self.N)]]
        self.fglob = 0.0
        self.change_dimensionality = True

    def fun(self, x, *args):
        self.nfev += 1
        u = sum(x ** 2)
        v = sum(cos(2 * pi * x))
        return (-20. * exp(-0.2 * sqrt(u / self.N))
                - exp(v / self.N) + 20. + exp(1.))

class Ackley02(Benchmark):

    r"""
    Ackley02 objective function.

    The Ackley02 [1]_ global optimization problem is a multimodal minimization
    problem defined as follows:

    .. math::

        f_{\text{Ackley02}(x) = -200 e^{-0.02 \sqrt{x_1^2 + x_2^2}}

    with :math:`x_i \in [-32, 32]` for :math:`i=1, 2`.

    *Global optimum*: :math:`f(x) = -200` for :math:`x = [0, 0]`

    .. [1] Jamil, M. & Yang, X.-S. A Literature Survey of Benchmark Functions
    For Global Optimization Problems Int. Journal of Mathematical Modelling
    and Numerical Optimisation, 2013, 4, 150-194.

    """
    def __init__(self, dimensions=2):
        Benchmark.__init__(self, dimensions)

        self._bounds = list(zip([-32.0] * self.N, [32.0] * self.N))
        self.global_optimum = [[0 for _ in range(self.N)]]
        self.fglob = -200.

    def fun(self, x, *args):
        self.nfev += 1
        return -200 * exp(-0.02 * sqrt(x[0] ** 2 + x[1] ** 2))

class Ackley03(Benchmark):

    r"""
    Ackley03 [1]_ objective function.

    The Ackley03 global optimization problem is a multimodal minimization
    problem defined as follows:

    .. math::

        f_{\text{Ackley03}}(x) = -200 e^{-0.02 \sqrt{x_1^2 + x_2^2}} +
            5e^{\cos(3x_1) + \sin(3x_2)}

    with :math:`x_i \in [-32, 32]` for :math:`i=1, 2`.

    *Global optimum*: :math:`f(x) = -195.62902825923879` for :math:`x
    = [-0.68255758, -0.36070859]`

    .. [1] Jamil, M. & Yang, X.-S. A Literature Survey of Benchmark Functions
    For Global Optimization Problems Int. Journal of Mathematical Modelling
    and Numerical Optimisation, 2013, 4, 150-194.

     TODO: I think the minus sign is missing in front of the first term in eqn3
      in [1]_.  This changes the global minimum
    """

    def __init__(self, dimensions=2):
        Benchmark.__init__(self, dimensions)

        self._bounds = list(zip([-32.0] * self.N, [32.0] * self.N))
        self.global_optimum = [[-0.68255758, -0.36070859]]
        self.fglob = -195.62902825923879

    def fun(self, x, *args):
        self.nfev += 1
        a = -200 * exp(-0.02 * sqrt(x[0] ** 2 + x[1] ** 2))
        a += 5 * exp(cos(3 * x[0]) + sin(3 * x[1]))
        return a

class Adjiman(Benchmark):

    r"""
    Adjiman objective function.

    The Adjiman [1]_ global optimization problem is a multimodal minimization
    problem defined as follows:

    .. math::

        f_{\text{Adjiman}}(x) = \cos(x_1)\sin(x_2) - \frac{x_1}{(x_2^2 + 1)}

    with, :math:`x_1 \in [-1, 2]` and :math:`x_2 \in [-1, 1]`.

    *Global optimum*: :math:`f(x) = -2.02181` for :math:`x = [2.0, 0.10578]`

    .. [1] Jamil, M. & Yang, X.-S. A Literature Survey of Benchmark Functions
    For Global Optimization Problems Int. Journal of Mathematical Modelling
    and Numerical Optimisation, 2013, 4, 150-194.
    """

    def __init__(self, dimensions=2):
        Benchmark.__init__(self, dimensions)

        self._bounds = ([-1.0, 2.0], [-1.0, 1.0])
        self.global_optimum = [[2.0, 0.10578]]
        self.fglob = -2.02180678

    def fun(self, x, *args):
        self.nfev += 1
        return cos(x[0]) * sin(x[1]) - x[0] / (x[1] ** 2 + 1)

class Alpine01(Benchmark):

    r"""
    Alpine01 objective function.

    The Alpine01 [1]_ global optimization problem is a multimodal minimization
    problem defined as follows:

    .. math::

        f_{\text{Alpine01}}(x) = \sum_{i=1}^{n} \lvert {x_i \sin \left( x_i
        \right) + 0.1 x_i} \rvert

    Here, :math:`n` represents the number of dimensions and :math:`x_i \in
    [-10, 10]` for :math:`i = 1, ..., n`.

    *Global optimum*: :math:`f(x) = 0` for :math:`x_i = 0` for
    :math:`i = 1, ..., n`

    .. [1] Jamil, M. & Yang, X.-S. A Literature Survey of Benchmark Functions
    For Global Optimization Problems Int. Journal of Mathematical Modelling
    and Numerical Optimisation, 2013, 4, 150-194.
    """

    def __init__(self, dimensions=2):
        Benchmark.__init__(self, dimensions)

        self._bounds = list(zip([-10.0] * self.N, [10.0] * self.N))
        self.global_optimum = [[0 for _ in range(self.N)]]
        self.fglob = 0.0
        self.change_dimensionality = True

    def fun(self, x, *args):
        self.nfev += 1

        return sum(abs(x * sin(x) + 0.1 * x))

class Alpine02(Benchmark):

    r"""
    Alpine02 objective function.

    The Alpine02 [1]_ global optimization problem is a multimodal minimization
    problem defined as follows:

    .. math::

        f_{\text{Alpine02}(x) = \prod_{i=1}^{n} \sqrt{x_i} \sin(x_i)

    Here, :math:`n` represents the number of dimensions and :math:`x_i \in [0,
    10]` for :math:`i = 1, ..., n`.

    *Global optimum*: :math:`f(x) = -6.1295` for :math:`x =
    [7.91705268, 4.81584232]` for :math:`i = 1, 2`

    .. [1] Jamil, M. & Yang, X.-S. A Literature Survey of Benchmark Functions
    For Global Optimization Problems Int. Journal of Mathematical Modelling
    and Numerical Optimisation, 2013, 4, 150-194.

    TODO: eqn 7 in [1]_ has the wrong global minimum value.
    """

    def __init__(self, dimensions=2):
        Benchmark.__init__(self, dimensions)

        self._bounds = list(zip([0.0] * self.N, [10.0] * self.N))
        self.global_optimum = [[7.91705268, 4.81584232]]
        self.fglob = -6.12950
        self.change_dimensionality = True

    def fun(self, x, *args):
        self.nfev += 1

        return prod(sqrt(x) * sin(x))

class AMGM(Benchmark):

    r"""
    AMGM objective function.

    The AMGM (Arithmetic Mean - Geometric Mean Equality) global optimization
    problem is a multimodal minimization problem defined as follows

    .. math::

        f_{\text{AMGM}}(x) = \left ( \frac{1}{n} \sum_{i=1}^{n} x_i -
         \sqrt[n]{ \prod_{i=1}^{n} x_i} \right )^2

    Here, :math:`n` represents the number of dimensions and :math:`x_i \in
    [0, 10]` for :math:`i = 1, ..., n`.

    *Global optimum*: :math:`f(x) = 0` for :math:`x_1 = x_2 = ... = x_n` for
    :math:`i = 1, ..., n`

    .. [1] Gavana, A. Global Optimization Benchmarks and AMPGO, retrieved 2015

    TODO: eqn 7 in [1]_ has the wrong global minimum value.
    """

    def __init__(self, dimensions=2):
        Benchmark.__init__(self, dimensions)

        self._bounds = list(zip([0.0] * self.N, [10.0] * self.N))
        self.global_optimum = [[1, 1]]
        self.fglob = 0.0
        self.change_dimensionality = True

    def fun(self, x, *args):
        self.nfev += 1

        f1 = sum(x)
        f2 = prod(x)
        f1 = f1 / self.N
        f2 = f2 ** (1.0 / self.N)
        f = (f1 - f2) ** 2

        return f

max_stars_repo_name: Ennosigaeon/scipy

max_issues_repo_name: Ennosigaeon/scipy

max_forks_repo_name: Ennosigaeon/scipy

reason: from numpy


text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Author: lfurushchev <furushchev@jsk.imi.i.u-tokyo.ac.jp>

import csv
import chainer
import numpy as np
import os.path as osp
from chainercv import utils
from epic_kitchen_utils import get_object_detection_images
from epic_kitchen_utils import parse_object_detection_annotation
from epic_kitchen_object_detection_labels import epic_kitchen_object_detection_label_names

class EpicKitchenObjectDetectionDataset(chainer.dataset.DatasetMixin):
    def __init__(self, data_dir="auto", anno_path="auto", split="train",
                 force_download=False, download_timeout=None,
                 skip_no_image=True, skip_no_bbox=False):
        if split not in ["train"]:
            raise ValueError("Split '%s' not available" % split)

        if data_dir == "auto":
            if anno_path == "auto":
                data_dir, anno_path, annotations = get_object_detection_images(
                    split, force_download=force_download, download_timeout=download_timeout)
            else:
                data_dir = get_object_detection_images(
                    split, force_download=force_download, download_timeout=download_timeout)[0]
                annotations = parse_object_detection_annotation(anno_path)
        else:
            if anno_path == "auto":
                raise ValueError("unknown annotation path")
            else:
                annotations = parse_object_detection_annotation(anno_path)

        self.data_dir = data_dir
        self.split = split
        self.annotations = self._filter_annotations(
            annotations, skip_no_image, skip_no_bbox)

    def __len__(self):
        return len(self.annotations)

    def _filter_annotations(self, annotations, skip_no_image=True, skip_no_bbox=False):
        valid_annos = []
        for anno in annotations:
            if skip_no_image:
                pid, vid = anno["participant_id"], anno["video_id"]
                fn = "%010d.jpg" % (anno["frame"])
                img_path = osp.join(self.data_dir, self.split, pid, vid, fn)
                if not osp.exists(img_path):
                    continue
            if skip_no_bbox:
                annos = anno["annotations"]
                if not annos:
                    continue
            valid_annos.append(anno)
        return valid_annos

    def get_example(self, i):
        dic = self.annotations[i]
        img = self._get_image(dic)
        annos = dic["annotations"]
        if len(annos) == 0:
            bbox = np.zeros((0, 4), dtype=np.float32)
            labels = []
        else:
            bbox = []
            labels = []
            for anno in annos:
                for top, left, height, width in anno["bounding_boxes"]:
                    bbox.append([top, left, top+height, left+width])
                    labels.append(anno["noun_class"])
            bbox = np.asarray(bbox, dtype=np.float32)

        labels = np.asarray(labels, dtype=np.int32)

        return tuple([img, bbox, labels])

    def _get_image(self, anno):
        pid, vid = anno["participant_id"], anno["video_id"]
        fn = "%010d.jpg" % (anno["frame"])
        img_path = osp.join(self.data_dir, self.split, pid, vid, fn)
        return utils.read_image(img_path, dtype=np.float32, color=True)

if __name__ == '__main__':
    from chainercv.visualizations import vis_bbox
    import matplotlib.pyplot as plt

    dataset = EpicKitchenObjectDetectionDataset(split="train")

    print "Loaded dataset. length:", len(dataset)
    print "Starting slideshow..."

    plt.ion()
    fig = plt.figure()
    ax = fig.add_subplot(111)
    plot_init = False
    for i in range(0, len(dataset), 10):
        try:
            img, bbox, label = dataset.get_example(i)
        except Exception as e:
            print i, e
        ax = vis_bbox(img, bbox, label, label_names=epic_kitchen_object_detection_label_names, ax=ax)
        if plot_init is False:
            plt.show()
            plot_init = True
        else:
            fig.canvas.draw()
            fig.canvas.flush_events()
        plt.pause(0.1)
        ax.clear()

max_stars_repo_name: otsubo/chainervr

max_issues_repo_name: otsubo/chainervr

max_forks_repo_name: otsubo/chainervr


text:
import numpy as np
import csv
from torch.backends import cudnn
import logging
import collections

from nn_model.replay_memory import ReplayBuffer
from nn_model.DDPG import DDPG
from nn_model.our_DDPG import DDPG_V2
from nn_model.TD3 import TD3
from nn_model.DQN import DQN
from nn_model.PPO import PPO
from nn_model.AC import AC
from nn_model.DDQN import DDQN
from nn_model.SAC import SAC
from nn_model.RVI import RVI

cudnn.benchmark = True

class Agent:

    def __init__(self, config, env):
        self.config = config
        self.lambda_U = self.config.lambda_U
        self.phi_U = self.config.phi_U
        self.capacity = self.config.capacity
        self.user_id = self.config.userId
        self.redundancy = self.config.redundancy
        self.logger = logging.getLogger("Agent for SOTA solution")

        # define environment
        self.env = env

        kwargs = {
            "state_dim": self.config.state_dim,
            "action_dim": self.config.action_dim,
            "max_action": self.config.max_action,
            "user_id": self.config.userId,
            "discount": self.config.discount,
            "batch_size": self.config.batch_size,
        }

        # Initialize policy
        if self.config.policy == "TD3":
            # Target policy smoothing is scaled wrt the action scale
            kwargs["policy_noise"] = self.config.policy_noise * self.config.max_action
            kwargs["noise_clip"] = self.config.noise_clip * self.config.max_action
            kwargs["policy_freq"] = self.config.policy_freq
            self.policy_model = TD3(**kwargs)
        elif self.config.policy == "OurDDPG":
            self.policy_model = DDPG_V2(**kwargs)
        elif self.config.policy == "DDPG":
            self.policy_model = DDPG(**kwargs)
        elif self.config.policy == "DQN":
            kwargs.pop("max_action")
            kwargs["learn_rate"] = self.config.lrRL
            self.policy_model = DQN(**kwargs)
        elif self.config.policy == "PPO":
            kwargs.pop("max_action")
            kwargs["n_epochs"] = self.config.epoch_rl
            self.policy_model = PPO(**kwargs)
        elif self.config.policy == "SAC":
            kwargs["beta"] = self.config.lrRL
            self.policy_model = SAC(**kwargs)
        elif self.config.policy == "AC":
            self.policy_model = AC(**kwargs)
        elif self.config.policy == "DDQN":
            kwargs.pop("max_action")
            kwargs["alpha"] = self.config.alpha
            kwargs["epsilon"] = self.config.epsilon
            kwargs["mem_size"] = self.config.mem_size
            self.policy_model = DDQN(**kwargs)
        elif self.config.policy == "RVI":
            kwargs["beta"] = self.config.lrRL
            # kwargs.pop("max_action")
            self.policy_model = RVI(**kwargs)

        # state, done = self.env.reset(), False
        self.threshold = 0
        self.out_reward = 200
        self.print_iteration = 1000
        self.episode_num = self.config.max_epoch_rl

        # define memory
        self.memory = ReplayBuffer(self.config.state_dim, self.config.action_dim)

        self.batch_size = self.config.batch_size

        self.cuda = self.config.cuda
        self.reward = 0
        self.record_time = 1

        log_filename = self.config.log_path + self.config.policy + '_' \
                       + str(self.config.videoId) + '_performance_' + str(self.user_id) + '.csv'
        f1 = open(log_filename, 'w', encoding='utf-8')
        self.csv_writer = csv.writer(f1)
        self.csv_writer.writerow(["accuracy", "precision", "recall", "predicted tile"])

        reward_filename = self.config.log_path + self.config.policy + '_reward_' \
                          + str(self.config.videoId) + '_' + str(self.user_id) + '.csv'
        f2 = open(reward_filename, 'w', encoding='utf-8')
        self.csv_writer_reward = csv.writer(f2)
        self.csv_writer_reward.writerow(["reward"])

    '''
    def load_checkpoint(self, index, file_name="checkpoint"):
        filename = self.config.model_path + file_name + str(index) + '.pth.tar'
        try:
            self.logger.info("Loading checkpoint '{}'".format(filename))
            checkpoint = torch.load(filename)

            self.episode_num = checkpoint['episode']
            self.policy_model.load_state_dict(checkpoint['state_dict'])
            # self.optim.load_state_dict(checkpoint['optimizer'])

            self.logger.info("Checkpoint loaded successfully from '{}' at (epoch {}) at (iteration {})\n"
                             .format(self.config.model_path, checkpoint['episode'], checkpoint['iteration']))
        except OSError as e:
            self.logger.info("No checkpoint exists from '{}'. Skipping...".format(self.config.model_path))
            self.logger.info("**First time to train**")

    def save_checkpoint(self, file_name="checkpoint", is_best=0):
        state = {
            'episode': self.episode_num,
            'state_dict': self.policy_model.state_dict(),
            # 'optimizer': self.optim.state_dict(),
        }
        # Save the state
        torch.save(state, self.config.model_path + file_name + str(self.user_id) + '.pth.tar')
        # If it is the best copy it to another file 'model_best.pth.tar'
        if is_best:
            shutil.copyfile(self.config.model_path + file_name,
                            self.config.model_path + 'model_best.pth.tar')
    '''

    def federate_learning(self, models):
        fl_model = models[0]
        worker_state_dict = [x.state_dict() for x in models]
        weight_keys = list(worker_state_dict[0].keys())
        fed_state_dict = collections.OrderedDict()
        for key in weight_keys:
            key_sum = 0
            for i in range(len(models)):
                key_sum = key_sum + worker_state_dict[i][key]
            fed_state_dict[key] = key_sum / len(models)
        # update fed weights to fl model
        fl_model.load_state_dict(fed_state_dict)
        return fl_model

    def load_checkpoint(self, policy_model):
        print('.... federate learning ....')
        if self.config.policy == "SAC":
            model_remote_actor = [policy_model.actor, self.policy_model.actor]
            model_remote_critic_1 = [policy_model.critic_1, self.policy_model.critic_1]
            model_remote_critic_2 = [policy_model.critic_2, self.policy_model.critic_2]
            model_remote_value = [policy_model.value, self.policy_model.value]
            model_remote_target_value = [policy_model.target_value, self.policy_model.target_value]

            self.policy_model = policy_model
            self.policy_model.actor = self.federate_learning(model_remote_actor)
            self.policy_model.critic_1 = self.federate_learning(model_remote_critic_1)
            self.policy_model.critic_2 = self.federate_learning(model_remote_critic_2)
            self.policy_model.value = self.federate_learning(model_remote_value)
            self.policy_model.target_value = self.federate_learning(model_remote_target_value)
        elif self.config.policy == "TD3":
            model_remote_actor = [policy_model.actor, self.policy_model.actor]
            model_remote_actor_target = [policy_model.actor_target, self.policy_model.actor_target]
            model_remote_critic = [policy_model.critic, self.policy_model.critic]
            model_remote_critic_target = [policy_model.critic_target, self.policy_model.critic_target]

            self.policy_model = policy_model
            self.policy_model.actor = self.federate_learning(model_remote_actor)
            self.policy_model.actor_target = self.federate_learning(model_remote_actor_target)
            self.policy_model.critic = self.federate_learning(model_remote_critic)
            self.policy_model.critic_target = self.federate_learning(model_remote_critic_target)
        else:
            self.policy_model = policy_model
        self.policy_model.reset_memory()

    def save_checkpoint(self, is_best=0):
        self.policy_model.save_models()

    def out_print(self, TP, TN, FP, FN, write):
        accuracy = (TP + TN) / (TP + TN + FP + FN)
        if (TP + FP) != 0:
            precision = TP / (TP + FP)
        else:
            precision = 0
        if (TP + FN) != 0:
            recall = TP / (TP + FN)
        else:
            recall = 0

        print("node ID:", str(self.user_id),
              "accuracy:", str(accuracy),
              "precision:", str(precision),
              "recall", str(recall),
              "predicted tile", str(TP + FP))

        if write:
            self.csv_writer.writerow([str(accuracy), str(precision), str(recall), str(TP + FP)])

    def out_print_reward(self, reward, write):

        if write:
            self.csv_writer_reward.writerow([str(reward)])

    def train_one_epoch(self, curr_state, next_state):

        episode_timesteps = 0
        episode_reward = 0
        action_out = []
        best_reward = -1e5
        best_action = np.random.randint(0, 2, 200)

        while True:
            TP, TN, FP, FN = 0, 0, 0, 0
            episode_timesteps += 1

            if self.config.policy == 'PPO':
                action, prob, val = self.policy_model.select_action(np.array(curr_state))
            else:
                action = (self.policy_model.select_action(np.array(curr_state))).clip(-1, 1)

            penalty = 0
            reward = 0
            for index, value in enumerate(action):
                if value > self.threshold:
                    tmp = 1
                else:
                    tmp = 0

                if tmp == next_state[index]:
                    if tmp == 1:
                        TP += 1
                        reward += 2
                        penalty += -2
                    else:
                        TN += 1
                        reward += 1
                        penalty += -1
                else:
                    if tmp == 1:
                        FP += 1
                        reward += -20
                        penalty += 30
                    else:
                        FN += 1
                        reward += -10
                        penalty += 15

            reward = (TP * TP - max(self.phi_U * ((TP + FP) - self.capacity), 0)) - self.lambda_U * penalty

            '''
            if torch.is_tensor(action) and self.cuda:
                action_np = action.cpu().detach().numpy()
            elif type(action) is np.ndarray:
                action_np = action
            else:
                action_np = action.detach().numpy()

            total_num_tiles = int(sum(next_state) * self.redundancy)
            actu_index_list = [k for k, x in enumerate(next_state) if x == 1]
            pred_index_list = np.argsort(action_np)[0:total_num_tiles].tolist()
            right_pred = len(set(actu_index_list) & set(pred_index_list))
            wrong_pred = total_num_tiles - right_pred

            reward = right_pred - 10 * wrong_pred
            '''
            """
            if reward > 0:
                self.record_time += 1
                self.reward = (self.reward + reward) / self.record_time
                reward = reward - self.reward
            """

            if best_reward < reward:
                best_reward = reward
                best_action = action

            # Store data in replay buffer
            if self.config.policy == 'PPO':
                self.policy_model.store_transition(curr_state, action, prob, val, reward, False)
            else:
                self.policy_model.store_transition(curr_state, action, next_state, reward, False)

            if episode_timesteps > self.batch_size:
                self.policy_model.train()

            episode_reward += reward

            self.out_print_reward(reward, True)

            if episode_timesteps % self.print_iteration == 0 and self.print_iteration != 0:
                print(f"Episode T: {episode_timesteps} "
                      f"Slot Reward: {best_reward:.1f} "
                      f"Total Reward: {episode_reward:.1f}")
                self.out_print(TP, TN, FP, FN, False)

            if max(self.episode_num, self.print_iteration) + 1 < episode_timesteps or best_reward > self.out_reward:
                for index, value in enumerate(best_action):
                    if value > self.threshold:
                        action_out.append(1)
                    else:
                        action_out.append(0)
                self.out_print(TP, TN, FP, FN, True)
                self.policy_model.reset_memory()
                break

        return action_out, reward

max_stars_repo_name: uglyghost/FedLive

max_issues_repo_name: uglyghost/FedLive

max_forks_repo_name: uglyghost/FedLive


text:
import os
import uuid
import pandas as pd
import numpy as np
import argparse
p_dir = os.path.dirname(os.path.realpath(__file__)) + "/"
code_dir = p_dir+"scripts/"
def run_matlab_code(seq):

	# step 1: copy matlab code to users folder because matlab scripts can't take parameters
	random_folder = "."+str(uuid.uuid4()).split("-")[-1]
	os.system("mkdir {0};cd {0};ln -s {1}* .".format(random_folder,code_dir))

	# step 2: run matlab
	command = """cd %s;matlab -nodisplay -nodesktop -nosplash -nojvm -r "weka_run(char('%s'),10,0.1,6); exit ; exit()" """%(random_folder,seq)
	print (command)
	os.system(command)

	# step3 parse job.out
	df = pd.read_csv("%s/job.out"%(random_folder))
	print (df.head())

	# step4 delete dir
	os.system("rm -r %s"%(random_folder))

	return df

def ML_features(seq):
	"""Return a list of 8 values for ML table"""
	df = run_matlab_code(seq)
	melting_mean = df['predicted'].mean()
	df['ratio'] = df['predicted']/melting_mean
	df['ratio'] = df['ratio'].apply(lambda x: np.log2(x))
	first_two = df['ratio'].tolist()[:2]
	last_two = df['ratio'].tolist()[-2:]
	ratio_max = df['ratio'].max()
	ratio_min = df['ratio'].min()
	ratio_std = df['ratio'].std()
	return first_two+last_two+[melting_mean,ratio_max,ratio_min,ratio_std]


text:
import os
from os.path import join

import cv2
import numpy as np
from PIL import Image

# Image Edit function(read, convert to Gray, Gap Fill, resize size*size(64), return img)
def image_to_normalized_square(img, size=64) :
    # convert to 1 channel
    if (img.ndim == 3) :
        img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    #record the original dimensions
    width = img.shape[0]
    height = img.shape[1]
    # finding the biggest dimension
    if width >= height:
        length = img.shape[0]
    else :
        length = img.shape[1]
    # create a square matrix, the size of our biggest dimension
    im_result = np.ones((length, length),  dtype = np.uint8) * 255
    wd = int((length - width) / 2)
    hd = int((length - height) / 2)
    # fill the data from both sides, inorder to have the image in the middle & keep the rest unchanged
    im_result[wd:wd + width, hd:hd + height] = img
    # resize to size*size
    im_result = cv2.resize(im_result, (size, size), interpolation=cv2.INTER_AREA)
    return im_result

def mean_block_feature(img):
    block_size = 5
    output_size = 64
    img_result = np.empty((0, output_size))
    for image in img:
        image = extract_mean_block_feature(image, block_size)
        image = image.reshape((1, output_size))
        img_result = np.append(img_result, image, axis=0)
    return img_result

def extract_mean_block_feature(img, block_size):
    if img.shape[0] != img.shape[1]:
        return

    block_num = img.shape[0] // block_size
    img_features = np.empty((0))
    for row in range(block_num):
        for col in range(block_num):
            row_from, row_to = (block_size * row, block_size * (row + 1))
            col_from, col_to = (block_size * col, block_size * (col + 1))
            img_features = np.append(img_features, np.mean(img[row_from:row_to, col_from:col_to]))
    return img_features

def connected_components(img):
    im_gray = cv2.cvtColor(img, cv2.COLOR_RGB2GRAY)
    _, im_bin = cv2.threshold(im_gray, 180, 255, cv2.THRESH_BINARY_INV)
    nlabels, labels, stats, centriods = cv2.connectedComponentsWithStats(im_bin, None, None, None, 8, cv2.CV_32S)
    im_rgb = 255 - cv2.cvtColor(im_bin, cv2.COLOR_GRAY2RGB)
    return nlabels, stats, im_rgb

def get_stat(index, stats):
    left = stats[index, cv2.CC_STAT_LEFT]
    top = stats[index, cv2.CC_STAT_TOP]
    right = stats[index, cv2.CC_STAT_WIDTH] + left
    bottom = stats[index, cv2.CC_STAT_HEIGHT] + top
    return left, top, right, bottom

def img_to_pdf(images_path, pdf_name, pdf_path):
    if not os.path.isdir(pdf_path):
        os.makedirs(pdf_path)

    imgList = []
    for img_path in images_path[1:]:
        img = Image.open(img_path, 'r')
        img = img.convert('RGB')
        imgList.append(img)

    img = Image.open(images_path[0], 'r')
    img = img.convert('RGB')
    img.save(join(pdf_path, pdf_name + ".pdf"), save_all=True, append_images=imgList)

def mergeImg(img_paths, save_filename):
    img1 = cv2.imread(img_paths[0])
    img2 = cv2.imread(img_paths[1])
    y, x, c = img1.shape
    shape = (y, 2*x, c)
    img = np.zeros(shape=shape, dtype=np.uint8)
    img[:, :x, :] = img2
    img[:, x:2*x, :] = img1
    cv2.imwrite(save_filename, img)

if __name__ == "__main__":
    print("Header File...")

max_stars_repo_name: sepehrnaserinia/NeuralNetworks

max_issues_repo_name: sepehrnaserinia/NeuralNetworks

max_forks_repo_name: sepehrnaserinia/NeuralNetworks


text:
#cloneId	cloneCount	cloneFraction	clonalSequence	bestVGene	bestDGene	bestJGene	bestVFamily	bestJFamily	aaSeqCDR3	nSeqCDR3	refPoints
import pytest
import os
import numpy as np
import pandas as pd
from tcrdist.repertoire import TCRrep
from tcrdist import mixcr

time_sensitive = True
@pytest.mark.skipif(time_sensitive, reason="DONT RUN THIS TIME SENSTIVE TEST")
def test_load():
    assert os.path.isfile('tcrdist/test_files/contracting_clones_M_alpha.tsv')
    fn =os.path.join('tcrdist','test_files','contracting_clones_M_alpha.tsv')
    assert os.path.isfile(fn)

@pytest.mark.skipif(time_sensitive, reason="DONT RUN THIS TIME SENSTIVE TEST")
def test_columns():
    fn =os.path.join('tcrdist','test_files','contracting_clones_M_alpha.tsv')
    df = pd.read_csv(fn, sep = "\t")
    assert np.all(df.columns == ['Rank', 'Read.count', 'Read.proportion', 'CDR3.nucleotide.sequence',
       'bestVGene', 'bestDGene', 'bestJGene', 'bestVFamily', 'bestJFamily',
       'CDR3.amino.acid.sequence', 'nSeqCDR3', 'refPoints', 'CD',
       'cdr3aa_length'])

    minervina_alpha_map = {'Rank' : 'clone_id',
                        'Read.count' : 'count',
                        'Read.proportion' : 'prop',
                        'CDR3.nucleotide.sequence': 'cdr3_a_nucseq',
                        'bestVGene' : 'v_b_gene' ,
                        'bestDGene' : 'd_b_gene' ,
                        'bestJGene' : 'j_b_gene' ,
                        'bestVFamily' : 'bestVFamily',
                        'bestJFamily' :'bestVFamily',
                        'CDR3.amino.acid.sequence' : 'cdr3_a_aa',
                        'nSeqCDR3' : 'n_seq_cdr3',
                        'refPoints' : 'ref_points',
                        'CD' : 'cd',
                        'cdr3aa_length' : 'cdr3_a_aa_len'}

    df = df.rename(columns=minervina_alpha_map)
    assert set(df.columns.to_list()) == set(minervina_alpha_map.values())
    dftcr = df[['clone_id','count','v_b_gene','j_b_gene','cdr3_a_aa','cdr3_a_nucseq']]

testfiles = [   ('contracting_clones_M_alpha.tsv', 'alpha'),
                ('contracting_clones_M_beta.tsv',  'beta'),
                ('contracting_clones_W_alpha.tsv', 'alpha'),
                ('contracting_clones_W_beta.tsv',  'beta'),
                ('expanding_clones_M_alpha.tsv',  'alpha'),
                ('expanding_clones_M_beta.tsv',  'beta'),
                ('expanding_clones_W_alpha.tsv',  'alpha'),
                ('expanding_clones_W_beta.tsv', 'beta')]

@pytest.mark.parametrize("f, my_chain", testfiles)
@pytest.mark.skipif(time_sensitive, reason="DONT RUN THIS TIME SENSTIVE TEST")
def test_convert_minervina_to_mixcr_run_tcrdist(f,my_chain):

    fn =os.path.join('tcrdist','test_files', f)
    df = pd.read_csv(fn, sep = "\t")
    df['bestVGene'] = df['bestVGene'].apply(lambda s : s + "*00")
    df['bestJGene'] = df['bestJGene'].apply(lambda s : s + "*00")

    map_minervina_to_mixcr = \
        {'Rank':'cloneId',
        'Read.count':'cloneCount',
        'Read.proportion':'cloneFraction',
        'bestVGene': 'allVHitsWithScore',
        'bestDGene': 'allDHitsWithScore',
        'bestJGene':'allJHitsWithScore',
        'CDR3.nucleotide.sequence':'nSeqCDR3',
        'CDR3.amino.acid.sequence':'aaSeqCDR3',
        'refPoints':'refPoints'}

    df = df.rename(columns = map_minervina_to_mixcr)
    # CREATE A FAUX MIXCR OUTPUT
    df.to_csv('dfmix.clns.txt', index = False, sep = "\t")
    # USE TCRDIST2 TOOL FOR PORTING MIXCR OUTPUTS
    dfmix = mixcr.mixcr_to_tcrdist2(chain = my_chain,
                                    organism = "human",
                                    clones_fn = 'dfmix.clns.txt')

    if my_chain == "alpha":
        assert set(dfmix.columns.to_list()) == set(['clone_id', 'count', 'v_a_gene', 'd_a_gene', 'j_a_gene','cdr3_a_nucseq', 'cdr3_a_aa'])

    assert df.shape[0] == dfmix.shape[0]
    dfmix = mixcr.remove_entries_with_invalid_vgene(dfmix,
                                                    chain = my_chain,
                                                    organism = "human")

    dfmix = mixcr.remove_entries_with_invalid_cdr3(dfmix, chain = my_chain)

    if my_chain == "alpha":
        tr = TCRrep(cell_df = dfmix, organism = "human", chains=['alpha'])
    elif my_chain == "beta":
        tr = TCRrep(cell_df = dfmix, organism = "human", chains=['beta'])

    tr.infer_cdrs_from_v_gene(chain = my_chain, imgt_aligned = True)
    tr.cell_df['subject'] = 'X'
    tr.cell_df['epitope'] = 'X'

    if my_chain == "alpha":
        tr.index_cols = ['clone_id', 'subject', 'epitope',
                        'v_a_gene','j_a_gene',
                        'cdr3_a_aa', 'cdr1_a_aa', 'cdr2_a_aa', 'pmhc_a_aa',
                        'cdr3_a_nucseq']

    elif my_chain == "beta":
        tr.index_cols = ['clone_id', 'subject', 'epitope',
                        'v_b_gene','j_b_gene',
                        'cdr3_b_aa', 'cdr1_b_aa', 'cdr2_b_aa', 'pmhc_b_aa',
                        'cdr3_b_nucseq']

    tr.deduplicate()

    if my_chain == "alpha":
        tr._tcrdist_legacy_method_alpha()
        assert isinstance(tr.cdr3_a_aa_pw, np.ndarray)
        assert isinstance(tr.paired_tcrdist, np.ndarray)
    elif my_chain == "beta":
        tr._tcrdist_legacy_method_beta()
        assert isinstance(tr.cdr3_b_aa_pw, np.ndarray)
        assert isinstance(tr.paired_tcrdist, np.ndarray)

@pytest.mark.skipif(time_sensitive, reason="DONT RUN THIS TIME SENSTIVE TEST")
def test_combine_betas_and_alphas():
    import pytest
    import os
    import numpy as np
    import pandas as pd
    from tcrdist.repertoire import TCRrep
    from tcrdist import mixcr
    import multiprocessing

    testfiles = [('contracting_clones_M_alpha.tsv', 'alpha', 'M', 'contracting'),
                ('contracting_clones_M_beta.tsv',  'beta', 'M', 'contracting'),
                ('contracting_clones_W_alpha.tsv', 'alpha', 'W', 'contracting'),
                ('contracting_clones_W_beta.tsv',  'beta', 'W','contracting'),
                ('expanding_clones_M_alpha.tsv',  'alpha', 'M', 'expanding'),
                ('expanding_clones_M_beta.tsv',  'beta', 'M', 'expanding'),
                ('expanding_clones_W_alpha.tsv',  'alpha', 'W', 'expanding'),
                ('expanding_clones_W_beta.tsv', 'beta', 'W', 'expanding')]

    betas = []
    alphas = []

    for f,my_chain,sub,group in testfiles:
        fn =os.path.join('tcrdist','test_files', f)
        df = pd.read_csv(fn, sep = "\t")
        df['bestVGene'] = df['bestVGene'].apply(lambda s : s + "*00")
        df['bestJGene'] = df['bestJGene'].apply(lambda s : s + "*00")
        print(df.columns)
        map_minervina_to_mixcr = \
            {'Rank':'cloneId',
            'Read.count':'cloneCount',
            'Read.proportion':'cloneFraction',
            'bestVGene': 'allVHitsWithScore',
            'bestDGene': 'allDHitsWithScore',
            'bestJGene':'allJHitsWithScore',
            'CDR3.nucleotide.sequence':'nSeqCDR3',
            'CDR3.amino.acid.sequence':'aaSeqCDR3',
            'refPoints':'refPoints'}

        df = df.rename(columns = map_minervina_to_mixcr)
        # CREATE A FAUX MIXCR OUTPUT
        print(df.columns)
        df.to_csv('dfmix.clns.txt', index = False, sep = "\t")
        dfmix = mixcr.mixcr_to_tcrdist2(chain = my_chain,
                                organism = "human",
                                clones_fn = 'dfmix.clns.txt')
        dfmix['CD'] = df['CD'].copy()
        dfmix['proportion'] = df['cloneFraction'].copy()
        dfmix = mixcr.remove_entries_with_invalid_vgene(dfmix,
                                                chain = my_chain,
                                                organism = "human")

        dfmix = mixcr.remove_entries_with_invalid_cdr3(dfmix, chain = my_chain)
        dfmix['source'] = f
        dfmix['subject'] = sub
        dfmix['epitope'] = 'X'
        dfmix['trajectory'] = group

        if my_chain == "alpha":
            alphas.append(dfmix)
        elif my_chain == "beta":
            betas.append(dfmix)

    betas_joined = pd.concat(betas)
    alpha_joined = pd.concat(alphas)

    testsets = [(alpha_joined, 'alpha'),
                (betas_joined, 'beta')]

    tcr_rep_results = dict()
    for dfmix, my_chain in testsets:
        if my_chain == "alpha":
            tr = TCRrep(cell_df = dfmix, organism = "human", chains=['alpha'])
        elif my_chain == "beta":
            tr = TCRrep(cell_df = dfmix, organism = "human", chains=['beta'])

        tr.infer_cdrs_from_v_gene(chain = my_chain, imgt_aligned = True)

        if my_chain == "alpha":
            tr.index_cols = ['clone_id', 'subject', 'epitope', 'trajectory',
                            'v_a_gene','j_a_gene',
                            'cdr3_a_aa', 'cdr1_a_aa', 'cdr2_a_aa', 'pmhc_a_aa',
                            'cdr3_a_nucseq', 'CD','source']

        elif my_chain == "beta":
            tr.index_cols = ['clone_id', 'subject', 'epitope', 'trajectory',
                            'v_b_gene','j_b_gene',
                            'cdr3_b_aa', 'cdr1_b_aa', 'cdr2_b_aa', 'pmhc_b_aa',
                            'cdr3_b_nucseq', 'CD','source']

        tr.deduplicate()
        my_processes = 10
        if my_processes > multiprocessing.cpu_count():
             my_processes = multiprocessing.cpu_count()

        if my_chain == "alpha":
            tr._tcrdist_legacy_method_alpha(processes = my_processes)
            assert isinstance(tr.cdr3_a_aa_pw, np.ndarray)
            assert isinstance(tr.paired_tcrdist, np.ndarray)
        elif my_chain == "beta":
            tr._tcrdist_legacy_method_beta(processes = my_processes)
            assert isinstance(tr.cdr3_b_aa_pw, np.ndarray)
            assert isinstance(tr.paired_tcrdist, np.ndarray)

        tcr_rep_results[my_chain] = tr
    return tcr_rep_results
        #import pickle
        #with open("minervina_dynamic_clones.p", 'wb') as p:
        #     pickle.dump(tcr_rep_results, p)


text:
# -*- coding: utf-8 -*-
"""utils.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1Rnkq2VXUAyNv_FW7tnuIkG61r0jCuCc7
"""

import torch
import random
import numpy as np
import os
import torchvision
import config
from torchvision.utils import save_image

def plot_to_tensorboard(writer, loss_critic, loss_gen, real, fake, tensorboard_step):
  writer.add_scalar("Loss Critic", loss_critic, global_step=tensorboard_step)

  with torch.no_grad():
    img_grid_real = torchvision.utils.make_grid(real[:8], normalize=True)
    img_grid_fake = torchvision.utils.make_grid(fake[:8], normalize=True)
    writer.add_image("Real", img_grid_real, global_step=tensorboard_step)
    writer.add_image("Fake", img_grid_fake, global_step=tensorboard_step)

def gradient_penalty(critic, real, fake, alpha, train_step, device="cpu"):
  BATCH_SIZE, C, H, W = real.shape
  beta = torch.randn((BATCH_SIZE, 1, 1, 1)).repeat(1, C, H, W).to(device)

  interpolated_images = real * beta + fake.detach() * (1 - beta)
  interpolated_images.requires_grad_(True)

  mixed_scores = critic(interpolated_images, alpha, train_step)

  gradient = torch.autograd.grad(
      inputs=interpolated_images,
      outputs=mixed_scores,
      grad_outputs=torch.ones_like(mixed_scores),
      create_graph=True,
      retain_graph=True,
  )[0]
  gradient = gradient.view(gradient.shape[0], -1)
  gradient_norm = gradient.norm(2, dim=1)
  gradient_penalty = torch.mean((gradient_norm - 1) ** 2)
  return gradient_penalty

def save_checkpoint(model, optimizer, filename="StyleGAN_checkpoint.pth.tar"):
  print("=> Saving checkpoint")
  checkpoint = {
      "state_dict": model.state_dict(),
      "optimizer": optimizer.state_dict(),
  }
  torch.save(checkpoint, filename)

def load_checkpoint(checkpoint_file, model, optimizer=None, lr=None):
  print("=> Loading checkpoint")
  checkpoint = torch.load(checkpoint_file, map_location="cpu")
  model.load_state_dict(checkpoint["state_dict"])
  if optimizer:
    optimizer.load_state_dict(checkpoint["optimizer"])
    for param_group in optimizer.param_groups:
      param_group["lr"] = lr

def generate_examples(gen, steps, epoch, n=100):
  gen.eval()
  alpha = 0.3
  for i in range(n):
    with torch.no_grad():
      noise = torch.zeros(1, config.L_DIM).to(config.DEVICE)
      img = gen(noise, steps, alpha)
      save_image(img * 0.5 * 0.5, f"/content/drive/MyDrive/BTP/savedExamples/saved_examples_{epoch+16}_{i}.png")

  gen.train()