text:
"""
This script gather several functions which are related to the profile
of galaxy clusters.
"""

import numpy as np
from scipy.special import gamma
import scipy.integrate as integrate
import scipy.interpolate as interpolate
import matplotlib.pyplot as plt

from ClusterModel.ClusterTools import map_tools

#===================================================
#========== Beta model
#===================================================
def beta_model(r3d_kpc, n0, r_c, beta):
    """
    Compute a beta model

    Parameters
    ----------
    - r3d_kpc: array of radius in kpc
    - r_c : core radius parameter
    - n_0  : normalization
    - beta : slope of the profile

    Outputs
    --------
    - beta model profile as a function of the input radius vector

    """

    return n0 * (1 + (r3d_kpc / r_c)**2)**(-3.0*beta/2.0)

#===================================================
#========== Beta model
#===================================================
def beta_model_derivative(r3d_kpc, n0, r_c, beta):
    """
    Compute the derivative of a beta model dn/dr

    Parameters
    ----------
    - r3d_kpc: array of radius in kpc
    - r_c : core radius parameter
    - n_0  : normalization
    - beta : slope of the profile

    Outputs
    --------
    - beta model derivative profile as a function of the input radius vector

    """

    return -3.0*n0*beta*r3d_kpc*(1+(r3d_kpc/r_c)**2)**(-3.0*beta/2.0-1.0)/r_c**2

#===================================================
#========== Simplified Vikhlinin model
#===================================================
def svm_model(r3d_kpc, n0, r_c, beta, r_s=1.0, gamma=3.0, epsilon=0.0, alpha=0.0):
    """
    Compute a Simplified Vikhlinin model

    Parameters
    ----------
    - r3d_kpc: array of radius in kpc
    - r_c : core radius parameter
    - n_0  : normalization
    - beta : slope of the profile
    - r_s  : characteristic radius parameter (kpc)
    - gamma : slope of the profile
    - epsilon : slope of the profile

    Outputs
    --------
    - SVM model profile as a function of the input radius vector

    """

    term1 = n0 * (1 + (r3d_kpc / r_c)**2)**(-3.0*beta/2.0)
    term2 = (r3d_kpc / r_c)**(-alpha/2.0)
    term3 = (1 + (r3d_kpc / r_s)**gamma)**(-epsilon/2.0/gamma)

    return term1 * term2 * term3

#===================================================
#========== Simplified Vikhlinin model
#===================================================
def svm_model_derivative(r3d_kpc, n0, r_c, beta, r_s=1.0, gamma=3.0, epsilon=0.0, alpha=0.0):
    """
    Compute the derivative of the Simplified Vikhlinin model

    Parameters
    ----------
    - r3d_kpc: array of radius in kpc
    - r_c : core radius parameter
    - n_0  : normalization
    - beta : slope of the profile
    - r_s  : characteristic radius parameter (kpc)
    - gamma : slope of the profile
    - epsilon : slope of the profile

    Outputs
    --------
    - SVM derivative model profile as a function of the input radius vector

    """

    t1 = n0 / (r3d_kpc*(r_c**2 + r3d_kpc**2))
    t2 = (r3d_kpc / r_c)**(-alpha/2.0)
    t3 = (1 + (r3d_kpc / r_c)**2)**(-3.0*beta/2.0)
    t4 = (1 + (r3d_kpc/r_s)**gamma)**(-epsilon/2.0/gamma - 1.0)
    t51 = -alpha * (r_c**2 + r3d_kpc**2) * ((r3d_kpc/r_s)**gamma + 1) / 2.0
    t52 = -3*beta*r3d_kpc**2*((r3d_kpc/r_s)**gamma + 1)
    t53 = -epsilon*(r_c**2 + r3d_kpc**2) / 2.0 * (r3d_kpc/r_s)**gamma

    return t1*t2*t3*t4*(t51+t52+t53)

#===================================================
#========== NFW model
#===================================================
def NFW_model(r3d_kpc, norm, rs):
    """
    Compute a NFW model

    Parameters
    ----------
    - r3d_kpc (kpc): array of radius
    - norm : the normalization
    - rs (kpc): characteristic radius parameter

    Outputs
    --------
    - NFW model profile as a function of the input radius vector

    """

    return norm / (r3d_kpc / rs) / (1 + (r3d_kpc / rs))**2

#===================================================
#========== gNFW model
#===================================================
def gNFW_model(r3d_kpc, P0, r_p, slope_a=1.33, slope_b=4.13, slope_c=0.31):
    """
    Compute a gNFW model

    Parameters
    ----------
    - r3d_kpc (kpc): array of radius
    - P0 : normalization
    - r_p (kpc): characteristic radius parameter
    - sope_a : intermediate slope parameter
    - sope_b : outer slope parameter
    - sope_c : inner slope parameter

    Outputs
    --------
    - gNFW model profile as a function of the input radius vector

    """

    return P0 * (r3d_kpc / r_p)**(-slope_c) * (1 + (r3d_kpc / r_p)**(slope_a))**((slope_c-slope_b)/slope_a)

#===================================================
#========== gNFW model
#===================================================
def gNFW_model_derivative(r3d_kpc, P0, r_p, slope_a=1.33, slope_b=4.13, slope_c=0.31):
    """
    Compute the derivative of the gNFW model (dP/dr)

    Parameters
    ----------
    - r3d_kpc (kpc): array of radius
    - P0 : normalization
    - r_p (kpc): characteristic radius parameter
    - sope_a : intermediate slope parameter
    - sope_b : outer slope parameter
    - sope_c : inner slope parameter

    Outputs
    --------
    - gNFW derivative model profile as a function of the input radius vector

    """

    t1 =  -P0 * (r3d_kpc / r_p)**(-slope_c)
    t2 = (1 + (r3d_kpc / r_p)**(slope_a))**((slope_c-slope_b-slope_a)/slope_a)
    t3 = (slope_b*(r3d_kpc/r_p)**slope_a + slope_c) / r3d_kpc

    return t1 * t2 * t3

#===================================================
#========== Project beta model
#===================================================
def proj_beta_model(r2d_kpc, n0, r_c, beta):
    """
    Compute a projected beta model:
    P(R) = \int n_e dl at given R

    Parameters
    ----------
    - r2d_kpc: array of projected radius at which to compute integration
    - n0  : normalization
    - r_c : core radius parameter
    - beta : slope of the profile

    Outputs
    --------
    - The projected profile in units of kpc times original profile

    """

    return np.sqrt(np.pi) * n0 * r_c * gamma(1.5*beta - 0.5) / gamma(1.5*beta) * (1 + (r2d_kpc/r_c)**2)**(0.5-1.5*beta)

#===================================================
#========== Integrated any model
#===================================================
def proj_any_model_naive(r3d, p3d, Rpmax=6000.0, Rpstep=10.0):
    """
    Compute the projection of any given profile in a
    naive way (simple steps and trapezoidal integration).
    P(R) = \int p(r) dl for values of R

    Parameters
    ----------
    - r3d (kpc): input radius array
    - p3d (arb): input profile array
    - Rpmax (kpc): projected range up to which we compute the profile
    - Rpstep (kpc): size of the steps along projection direction

    Outputs
    --------
    - r2d (kpc): output projected radius array
    - p2d (arb kpc): output projected profile

    """

    rproj_Nstep = int(Rpmax/Rpstep)
    r2d = np.linspace(Rpstep, Rpmax, num=rproj_Nstep)
    p2d = np.empty(rproj_Nstep)

    for istep in range(0,rproj_Nstep):
        w = r3d > r2d[istep]
        r3d_w = r3d[w]
        p3d_w = p3d[w]
        integ = 2.0 * p3d_w * r3d_w / np.sqrt((r3d_w**2.0 - r2d[istep]**2.0))
        p2d[istep] = np.trapz(integ, x=r3d_w)

    return r2d, p2d

#===================================================
#========== Integrated any model
#===================================================
def proj_any_model(r3d, p3d, Npt=100, Rmax=5000.0, Rpmax=6000.0, Rpstep=10.0, Rp_input=None):
    """
    Compute the projection of any given profile. This is done by
    creating a lin-log mapping of the distance from the center and summing
    the pixels directly. It works very well in terms of speed and is very accurate.
    P(R) = \int p(r) dl for values of R

    Parameters
    ----------
    - r3d (kpc): input radius array
    - p3d (arb): input profile array
    - Npt : number of point in the integration line of sight
    - Rmax (kpc): limit radius for integration
    - Rpmax (kpc): projected range up to which we compute the profile
    - Rpstep (kpc): size of the steps along projection direction
    - Rp_input (kpc): input vector for the projected radius

    Outputs
    --------
    - r2d (kpc): output projected radius array
    - p2d (arb kpc): output projected profile

    """

    #---------- Line-of-Sight grid
    los_bin = np.logspace(0, np.log10(Rmax), num=Npt)
    los_bin = np.concatenate((np.array([0]), los_bin), axis=None)

    los_step = los_bin - np.roll(los_bin, 1)              # Steps
    los_cbin = (los_bin + np.roll(los_bin, 1))/2.0        # Center of bins
    los_step = los_step[1:]                               # remove first bin which does not exist
    los_cbin = los_cbin[1:]                               # remove first bin which does not exist

    #---------- Projected radius grid
    if Rp_input is not None :
        Rproj = Rp_input
        rproj_Nstep = len(Rproj)
    else:
        rproj_Nstep = int(Rpmax/Rpstep)                       # Number of bins along L.o.S
        Rproj = np.linspace(Rpstep, Rpmax, num=rproj_Nstep)   # Vector going up to rproj_out with rp_step steps

    #---------- Compute the radius grid
    los_cbin_map = np.tile(los_cbin, [rproj_Nstep,1])     # Line of sight bin central value map
    los_step_map = np.tile(los_step, [rproj_Nstep,1])     # Line of sight bin size map
    Rproj_map    = np.tile(Rproj, [Npt,1]).T              # Projected radius map
    r_map = np.sqrt(Rproj_map**2 + los_cbin_map**2)       # Physical radius from the center map

    #---------- Interpolate the profile to the grid
    P_Rl = map_tools.profile2map(p3d, r3d, r_map)

    #---------- Integration
    integrant = P_Rl * los_step_map
    lim_out = (r_map > Rmax)
    integrant[lim_out] = 0

    Pproj = 2*np.sum(integrant, axis=1)

    return Rproj, Pproj

#===================================================
#========== Volume integrate of beta model
#===================================================
def get_volume_beta_model(Rmax_kpc, n0, r_c, beta):
    """
    Compute the volume integrated beta model

    Parameters
    ----------
    - Rmax_kpc (kpc): radius at which to stop integration
    - n0  : normalization
    - r_c : core radius parameter
    - beta : slope of the profile

    Outputs
    --------
    - The integrated profile in units of kpc^3 times the original profile

    """

    # Defines a function which returns the integrand
    def get_volume_beta_model_integrand(r3d_kpc, n0, r_c, beta):
        return 4.0*np.pi * r3d_kpc**2 * n0 * (1 + (r3d_kpc / r_c)**2)**(-3*beta/2.0)

    # Integrate the integrand
    output = integrate.quad(get_volume_beta_model_integrand, 0.0, Rmax_kpc, args=(n0, r_c, beta))[0]

    return output

#===================================================
#========== Volume any profile
#===================================================
def get_volume_any_model(r3d, p3d, Rmax_kpc, Npt=1000):
    """
    Compute the volume integrated profile:
    \int_0^Rmax 4pi r^2 p(r) dr

    Parameters
    ----------
    - Rmax_kpc (kpc): radius at which to stop integration
    - r3d (kpc) : input radius in kpc
    - p3d : input profile

    Outputs
    --------
    - The integrated profile in units of kpc^3 times the original profile

    Note
    --------
    The error can be large if the input profile is not well sampled. Use ~1000
    points in r3d/p3d for safety, and well sampled until Rmax. See also
    define_safe_radius_array() for this.

    """

    if Rmax_kpc > 1.001*np.amax(r3d):
        print('!!!!!!!!!! WARNING: you try to integrate beyond function limits !!!!!!!!!!')

    # interpolate to the location of interest
    r3d_i = np.logspace(np.log10(np.amin(r3d)), np.log10(Rmax_kpc), num=Npt)
    r3d_i = np.concatenate((np.array([0]), r3d_i), axis=None)

    itpl = interpolate.interp1d(r3d, p3d, kind='cubic', fill_value='extrapolate')
    p3d_i = itpl(r3d_i)

    # integrate
    output = np.trapz(4.0*np.pi * r3d_i**2 * p3d_i, r3d_i)

    return output

#===================================================
#========== Surface any profile
#===================================================
def get_surface_any_model(r2d, p2d, Rmax_kpc, Npt=1000):
    """
    Compute the surface integrated profile:
    \int_0^Rmax 2pi r p(r) dr

    Parameters
    ----------
    - Rmax_kpc (kpc): radius at which to stop integration
    - r2d (kpc) : input radius in kpc
    - p2d : input profile

    Outputs
    --------
    - The integrated profile in units of kpc^2 times the original profile

    Note
    --------

    """

    if Rmax_kpc > 1.001*np.amax(r2d):
        print('!!!!!!!!!! WARNING: you try to integrate beyond function limits !!!!!!!!!!')

    # interpolate to the location of interest
    r2d_i = np.logspace(np.log10(np.amin(r2d)), np.log10(Rmax_kpc), num=Npt)
    r2d_i = np.concatenate((np.array([0]), r2d_i), axis=None)

    itpl = interpolate.interp1d(r2d, p2d, kind='cubic', fill_value='extrapolate')
    p2d_i = itpl(r2d_i)

    # integrate
    output = np.trapz(2.0*np.pi * r2d_i * p2d_i, r2d_i)

    return output

#===================================================
#========== Define radius vector
#===================================================
def define_safe_radius_array(r_in, Rmin=1.0, Rmax=None, Nptmin=1000):
    """
    Defines a radius that safely sample the profiles for integration, based
    on the radius at which we want to get the samples

    Parameters
    ----------
    - r_in (kpc): input sampling radius
    - Rmin (kpc) : minimal for the sampling
    - Nptmin : minimal number of points in the sampling

    Outputs
    --------
    - r_out (kpc): radius used to sample the profile we want to integrate

    """

    # array case
    if len(r_in) > 1:
        logmax = np.amax(np.log10(r_in))
        logmin = np.amin(np.log10(r_in))
        Npt_per_decade = len(r_in) / (logmax - logmin)

        # We want at least to cover one decade
        if (logmax - logmin) < 1:
            logmin = logmax - 1

        # maximum radius remains the same, minimum radius should go to 1 kpc at least
        if Rmax is not None :
            if logmax < np.log10(Rmax):
                logmax_value = np.log10(Rmax)
            else:
                logmax_value = logmax
        else:
            logmax_value = logmax

        if logmin > np.log10(Rmin):
            logmin_value = np.log10(Rmin)
        else:
            logmin_value = logmin

        # keep the same number of point per decade, with at leat Nptmin points
        Npt = int(Npt_per_decade * (logmax_value - logmin_value))
        if Npt < Nptmin: Npt = Nptmin

    # case a scalar is given
    else:
        r_in = float(r_in)
        Npt = Nptmin

        if Rmax is not None :
            if r_in < Rmax:
                logmax_value = np.log10(Rmax)
            else:
                logmax_value = np.log10(r_in)
        else:
            logmax_value = np.log10(r_in)

        if r_in/10.0 > Rmin:
            logmin_value = np.amin(np.log10(Rmin))
        else:
            logmin_value = np.log10(r_in/10.0)

    # Defines the radius
    r_out = np.logspace(logmin_value, logmax_value, Npt)

    return r_out

#===================================================
#========== Main function
#===================================================
if __name__ == "__main__":

    print('========== Cluster profile module ===========')
    print('This gather modules dedicated to compute     ')
    print('profile related quantities.                  ')
    print('                                             ')
    print('Here we present some example and checks.     ')
    print('=============================================')

    # Cluster properties
    r3d = np.logspace(0,5, 100)    # radius (kpc)
    r_c = 200.0
    beta = 0.75
    n0 = 1.0

    # Get profiles in different ways
    n_r = beta_model(r3d, n0, r_c, beta)

    R2d, n_R1 = proj_any_model(r3d, n_r, Npt=1000, Rmax=10000.0, Rpmax=5000.0, Rpstep=10.0)
    n_R2      = proj_beta_model(R2d, n0, r_c, beta)

    vol1 = np.zeros(len(r3d))
    vol2 = np.zeros(len(r3d))
    for i in range(len(r3d)):
        vol1[i] = get_volume_beta_model(r3d[i], n0, r_c, beta)
        vol2[i] = get_volume_any_model(r3d, n_r, r3d[i], Npt=10000)

    # Plot the projected profile
    fig1 = plt.figure(1)
    plt.loglog(R2d, n_R1, label='Numerical projection')
    plt.loglog(R2d, n_R2, label='Analytical projection')
    plt.xlabel('Projected radius (kpc)')
    plt.ylabel('Projected profile')
    plt.title('Projected profile')
    plt.legend(loc='upper right')

    # Plot the relative error
    fig2 = plt.figure(2)
    plt.loglog(R2d, np.abs((n_R1 - n_R2)/n_R2) * 100)
    plt.xlabel('Projected radius (kpc)')
    plt.ylabel('Relative error (%)')
    plt.title('Error on projected profile')

    # Plot the volume integral
    fig3 = plt.figure(3)
    plt.loglog(r3d, vol1, label='beta model integration')
    plt.loglog(r3d, vol2, label='general integration')
    plt.xlabel('radius (kpc)')
    plt.ylabel('Volume (<r)')
    plt.title('Volume integration')

    # Plot the relative error
    fig3 = plt.figure(4)
    plt.loglog(r3d, np.abs((vol1-vol2)/vol1)*100)
    plt.xlabel('radius (kpc)')
    plt.ylabel('Relative error (%)')
    plt.title('Volume integration error')

    plt.show()

max_stars_repo_name: deimqs/ClusterModel

max_issues_repo_name: deimqs/ClusterModel

max_forks_repo_name: deimqs/ClusterModel

reason: import numpy,import scipy,from scipy


text:
# Copyright (c) 2020 PaddlePaddle Authors. All Rights Reserved.
#
# 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.

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

import os
import glob
import numpy as np
import argparse

from PIL import Image
from scipy.misc import imsave

import paddle.fluid as fluid
from check import check_gpu, check_version

from model import Model, Input, set_device
from cyclegan import Generator, GeneratorCombine

def main():
    place = set_device(FLAGS.device)
    fluid.enable_dygraph(place) if FLAGS.dynamic else None

    # Generators
    g_AB = Generator()
    g_BA = Generator()
    g = GeneratorCombine(g_AB, g_BA, is_train=False)

    im_shape = [-1, 3, 256, 256]
    input_A = Input(im_shape, 'float32', 'input_A')
    input_B = Input(im_shape, 'float32', 'input_B')
    g.prepare(inputs=[input_A, input_B])
    g.load(FLAGS.init_model, skip_mismatch=True, reset_optimizer=True)

    out_path = FLAGS.output + "/single"
    if not os.path.exists(out_path):
        os.makedirs(out_path)
    for f in glob.glob(FLAGS.input):
        image_name = os.path.basename(f)
        image = Image.open(f).convert('RGB')
        image = image.resize((256, 256), Image.BICUBIC)
        image = np.array(image) / 127.5 - 1

        image = image[:, :, 0:3].astype("float32")
        data = image.transpose([2, 0, 1])[np.newaxis, :]

        if FLAGS.input_style == "A":
            _, fake, _, _ = g.test([data, data])

        if FLAGS.input_style == "B":
            fake, _, _, _ = g.test([data, data])

        fake = np.squeeze(fake[0]).transpose([1, 2, 0])

        opath = "{}/fake{}{}".format(out_path, FLAGS.input_style, image_name)
        imsave(opath, ((fake + 1) * 127.5).astype(np.uint8))
        print("transfer {} to {}".format(f, opath))

if __name__ == "__main__":
    parser = argparse.ArgumentParser("CycleGAN inference")
    parser.add_argument(
        "-d", "--dynamic", action='store_false', help="Enable dygraph mode")
    parser.add_argument(
        "-p",
        "--device",
        type=str,
        default='gpu',
        help="device to use, gpu or cpu")
    parser.add_argument(
        "-i",
        "--input",
        type=str,
        default='./image/testA/123_A.jpg',
        help="input image")
    parser.add_argument(
        "-o",
        '--output',
        type=str,
        default='output',
        help="The test result to be saved to.")
    parser.add_argument(
        "-m",
        "--init_model",
        type=str,
        default='checkpoint/199',
        help="The init model file of directory.")
    parser.add_argument(
        "-s", "--input_style", type=str, default='A', help="A or B")
    FLAGS = parser.parse_args()
    print(FLAGS)
    check_gpu(str.lower(FLAGS.device) == 'gpu')
    check_version()
    main()

ext: py

max_stars_repo_name: LielinJiang/hapi

max_issues_repo_name: LielinJiang/hapi

max_forks_repo_name: LielinJiang/hapi

reason: import numpy,from scipy


text:
"""
Outdated module
"""

import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle
import numpy as np
from pyhdx.support import autowrap
import warnings

def plot_residue_map(pm, scores=None, ax=None, cmap='jet', bad='k', cbar=True, **kwargs): # pragma: no cover
    """
    FUNCTION IS MOST LIKELY OUT OF DATE

    Parameters
    ----------
    pm
    scores
    ax
    cmap
    bad
    cbar
    kwargs

    Returns
    -------

    """

    warnings.warn("This function will be removed", DeprecationWarning)

    img = (pm.X > 0).astype(float)
    if scores is not None:
        img *= scores[:, np.newaxis]
    elif pm.rfu is not None:
        img *= pm.rfu[:, np.newaxis]

    ma = np.ma.masked_where(img == 0, img)
    cmap = mpl.cm.get_cmap(cmap)
    cmap.set_bad(color=bad)

    ax = plt.gca() if ax is None else ax
    ax.set_facecolor(bad)

    im = ax.imshow(ma, cmap=cmap, **kwargs)
    if cbar:
        cbar = plt.colorbar(im, ax=ax)
        cbar.set_label('Uptake (%)')

    ax.set_xlabel('Residue number')
    ax.set_ylabel('Peptide index')

def make_kinetics_figure(pm_dict, cmap='cool'): # pragma: no cover
    """
    FUNCTION IS MOST LIKELY OUT OF DATE

    :param pm_dict: dictionary of PeptideMeasuements
    :param times: array_like of
    :param cmap: optional string indicating which colormap to use
    :return:
    """

    """returns matplotlib figure for visualization of kinetics"""

    warnings.warn("This function will be removed", DeprecationWarning)

    fig, (ax1, ax2) = plt.subplots(2, figsize=(8, 4), sharex=True, gridspec_kw={'hspace': 0})
    times = [p.exposure for p in pm_dict.values()]
    s = pm_dict[next(iter(pm_dict))]

    norm = mpl.colors.Normalize(vmin=0, vmax=np.max(times))
    normed_times = norm(times)
    r_number = np.arange(s.start, s.stop + 1)

    colors = mpl.cm.get_cmap(cmap, len(pm_dict))(normed_times)

    for c, (k, v) in zip(colors, pm_dict.items()):
        ax1.plot(r_number, v.scores_average, color=c, marker='.', linestyle='')

    ax1.set_ylabel('Score (%)')

    ax2.plot(r_number, np.ones_like(r_number), marker='.', linestyle='', color='k')

    # lc = LineCollection(self._lc_data, colors='r')
    # ax2.add_collection(lc)
    #        ax2.axhline(self.rate_max, color='r', autolim=False)
    ax2.set_yscale('log')
    ax2.set_xlabel('Residue number')
    ax2.set_ylabel('Rate constant\n (min$^{-1})$')

    fig.align_ylabels()  # todo doesnt seem to work

    fig.subplots_adjust(right=0.85)

    cbar_ax = fig.add_axes([0.87, 0.05, 0.02, 0.9])

    cb1 = mpl.colorbar.ColorbarBase(cbar_ax, cmap=mpl.cm.get_cmap(cmap),
                                    norm=norm,
                                    orientation='vertical')
    cb1.set_label('Time (min)')
    #   plt.tight_layout()

    return fig, (ax1, ax2, cbar_ax)

def plot_peptides(pm, ax, wrap=None, color=True, labels=False, cbar=False, intervals='corrected', cmap='jet', **kwargs):
    """

    TODO: needs to be checked if intervals (start, end) are still accurately taking inclusive, exclusive into account
    Plots peptides as rectangles in the provided axes

    Parameters
    ----------
    pm
    wrap
    ax
    color
    labels
    cmap
    kwargs

    Returns
    -------

    """

    wrap = wrap or autowrap(pm.data['start'], pm.data['end'])
    rect_kwargs = {'linewidth': 1, 'linestyle': '-', 'edgecolor': 'k'}
    rect_kwargs.update(kwargs)

    cmap = mpl.cm.get_cmap(cmap)
    norm = mpl.colors.Normalize(vmin=0, vmax=1)
    i = -1

    for p_num, e in enumerate(pm.data):
        if i < -wrap:
            i = -1

        if color:
            c = cmap(norm(e['rfu']))
        else:
            c = '#707070'

        if intervals == 'corrected':
            start, end = 'start', 'end'
        elif intervals == 'original':
            start, end = '_start', '_end'
        else:
            raise ValueError(f"Invalid value '{intervals}' for keyword 'intervals', options are 'corrected' or 'original'")

        width = e[end] - e[start]
        rect = Rectangle((e[start] - 0.5, i), width, 1, facecolor=c, **rect_kwargs)
        ax.add_patch(rect)
        if labels:
            rx, ry = rect.get_xy()
            cy = ry
            cx = rx
            ax.annotate(str(p_num), (cx, cy), color='k', fontsize=6, va='bottom', ha='right')

        i -= 1

    if cbar:
        scalar_mappable = mpl.cm.ScalarMappable(norm=norm, cmap=cmap)
        plt.colorbar(scalar_mappable, label='Percentage D')

    ax.set_ylim(-wrap, 0)
    end = pm.interval[1]
    ax.set_xlim(0, end)
    ax.set_yticks([])

max_stars_repo_name: sajetan/PyHDX

max_issues_repo_name: sajetan/PyHDX

max_forks_repo_name: sajetan/PyHDX

reason: import numpy


text:
import numpy

def NoConstraint(lower, upper):
    """No constraint."""

    def cons(U):
        return U

    return cons

def Random(lower, upper):
    """
    Random constraint.

    Solutions that are in the infeasible space are resampled.

    """

    def cons(U):
        return numpy.where(
            numpy.logical_or(U < lower, U > upper),
            numpy.random.uniform(lower, upper, U.shape),
            U,
        )

    return cons

_constraints_map = {
    None: NoConstraint,
    "Random": Random,
}

max_stars_repo_name: keurfonluu/stochopy

max_issues_repo_name: keurfonluu/stochopy

max_forks_repo_name: keurfonluu/stochopy


text:
# coding=utf-8

#     National Oceanic and Atmospheric Administration (NOAA)
#     Alaskan Fisheries Science Center (AFSC)
#     Resource Assessment and Conservation Engineering (RACE)
#     Midwater Assessment and Conservation Engineering (MACE)

#  THIS SOFTWARE AND ITS DOCUMENTATION ARE CONSIDERED TO BE IN THE PUBLIC DOMAIN
#  AND THUS ARE AVAILABLE FOR UNRESTRICTED PUBLIC USE. THEY ARE FURNISHED "AS IS."
#  THE AUTHORS, THE UNITED STATES GOVERNMENT, ITS INSTRUMENTALITIES, OFFICERS,
#  EMPLOYEES, AND AGENTS MAKE NO WARRANTY, EXPRESS OR IMPLIED, AS TO THE USEFULNESS
#  OF THE SOFTWARE AND DOCUMENTATION FOR ANY PURPOSE. THEY ASSUME NO RESPONSIBILITY
#  (1) FOR THE USE OF THE SOFTWARE AND DOCUMENTATION; OR (2) TO PROVIDE TECHNICAL
#  SUPPORT TO USERS.

import numpy as np

class NMEAData(object):
    """The nmea_data class provides storage for and parsing of NMEA data commonly
    collected along with sonar data.
    """

    def __init__(self):

        self.CHUNK_SIZE = 500

        # Create a counter to keep track of the number of datagrams, This is
        # used to inform the array sizes.
        self.n_raw = 0

        # Create arrays to store raw NMEA data as well as times, talkers,
        # and message IDs.
        self.raw_datagrams = np.empty(self.CHUNK_SIZE, dtype=object)
        self.nmea_times = np.empty(self.CHUNK_SIZE, dtype='datetime64[ms]')
        self.talkers = np.empty(self.CHUNK_SIZE, dtype='U2')
        self.messages = np.empty(self.CHUNK_SIZE, dtype='U3')

        # Create a couple of lists to store the unique talkers and message IDs.
        self.talker_ids = []
        self.message_ids = []

        # nmea_definitions define the NMEA message(s) and pynmea2.NMEASentence
        # attributes of those messages that the NMEA interpolation routine
        # will process. These definitions can also be used to define meta-types
        # which are data types that can be contained within multiple message
        # types and allow the user to request the meta-type without knowing
        # ahead of time if specific messages are contained within the data.
        self.nmea_definitions = {}

        self.nmea_definitions['GGA'] = {'message': ['GGA'],
                                        'fields': ['latitude', 'longitude']}
        self.nmea_definitions['GLL'] = {'message': ['GLL'],
                                        'fields': ['latitude', 'longitude']}
        self.nmea_definitions['RMC'] = {'message': ['RMC'],
                                        'fields': ['latitude', 'longitude']}
        self.nmea_definitions['HDT'] = {'message': ['HDT'],
                                        'fields': ['heading_true']}
        self.nmea_definitions['VTG'] = {'message': ['VTG'],
                                        'fields': ['true_track',
                                                   'spd_over_grnd_kts']}
        # Meta-types
        #
        #  Define a "position" meta-type that extracts lat/lon data from GGA, RMC
        #  and/or GLL datagrams if they are present.
        self.nmea_definitions['position'] = {'message': ['GGA','RMC', 'GLL'],
                                             'fields': ['latitude',
                                                        'longitude']}
        #  for these meta-types we don't have multiple messages and are simply
        #  using the meta-type to define an easy to remember label like "speed"
        #  or "distance" that you can use instead of 'VTG' or 'VLW'.
        self.nmea_definitions['speed']    = {'message': ['VTG'],
                                             'fields': ['spd_over_grnd_kts']}
        self.nmea_definitions['attitude'] = {'message': ['SHR'],
                                             'fields': ['heave', 'pitch',
                                                        'roll']}
        self.nmea_definitions['distance'] = {'message':['VLW'],
                                             'fields': ['trip_distance_nmi']}

    def add_datagram(self, time, text, allow_duplicates=False):
        """
        Add NMEA datagrams to this object.

        add_datagram adds a NMEA datagram to the class. It adds it to the
        raw_datagram list as well as parsing the header and adding the
        talker + mesage ID to the type_index dictionary.
        Args:
            time (datetime64): Timestamp of NMEA datagram (this is likely
                different than the timestamp within the NMEA string itself).
            text (str): The raw NMEA string.
            allow_duplicates (bool): When False, NMEA datagrams that share
                the same timestamp, talker ID, and message ID with an
                existing datagram will be discarded.

        """
        # Parse the NMEA message header
        header = str(text[1:6].upper())

        # Verify we have a plausible header and then process.
        if header.isalpha() and len(header) == 5:

            #  check if we're allowing duplicates and if this is one. We need
            #  to do this since .out files can contain duplicate NMEA data.
            if (not allow_duplicates):
                dup_idx = (self.nmea_times == time)
                if (np.any(dup_idx)):
                    #  We have a time match - check the talker and message id
                    my_talker = self.talkers[dup_idx]
                    my_message = self.messages[dup_idx]
                    if ((header[0:2] in my_talker) and (header[2:6] in my_message)):
                        #  this is the same - discard it
                        return

            # Increment datagram counter.
            self.n_raw += 1

            # Check if we need to resize our arrays. If so, resize arrays.
            if self.n_raw > self.nmea_times.shape[0]:
                self._resize_arrays(self.nmea_times.shape[0] +
                                    self.CHUNK_SIZE)

            # Add this datagram and associated data to our data arrays and
            # then Add the talker and message ID to our list of unique talkers
            # and messages.
            self.raw_datagrams[self.n_raw-1] = text
            self.nmea_times[self.n_raw-1] = time
            self.talkers[self.n_raw-1] = header[0:2]
            self.messages[self.n_raw-1] = header[2:6]

            if not header[0:2] in self.talker_ids:
                self.talker_ids.append(header[0:2])
            if not header[2:5] in self.message_ids:
                self.message_ids.append(header[2:5])

    def _resize_arrays(self, new_size):
        """
        Resize arrays if needed to hold more data.

        _resize_arrays expands our data arrays and is called when said arrays
        are filled with data and more data need to be added.

        Args:
            new_size (int): New size for arrays, Since these are all 1d
            arrays the value is simply an integer.

        """

        self.nmea_times = np.resize(self.nmea_times,(new_size))
        self.raw_datagrams = np.resize(self.raw_datagrams,(new_size))
        self.talkers = np.resize(self.talkers,(new_size))
        self.messages = np.resize(self.messages,(new_size))

    def trim(self):
        """
        Trim arrays to proper size after all data are added.

        trim is called when one is done adding data to the object. It
        removes empty elements of the data arrays.
        """

        self._resize_arrays(self.n_raw)

    def __str__(self):
        """
        Reimplemented string method that provides some basic info about the
        nmea_data object.

        """

        #  print the class and address
        msg = str(self.__class__) + " at " + str(hex(id(self))) + "\n"

        #  print some more info about the nmea_data instance
        if (self.n_raw > 0):
            msg = "{0}      NMEA data start time: {1}\n".format(
                                                       msg, self.nmea_times[0])
            msg = "{0}        NMEA data end time: {1}\n".format(
                                              msg,self.nmea_times[self.n_raw-1])
            msg = "{0}  number of NMEA datagrams: {1}\n".format(msg, self.n_raw)
            msg = "{0}         unique talker IDs: {1}\n".format(
                                               msg, (','.join(self.talker_ids)))
            msg = "{0}        unique message IDs: {1}\n".format(
                                              msg, (','.join(self.message_ids)))
            #  TODO: add reporting of numbers of individual message IDs
        else:
            msg = msg + ("  nmea_data object contains no data\n")

        return msg


text:
# Differentiable plasticity: maze exploration task.
#
# Copyright (c) 2018 Uber Technologies, Inc.
#
# 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.

# (Note: this contains an EXTREMELY EXPERIMENTAL implementation of an LSTM with plastic connections)

import argparse
import pdb
import torch
import torch.nn as nn
from torch.autograd import Variable
import numpy as np
from numpy import random
import torch.nn.functional as F
from torch import optim
from torch.optim import lr_scheduler
import random
import sys
import pickle
import time
import os
import platform

# Uber-only:
import OpusHdfsCopy
from OpusHdfsCopy import transferFileToHdfsDir, checkHdfs

import numpy as np
#import matplotlib.pyplot as plt
import glob

np.set_printoptions(precision=4)

ETA = .02  # Not used

ADDINPUT = 4 # 1 input for the previous reward, 1 input for numstep, 1 for whether currently on reward square, 1 "Bias" input

NBACTIONS = 4  # U, D, L, R

RFSIZE = 3 # Receptive field size

TOTALNBINPUTS =  RFSIZE * RFSIZE + ADDINPUT + NBACTIONS

##ttype = torch.FloatTensor;    # For CPU
ttype = torch.cuda.FloatTensor; # Gor GPU

class Network(nn.Module):
    def __init__(self, params):
        super(Network, self).__init__()
        self.rule = params['rule']
        self.type = params['type']
        self.softmax= torch.nn.functional.softmax
        if params['activ'] == 'tanh':
            self.activ = F.tanh
        elif params['activ'] == 'selu':
            self.activ = F.selu
        else:
            raise ValueError('Must choose an activ function')
        if params['type'] == 'lstm':
            self.lstm = torch.nn.LSTM(TOTALNBINPUTS, params['hiddensize']).cuda()
        elif params['type'] == 'rnn':
            self.i2h = torch.nn.Linear(TOTALNBINPUTS, params['hiddensize']).cuda()
            self.w =  torch.nn.Parameter((.01 * torch.rand(params['hiddensize'], params['hiddensize'])).cuda(), requires_grad=True)
        elif params['type'] == 'homo':
            self.i2h = torch.nn.Linear(TOTALNBINPUTS, params['hiddensize']).cuda()
            self.w =  torch.nn.Parameter((.01 * torch.rand(params['hiddensize'], params['hiddensize'])).cuda(), requires_grad=True)
            self.alpha = torch.nn.Parameter((.01 * torch.ones(1)).cuda(), requires_grad=True) # Homogenous plasticity: everyone has the same alpha
            self.eta = torch.nn.Parameter((.01 * torch.ones(1)).cuda(), requires_grad=True)   # Everyone has the same eta
        elif params['type'] == 'plastic':
            self.i2h = torch.nn.Linear(TOTALNBINPUTS, params['hiddensize']).cuda()
            self.w =  torch.nn.Parameter((.01 * torch.rand(params['hiddensize'], params['hiddensize'])).cuda(), requires_grad=True)
            self.alpha =  torch.nn.Parameter((.01 * torch.rand(params['hiddensize'], params['hiddensize'])).cuda(), requires_grad=True)
            self.eta = torch.nn.Parameter((.01 * torch.ones(1)).cuda(), requires_grad=True)  # Everyone has the same eta
        elif params['type'] == 'lstmplastic':   # LSTM with plastic connections. HIGHLY EXPERIMENTAL, NOT DEBUGGED
            self.h2f = torch.nn.Linear(params['hiddensize'], params['hiddensize']).cuda()
            self.h2i = torch.nn.Linear(params['hiddensize'], params['hiddensize']).cuda()
            self.h2opt = torch.nn.Linear(params['hiddensize'], params['hiddensize']).cuda()

            # Plasticity only in the recurrent connections, h to c.
            #self.h2c = torch.nn.Linear(params['hiddensize'], params['hiddensize']).cuda()  # This is replaced by the plastic connection matrices below
            self.w =  torch.nn.Parameter((.01 * torch.rand(params['hiddensize'], params['hiddensize'])).cuda(), requires_grad=True)
            self.alpha =  torch.nn.Parameter((.01 * torch.rand(params['hiddensize'], params['hiddensize'])).cuda(), requires_grad=True)
            self.eta = torch.nn.Parameter((.01 * torch.ones(1)).cuda(), requires_grad=True)  # Everyone has the same eta

            self.x2f = torch.nn.Linear(TOTALNBINPUTS, params['hiddensize']).cuda()
            self.x2opt = torch.nn.Linear(TOTALNBINPUTS, params['hiddensize']).cuda()
            self.x2i = torch.nn.Linear(TOTALNBINPUTS, params['hiddensize']).cuda()
            self.x2c = torch.nn.Linear(TOTALNBINPUTS, params['hiddensize']).cuda()
        elif params['type'] == 'lstmmanual':   # An LSTM implemented "by hand", to ensure maximum simlarity with the plastic LSTM
            self.h2f = torch.nn.Linear(params['hiddensize'], params['hiddensize']).cuda()
            self.h2i = torch.nn.Linear(params['hiddensize'], params['hiddensize']).cuda()
            self.h2opt = torch.nn.Linear(params['hiddensize'], params['hiddensize']).cuda()
            self.h2c = torch.nn.Linear(params['hiddensize'], params['hiddensize']).cuda()
            self.x2f = torch.nn.Linear(TOTALNBINPUTS, params['hiddensize']).cuda()
            self.x2opt = torch.nn.Linear(TOTALNBINPUTS, params['hiddensize']).cuda()
            self.x2i = torch.nn.Linear(TOTALNBINPUTS, params['hiddensize']).cuda()
            self.x2c = torch.nn.Linear(TOTALNBINPUTS, params['hiddensize']).cuda()
            ##fgt = F.sigmoid(self.x2f(input) + self.h2f(hidden[0]))
            ##ipt = F.sigmoid(self.x2i(input) + self.h2i(hidden[0]))
            ##opt = F.sigmoid(self.x2o(input) + self.h2o(hidden[0]))
            ##cell = torch.mul(fgt, hidden[1]) + torch.mul(ipt, F.tanh(self.x2c(input) + self.h2c(hidden[0])))
            ##h = torch.mul(opt, cell)
            ##hidden = (h, cell)
        else:
            raise ValueError("Which network type?")
        self.h2o = torch.nn.Linear(params['hiddensize'], NBACTIONS).cuda()  # From hidden to action output
        self.h2v = torch.nn.Linear(params['hiddensize'], 1).cuda()          # From hidden to value prediction (for A3C)
        self.params = params

        # Notice that the vectors are row vectors, and the matrices are transposed wrt the usual order, following apparent pytorch conventions
        # Each *column* of w targets a single output neuron

    def forward(self, input, hidden, hebb):
        if self.type == 'lstm':
            hactiv, hidden = self.lstm(input.view(1, 1, -1), hidden)  # hactiv is just the h. hidden is the h and the cell state, in a tuple
            hactiv = hactiv.view(1, -1)

        elif self.type == 'rnn':
            hactiv = self.activ(self.i2h(input) + hidden.mm(self.w))
            hidden = hactiv

        # Draft for a "manual" lstm:
        elif self.type== 'lstmmanual':
            # hidden[0] is the previous h state. hidden[1] is the previous c state
            fgt = F.sigmoid(self.x2f(input) + self.h2f(hidden[0]))
            ipt = F.sigmoid(self.x2i(input) + self.h2i(hidden[0]))
            opt = F.sigmoid(self.x2opt(input) + self.h2opt(hidden[0]))
            cell = torch.mul(fgt, hidden[1]) + torch.mul(ipt, F.tanh(self.x2c(input) + self.h2c(hidden[0])))
            hactiv = torch.mul(opt, F.tanh(cell))
            #pdb.set_trace()
            hidden = (hactiv, cell)
            if np.isnan(np.sum(hactiv.data.cpu().numpy())) or np.isnan(np.sum(hidden[1].data.cpu().numpy())) :
                raise ValueError("Nan detected !")

        elif self.type== 'lstmplastic':
            fgt = F.sigmoid(self.x2f(input) + self.h2f(hidden[0]))
            ipt = F.sigmoid(self.x2i(input) + self.h2i(hidden[0]))
            opt = F.sigmoid(self.x2opt(input) + self.h2opt(hidden[0]))
            #cell = torch.mul(fgt, hidden[1]) + torch.mul(ipt, F.tanh(self.x2c(input) + self.h2c(hidden[0])))

            #Need to think what the inputs and outputs should be for the
            #plasticity. It might be worth introducing an additional stage
            #consisting of whatever is multiplied by ift and then added to the
            #cell state, rather than the full cell state.... But we can
            #experiment both!

            #cell = torch.mul(fgt, hidden[1]) + torch.mul(ipt, F.tanh(self.x2c(input) + hidden[0].mm(self.w + torch.mul(self.alpha, hebb)))) #  self.h2c(hidden[0])))
            inputtocell =  F.tanh(self.x2c(input) + hidden[0].mm(self.w + torch.mul(self.alpha, hebb)))
            cell = torch.mul(fgt, hidden[1]) + torch.mul(ipt, inputtocell) #  self.h2c(hidden[0])))

            if self.rule == 'hebb':
                #hebb = (1 - self.eta) * hebb + self.eta * torch.bmm(hidden[0].unsqueeze(2), cell.unsqueeze(1))[0]
                hebb = (1 - self.eta) * hebb + self.eta * torch.bmm(hidden[0].unsqueeze(2), inputtocell.unsqueeze(1))[0]
            elif self.rule == 'oja':
                # NOTE: NOT SURE ABOUT THE OJA VERSION !!
                hebb = hebb + self.eta * torch.mul((hidden[0][0].unsqueeze(1) - torch.mul(hebb , inputtocell[0].unsqueeze(0))) , inputtocell[0].unsqueeze(0))  # Oja's rule. Remember that yin, yout are row vectors (dim (1,N)). Also, broadcasting!
                #hebb = hebb + self.eta * torch.mul((hidden[0].unsqueeze(1) - torch.mul(hebb , hactiv[0].unsqueeze(0))) , hactiv[0].unsqueeze(0))  # Oja's rule. Remember that yin, yout are row vectors (dim (1,N)). Also, broadcasting!
            hactiv = torch.mul(opt, F.tanh(cell))
            #pdb.set_trace()
            hidden = (hactiv, cell)
            if np.isnan(np.sum(hactiv.data.cpu().numpy())) or np.isnan(np.sum(hidden[1].data.cpu().numpy())) :
                raise ValueError("Nan detected !")

        elif self.type == 'plastic':
            hactiv = self.activ(self.i2h(input) + hidden.mm(self.w + torch.mul(self.alpha, hebb)))
            if self.rule == 'hebb':
                hebb = (1 - self.eta) * hebb + self.eta * torch.bmm(hidden.unsqueeze(2), hactiv.unsqueeze(1))[0]
            elif self.rule == 'oja':
                hebb = hebb + self.eta * torch.mul((hidden[0].unsqueeze(1) - torch.mul(hebb , hactiv[0].unsqueeze(0))) , hactiv[0].unsqueeze(0))  # Oja's rule. Remember that yin, yout are row vectors (dim (1,N)). Also, broadcasting!
            else:
                raise ValueError("Must specify learning rule ('hebb' or 'oja')")
            hidden = hactiv

        elif self.type == 'homo':
            hactiv = self.activ(self.i2h(input) + hidden.mm(self.w + self.alpha * hebb))
            if self.rule == 'hebb':
                hebb = (1 - self.eta) * hebb + self.eta * torch.bmm(hidden.unsqueeze(2), hactiv.unsqueeze(1))[0]
            elif self.rule == 'oja':
                hebb = hebb + self.eta * torch.mul((hidden[0].unsqueeze(1) - torch.mul(hebb , hactiv[0].unsqueeze(0))) , hactiv[0].unsqueeze(0))  # Oja's rule. Remember that yin, yout are row vectors (dim (1,N)). Also, broadcasting!
            else:
                raise ValueError("Must specify learning rule ('hebb' or 'oja')")
            hidden = hactiv

        activout = self.softmax(self.h2o(hactiv))   # Action selection
        valueout = self.h2v(hactiv)                 # Value prediction (for A3C)

        return activout, valueout, hidden, hebb

    def initialZeroHebb(self):
        return Variable(torch.zeros(self.params['hiddensize'], self.params['hiddensize']) , requires_grad=False).cuda()

    def initialZeroState(self):
        if self.params['type'] == 'lstm':
            return (Variable(torch.zeros(1, 1, self.params['hiddensize']), requires_grad=False).cuda() , Variable(torch.zeros(1, 1, self.params['hiddensize']), requires_grad=False ).cuda() )
        elif self.params['type'] == 'lstmmanual' or self.params['type'] == 'lstmplastic':
            return (Variable(torch.zeros(1, self.params['hiddensize']), requires_grad=False).cuda() , Variable(torch.zeros(1, self.params['hiddensize']), requires_grad=False ).cuda() )
        elif self.params['type'] == 'rnn' or self.params['type'] == 'plastic' or self.params['type'] == 'homo':
            return Variable(torch.zeros(1, self.params['hiddensize']), requires_grad=False ).cuda()
        else:
            raise ValueError("Which type?")

def train(paramdict):
    #params = dict(click.get_current_context().params)
    print("Starting training...")
    params = {}
    #params.update(defaultParams)
    params.update(paramdict)
    print("Passed params: ", params)
    print(platform.uname())
    #params['nbsteps'] = params['nbshots'] * ((params['prestime'] + params['interpresdelay']) * params['nbclasses']) + params['prestimetest']  # Total number of steps per episode
    suffix = "maze_"+"".join([str(x)+"_" if pair[0] is not 'nbsteps' and pair[0] is not 'rngseed' and pair[0] is not 'save_every' and pair[0] is not 'test_every' else '' for pair in sorted(zip(params.keys(), params.values()), key=lambda x:x[0] ) for x in pair])[:-1] + "_rngseed_" + str(params['rngseed'])   # Turning the parameters into a nice suffix for filenames

    # Initialize random seeds (first two redundant?)
    print("Setting random seeds")
    np.random.seed(params['rngseed']); random.seed(params['rngseed']); torch.manual_seed(params['rngseed'])

    print("Initializing network")
    net = Network(params)
    print ("Shape of all optimized parameters:", [x.size() for x in net.parameters()])
    allsizes = [torch.numel(x.data.cpu()) for x in net.parameters()]
    print ("Size (numel) of all optimized elements:", allsizes)
    print ("Total size (numel) of all optimized elements:", sum(allsizes))

    print("Initializing optimizer")
    optimizer = torch.optim.Adam(net.parameters(), lr=1.0*params['lr'], eps=1e-4)
    #scheduler = torch.optim.lr_scheduler.StepLR(optimizer, gamma=params['gamma'], step_size=params['steplr'])

    LABSIZE = params['labsize']
    lab = np.ones((LABSIZE, LABSIZE))
    CTR = LABSIZE // 2

    # Simple cross maze
    #lab[CTR, 1:LABSIZE-1] = 0
    #lab[1:LABSIZE-1, CTR] = 0

    # Double-T maze
    #lab[CTR, 1:LABSIZE-1] = 0
    #lab[1:LABSIZE-1, 1] = 0
    #lab[1:LABSIZE-1, LABSIZE - 2] = 0

    # Grid maze
    lab[1:LABSIZE-1, 1:LABSIZE-1].fill(0)
    for row in range(1, LABSIZE - 1):
        for col in range(1, LABSIZE - 1):
            if row % 2 == 0 and col % 2 == 0:
                lab[row, col] = 1
    lab[CTR,CTR] = 0 # Not really necessary, but nicer to not start on a wall, and perhaps helps localization by introducing a detectable irregularity in the center?

    all_losses = []
    all_losses_objective = []
    all_losses_eval = []
    all_losses_v = []
    lossbetweensaves = 0
    nowtime = time.time()

    print("Starting episodes...")
    sys.stdout.flush()

    pos = 0
    hidden = net.initialZeroState()
    hebb = net.initialZeroHebb()

    # Starting episodes!

    for numiter in range(params['nbiter']):

        PRINTTRACE = 0
        if (numiter+1) % (1 + params['print_every']) == 0:
            PRINTTRACE = 1

        # Note: it doesn't matter if the reward is on the center (reward is only computed after an action is taken). All we need is not to put it on a wall or pillar (lab=1)
        rposr = 0; rposc = 0
        if params['rp'] == 0:
            # If we want to constrain the reward to fall on the periphery of the maze
            while lab[rposr, rposc] == 1:
                rposr = np.random.randint(1, LABSIZE - 1)
                rposc = np.random.randint(1, LABSIZE - 1)
        elif params['rp'] == 1:
            while lab[rposr, rposc] == 1 or (rposr != 1 and rposr != LABSIZE -2 and rposc != 1 and rposc != LABSIZE-2):
                rposr = np.random.randint(1, LABSIZE - 1)
                rposc = np.random.randint(1, LABSIZE - 1)
        #print("Reward pos:", rposr, rposc)

        # Agent always starts an episode from the center
        posc = CTR
        posr = CTR

        optimizer.zero_grad()
        loss = 0
        lossv = 0
        hidden = net.initialZeroState()
        hebb = net.initialZeroHebb()

        reward = 0.0
        rewards = []
        vs = []
        logprobs = []
        sumreward = 0.0
        dist = 0

        for numstep in range(params['eplen']):

            inputsN = np.zeros((1, TOTALNBINPUTS), dtype='float32')
            inputsN[0, 0:RFSIZE * RFSIZE] = lab[posr - RFSIZE//2:posr + RFSIZE//2 +1, posc - RFSIZE //2:posc + RFSIZE//2 +1].flatten()

            inputs = torch.from_numpy(inputsN).cuda()
            # Previous chosen action
            #inputs[0][numactionchosen] = 1
            inputs[0][-1] = 1 # Bias neuron
            inputs[0][-2] = numstep
            inputs[0][-3] = reward

            # Running the network
            y, v, hidden, hebb = net(Variable(inputs, requires_grad=False), hidden, hebb)  # y  should output probabilities; v is the value prediction

            distrib = torch.distributions.Categorical(y)
            actionchosen = distrib.sample()  # sample() returns a Pytorch tensor of size 1; this is needed for the backprop below
            numactionchosen = actionchosen.data[0]    # Turn to scalar

            # Target position, based on the selected action
            tgtposc = posc
            tgtposr = posr
            if numactionchosen == 0:  # Up
                tgtposr -= 1
            elif numactionchosen == 1:  # Down
                tgtposr += 1
            elif numactionchosen == 2:  # Left
                tgtposc -= 1
            elif numactionchosen == 3:  # Right
                tgtposc += 1
            else:
                raise ValueError("Wrong Action")

            reward = 0.0
            if lab[tgtposr][tgtposc] == 1:
                reward = -.1
            else:
                dist += 1
                posc = tgtposc
                posr = tgtposr

            # Did we hit the reward location ? Increase reward and teleport!
            # Note that it doesn't matter if we teleport onto the reward, since reward hitting is only evaluated after the (obligatory) move
            if rposr == posr and rposc == posc:
                reward += 10
                if params['randstart'] == 1:
                    posr = np.random.randint(1, LABSIZE - 1)
                    posc = np.random.randint(1, LABSIZE - 1)
                    while lab[posr, posc] == 1:
                        posr = np.random.randint(1, LABSIZE - 1)
                        posc = np.random.randint(1, LABSIZE - 1)
                else:
                    posr = CTR
                    posc = CTR

            # Store the obtained reward, value prediction, and log-probabilities, for this time step
            rewards.append(reward)
            sumreward += reward
            vs.append(v)
            logprobs.append(distrib.log_prob(actionchosen))

            # A3C/A2C has an entropy reward on the output probabilities, to
            # encourage exploration. Our version of PyTorch does not have an
            # entropy() function for Distribution, so we use a penalty on the
            # sum of squares instead, which has the same basic property
            # (discourages concentration). It really does help!
            loss += params['bentropy'] * y.pow(2).sum()

            #if PRINTTRACE:
            #    print("Probabilities:", y.data.cpu().numpy(), "Picked action:", numactionchosen, ", got reward", reward)

        # Do the A2C ! (essentially copied from V. Mnih, https://arxiv.org/abs/1602.01783, Algorithm S3)
        R = 0
        gammaR = params['gr']
        for numstepb in reversed(range(params['eplen'])) :
            R = gammaR * R + rewards[numstepb]
            lossv += (vs[numstepb][0] - R).pow(2)
            loss -= logprobs[numstepb] * (R - vs[numstepb].data[0][0])

        if PRINTTRACE:
            print("lossv: ", lossv.data.cpu().numpy()[0])
            print ("Total reward for this episode:", sumreward, "Dist:", dist)

        # Do we want to squash rewards for stabilization?
        if params['squash'] == 1:
            if sumreward < 0:
                sumreward = -np.sqrt(-sumreward)
            else:
                sumreward = np.sqrt(sumreward)
        elif params['squash'] == 0:
            pass
        else:
            raise ValueError("Incorrect value for squash parameter")

        # Mixing the reward loss and the value-prediction loss
        loss += params['blossv'] * lossv
        loss /= params['eplen']
        loss.backward()

        #scheduler.step()
        optimizer.step()
        #torch.cuda.empty_cache()

        lossnum = loss.data[0]
        lossbetweensaves += lossnum
        if (numiter + 1) % 10 == 0:
            all_losses_objective.append(lossnum)
            all_losses_eval.append(sumreward)
            all_losses_v.append(lossv.data[0])

        # Algorithm done. Now print statistics and save files.

        if (numiter+1) % params['print_every'] == 0:

            print(numiter, "====")
            print("Mean loss: ", lossbetweensaves / params['print_every'])
            lossbetweensaves = 0
            previoustime = nowtime
            nowtime = time.time()
            print("Time spent on last", params['print_every'], "iters: ", nowtime - previoustime)
            if params['type'] == 'plastic' or params['type'] == 'lstmplastic':
                print("ETA: ", net.eta.data.cpu().numpy(), "alpha[0,1]: ", net.alpha.data.cpu().numpy()[0,1], "w[0,1]: ", net.w.data.cpu().numpy()[0,1] )
            elif params['type'] == 'rnn':
                print("w[0,1]: ", net.w.data.cpu().numpy()[0,1] )

        if (numiter+1) % params['save_every'] == 0:
            print("Saving files...")
            losslast100 = np.mean(all_losses_objective[-100:])
            print("Average loss over the last 100 episodes:", losslast100)
            print("Saving local files...")
            with open('params_'+suffix+'.dat', 'wb') as fo:
                    pickle.dump(params, fo)
            with open('lossv_'+suffix+'.txt', 'w') as thefile:
                for item in all_losses_v:
                        thefile.write("%s\n" % item)
            with open('loss_'+suffix+'.txt', 'w') as thefile:
                for item in all_losses_eval:
                        thefile.write("%s\n" % item)
            torch.save(net.state_dict(), 'torchmodel_'+suffix+'.dat')
            # Uber-only
            print("Saving HDFS files...")
            if checkHdfs():
                print("Transfering to HDFS...")
                transferFileToHdfsDir('loss_'+suffix+'.txt', '/ailabs/tmiconi/gridlab/')
                transferFileToHdfsDir('torchmodel_'+suffix+'.dat', '/ailabs/tmiconi/gridlab/')
                transferFileToHdfsDir('params_'+suffix+'.dat', '/ailabs/tmiconi/gridlab/')
            #print("Saved!")

if __name__ == "__main__":
    parser = argparse.ArgumentParser()
    parser.add_argument("--rngseed", type=int, help="random seed", default=0)
    #parser.add_argument("--clamp", type=float, help="maximum (absolute value) gradient for clamping", default=1000000.0)
    parser.add_argument("--bentropy", type=float, help="coefficient for the A2C 'entropy' reward (really Simpson index concentration measure)", default=0.1)
    parser.add_argument("--blossv", type=float, help="coefficient for the A2C value prediction loss", default=.03)
    parser.add_argument("--labsize", type=int, help="size of the labyrinth; must be odd", default=9)
    parser.add_argument("--randstart", type=int, help="when hitting reward, should we teleport to random location (1) or center (0)?", default=1)
    parser.add_argument("--rp", type=int, help="whether the reward should be on the periphery", default=0)
    parser.add_argument("--squash", type=int, help="squash reward through signed sqrt (1 or 0)", default=0)
    #parser.add_argument("--nbarms", type=int, help="number of arms", default=2)
    #parser.add_argument("--nbseq", type=int, help="number of sequences between reinitializations of hidden/Hebbian state and position", default=3)
    parser.add_argument("--activ", help="activ function ('tanh' or 'selu')", default='tanh')
    parser.add_argument("--rule", help="learning rule ('hebb' or 'oja')", default='oja')
    parser.add_argument("--type", help="network type ('rnn' or 'plastic')", default='rnn')
    parser.add_argument("--gr", type=float, help="gammaR: discounting factor for rewards", default=.9)
    parser.add_argument("--lr", type=float, help="learning rate (Adam optimizer)", default=1e-4)
    parser.add_argument("--eplen", type=int, help="length of episodes", default=250)
    parser.add_argument("--hiddensize", type=int, help="size of the recurrent (hidden) layer", default=200)
    #parser.add_argument("--steplr", type=int, help="duration of each step in the learning rate annealing schedule", default=100000000)
    #parser.add_argument("--gamma", type=float, help="learning rate annealing factor", default=0.3)
    parser.add_argument("--nbiter", type=int, help="number of learning cycles", default=1000000)
    parser.add_argument("--save_every", type=int, help="number of cycles between successive save points", default=200)
    parser.add_argument("--print_every", type=int, help="number of cycles between successive printing of information", default=100)
    #parser.add_argument("--", type=int, help="", default=1e-4)
    args = parser.parse_args(); argvars = vars(args); argdict =  { k : argvars[k] for k in argvars if argvars[k] != None }
    #train()
    train(argdict)

max_stars_repo_name: ProjectAGI/differentiable-plasticity

max_issues_repo_name: ProjectAGI/differentiable-plasticity

max_forks_repo_name: ProjectAGI/differentiable-plasticity

reason: import numpy,from numpy


text:
import torch
from torch.nn import functional as F

import matplotlib.pylab as plt
import seaborn as sea
sea.set_style("whitegrid")

from DeepProteinDocking2D.Models.BruteForce.utility_functions import read_pkl
from DeepProteinDocking2D.Models.BruteForce.validation_metrics import RMSD
from DeepProteinDocking2D.Models.BruteForce.utility_functions import plot_assembly
import numpy as np

class TorchDockingFilter:
    def __init__(self):
        self.dim = 100
        self.num_angles = 360
        self.angles = torch.from_numpy(np.linspace(-np.pi, np.pi, num=self.num_angles)).cuda()

    def encode_transform(self, gt_rot, gt_txy):
        empty_3D = torch.zeros([self.num_angles, self.dim, self.dim], dtype=torch.double).cuda()
        deg_index_rot = (((gt_rot * 180.0/np.pi) + 180.0) % self.num_angles).type(torch.long)
        centered_txy = gt_txy.type(torch.long)

        empty_3D[deg_index_rot, centered_txy[0], centered_txy[1]] = 1
        target_flatindex = torch.argmax(empty_3D.flatten()).cuda()
        # print(gt_rot, gt_txy)
        # print(deg_index_rot, centered_txy)
        # print(ConcaveTrainer().extract_transform(empty_3D))
        return target_flatindex

    def extract_transform(self, pred_score):
        pred_argmax = torch.argmax(pred_score)
        pred_rot = ((pred_argmax / self.dim ** 2) * np.pi / 180.0) - np.pi

        XYind = torch.remainder(pred_argmax, self.dim ** 2)
        pred_X = XYind // self.dim
        pred_Y = XYind % self.dim
        if pred_X > self.dim//2:
            pred_X = pred_X - self.dim
        if pred_Y > self.dim//2:
            pred_Y = pred_Y - self.dim
        return pred_rot, torch.stack((pred_X, pred_Y), dim=0)

    @staticmethod
    def make_boundary(grid_shape):
        epsilon = 1e-5
        sobel_top = torch.tensor([[1, 2, 1], [0, 0, 0], [-1, -2, -1]], dtype=torch.float, requires_grad=True).cuda()
        sobel_left = sobel_top.t()

        feat_top = F.conv2d(grid_shape.unsqueeze(0).unsqueeze(0), weight=sobel_top.unsqueeze(0).unsqueeze(0), padding=1)
        feat_left = F.conv2d(grid_shape.unsqueeze(0).unsqueeze(0), weight=sobel_left.unsqueeze(0).unsqueeze(0), padding=1)

        top = feat_top.squeeze() + epsilon
        right = feat_left.squeeze() + epsilon
        boundary = torch.sqrt(top ** 2 + right ** 2)
        feat_stack = torch.stack([grid_shape, boundary], dim=0)

        return feat_stack

    @staticmethod
    def rotate(repr, angle):
        alpha = angle.detach()
        T0 = torch.stack([torch.cos(alpha), -torch.sin(alpha), torch.zeros_like(alpha)], dim=1)
        T1 = torch.stack([torch.sin(alpha), torch.cos(alpha), torch.zeros_like(alpha)], dim=1)
        R = torch.stack([T0, T1], dim=1)
        curr_grid = F.affine_grid(R, size=repr.size(), align_corners=True).type(torch.float)
        return F.grid_sample(repr, curr_grid, align_corners=True)

    def dock_global(self, receptor, ligand, weight_bulk=1.0, weight_bound=1.0, weight_crossterm1=1.0, weight_crossterm2=1.0, debug=False):
        initbox_size = receptor.shape[-1]
        # print(receptor.shape)
        pad_size = initbox_size // 2

        if initbox_size % 2 == 0:
            receptor = F.pad(receptor, pad=([pad_size, pad_size, pad_size, pad_size]), mode='constant', value=0)
            ligand = F.pad(ligand, pad=([pad_size, pad_size, pad_size, pad_size]), mode='constant', value=0)
        else:
            receptor = F.pad(receptor, pad=([pad_size, pad_size+1, pad_size, pad_size+1]), mode='constant', value=0)
            ligand = F.pad(ligand, pad=([pad_size, pad_size+1, pad_size, pad_size+1]), mode='constant', value=0)
        # print('padded shape', receptor.shape)

        f_rec = receptor.unsqueeze(0).repeat(self.num_angles,1,1,1)
        f_lig = ligand.unsqueeze(0).repeat(self.num_angles,1,1,1)
        rot_lig = TorchDockingFilter().rotate(f_lig, self.angles)

        if debug:
            with torch.no_grad():
                step = 30
                for i in range(rot_lig.shape[0]):
                    if i % step == 0:
                        plt.title('Torch '+str(i)+' degree rotation')
                        plt.imshow(rot_lig[i,0,:,:].detach().cpu())
                        plt.show()

        score = TorchDockingFilter().CE_dock_translations(f_rec, rot_lig, weight_bulk, weight_bound, weight_crossterm1, weight_crossterm2)

        return score

    def CE_dock_translations(self, receptor, ligand, weight_bulk=2.8, weight_bound=3.0, weight_crossterm1=-0.3, weight_crossterm2=-0.3):
        box_size = receptor.shape[-1]

        receptor_bulk, receptor_bound = torch.chunk(receptor, chunks=2, dim=1)
        ligand_bulk, ligand_bound = torch.chunk(ligand, chunks=2, dim=1)
        receptor_bulk = receptor_bulk.squeeze()
        receptor_bound = receptor_bound.squeeze()
        ligand_bulk = ligand_bulk.squeeze()
        ligand_bound = ligand_bound.squeeze()
        # print(receptor_bulk.shape)

        signal_dim = 2
        # Bulk score
        cplx_rec = torch.rfft(receptor_bulk, signal_ndim=signal_dim)
        cplx_lig = torch.rfft(ligand_bulk, signal_ndim=signal_dim)
        re = cplx_rec[:, :, :, 0] * cplx_lig[:, :, :, 0] + cplx_rec[:, :, :, 1] * cplx_lig[:, :, :, 1]
        im = -cplx_rec[:, :, :, 0] * cplx_lig[:, :, :, 1] + cplx_rec[:, :, :, 1] * cplx_lig[:, :, :, 0]
        cconv = torch.stack([re, im], dim=3)
        trans_bulk = torch.irfft(cconv, signal_dim, signal_sizes=(box_size, box_size))

        # print(cplx_rec.shape, cplx_lig.shape)

        # Boundary score
        cplx_rec = torch.rfft(receptor_bound, signal_ndim=signal_dim)
        cplx_lig = torch.rfft(ligand_bound, signal_ndim=signal_dim)
        re = cplx_rec[:, :, :, 0] * cplx_lig[:, :, :, 0] + cplx_rec[:, :, :, 1] * cplx_lig[:, :, :, 1]
        im = -cplx_rec[:, :, :, 0] * cplx_lig[:, :, :, 1] + cplx_rec[:, :, :, 1] * cplx_lig[:, :, :, 0]
        cconv = torch.stack([re, im], dim=3)
        trans_bound = torch.irfft(cconv, signal_dim, signal_sizes=(box_size, box_size))
        # print(trans_bound.shape)

        # Boundary - bulk score
        cplx_rec = torch.rfft(receptor_bound, signal_ndim=signal_dim)
        cplx_lig = torch.rfft(ligand_bulk, signal_ndim=signal_dim)
        re = cplx_rec[:, :, :, 0] * cplx_lig[:, :, :, 0] + cplx_rec[:, :, :, 1] * cplx_lig[:, :, :, 1]
        im = -cplx_rec[:, :, :, 0] * cplx_lig[:, :, :, 1] + cplx_rec[:, :, :, 1] * cplx_lig[:, :, :, 0]
        cconv = torch.stack([re, im], dim=3)
        trans_bulk_bound = torch.irfft(cconv, signal_dim, signal_sizes=(box_size, box_size))

        # Bulk - boundary score
        cplx_rec = torch.rfft(receptor_bulk, signal_ndim=signal_dim)
        cplx_lig = torch.rfft(ligand_bound, signal_ndim=signal_dim)
        re = cplx_rec[:, :, :, 0] * cplx_lig[:, :, :, 0] + cplx_rec[:, :, :, 1] * cplx_lig[:, :, :, 1]
        im = -cplx_rec[:, :, :, 0] * cplx_lig[:, :, :, 1] + cplx_rec[:, :, :, 1] * cplx_lig[:, :, :, 0]
        cconv = torch.stack([re, im], dim=3)
        trans_bound_bulk = torch.irfft(cconv, signal_dim, signal_sizes=(box_size, box_size))

        ## cross-term score maximizing
        score = weight_bound * trans_bound + weight_crossterm1 * trans_bulk_bound + weight_crossterm2 * trans_bound_bulk - weight_bulk * trans_bulk

        # print(score.shape)

        return score

    @staticmethod
    def check_FFT_predictions(pred_score, receptor, ligand, gt_rot, gt_txy):
        print('*'*50)

        gt_rot = torch.tensor(gt_rot, dtype=torch.float).cuda()
        gt_txy = torch.tensor(gt_txy, dtype=torch.float).cuda()
        # print(gt_rot)

        pred_rot, pred_txy = TorchDockingFilter().extract_transform(pred_score)
        rmsd_out = RMSD(ligand, gt_rot, gt_txy, pred_rot, pred_txy).calc_rmsd()
        print('extracted predicted indices', pred_rot, pred_txy)
        print('gt indices', gt_rot, gt_txy)
        print('RMSD', rmsd_out.item())
        print()

        pair = plot_assembly(receptor.detach().cpu(), ligand.detach().cpu().numpy(), pred_rot.detach().cpu().numpy(),
                             pred_txy.detach().cpu().numpy(), gt_rot.detach().cpu().numpy(), gt_txy.detach().cpu().numpy())
        plt.imshow(pair.transpose())
        plt.show()

if __name__ == '__main__':

    data = read_pkl('toy_concave_data/docking_data_train')

    for i in range(len(data)):
        receptor, ligand, gt_txy, gt_rot = data[i]

        init_dim = receptor.shape[-1]
        ### print(receptor.shape)
        if init_dim > 51:
            receptor = receptor[25:75, 25:75]
            ligand = ligand[25:75, 25:75]

        receptor = torch.tensor(receptor, dtype=torch.float).cuda()
        ligand = torch.tensor(ligand, dtype=torch.float).cuda()
        receptor_stack = TorchDockingFilter().make_boundary(receptor)
        ligand_stack = TorchDockingFilter().make_boundary(ligand)
        pred_score = TorchDockingFilter().dock_global(receptor_stack, ligand_stack, debug=False)

        TorchDockingFilter().check_FFT_predictions(pred_score, receptor, ligand, gt_rot, gt_txy)


text:
import numpy as np
from scipy import linalg as sp_la

def tile_add_overlap(array, n, m=1):
    if n <= 1:
        return array
    else:
        add_on = np.copy(array[:-m])
        add_on[:m] += array[-m:]
        last = np.copy(array)
        last[:m] += array[-m:]
        return np.hstack((array[:-m], np.tile(add_on, n-2), last))

def overlapping_vector(vector, reps, overlap_size):
    n_vector = len(vector)
    n_result = reps*n_vector - (reps-1)*overlap_size
    # n_result = reps * (n_vector - overlap_size) + overlap_size
    result = np.zeros(n_result)
    non_overlap_size = int(n_vector-overlap_size)
    for i in range(reps):
        start_id = i * non_overlap_size
        end_id = start_id + n_vector
        result[start_id:end_id] += vector
    return result

def block_diag_overlap(block_list, overlap):
    m_sblks = [block.shape[0] for block in block_list]
    n_sblks = [block.shape[1] for block in block_list]
    n_blocks = len(block_list)
    m_ovl = overlap[0]
    n_ovl = overlap[1]
    m_final = np.sum(np.asarray(m_sblks)) - (n_blocks - 1) * m_ovl
    n_final = np.sum(np.asarray(n_sblks)) - (n_blocks - 1) * n_ovl
    block_array = np.zeros((m_final, n_final))
    for i, block in enumerate(block_list):
        x_id_start = i * m_sblks[i - 1] - i * m_ovl
        x_id_end = (i + 1) * m_sblks[i] - i * m_ovl
        y_id_start = i * n_sblks[i - 1] - i * n_ovl
        y_id_end = (i + 1) * n_sblks[i] - i * n_ovl
        block_array[x_id_start:x_id_end, y_id_start:y_id_end] += block[:, :]
    return block_array

def build_1d_conductance_matrix(cond_vector, offset=1):
    n_layer = len(cond_vector) + 1
    center_diag = overlapping_vector(cond_vector, 2, n_layer-2)
    center_diag *= -1.0
    off_diag = np.copy(cond_vector)
    return np.diag(center_diag, k=0) \
        + np.diag(off_diag, k=-offset) \
        + np.diag(off_diag, k=offset)

def build_z_cell_conductance_matrix(cond_vector, n_ele):
    list_mat = []
    for i in range(n_ele):
        list_mat.append(build_1d_conductance_matrix(cond_vector[:, i]))
    return sp_la.block_diag(*list_mat)

def build_x_cell_conductance_matrix(cond_vector, n_ele, n_layer=None):
    if not n_ele > 1:
        raise ValueError('x-conductance matrix can only be built for n_ele > 1')
    if n_layer is None:
        n_layer = len(cond_vector)
    center_diag = np.concatenate([cond_vector[:, i] for i in range(n_ele)])
    center_diag[n_layer:-n_layer] *= 2.0
    center_diag *= -1.0
    off_diag = np.concatenate([cond_vector[:, i] for i in range(n_ele-1)])
    return np.diag(center_diag, k=0) \
        + np.diag(off_diag, k=-n_layer) \
        + np.diag(off_diag, k=n_layer)

def build_cell_conductance_matrix(x_cond_vector, z_cond_vector, n_ele):
    z_cond_mtx = build_z_cell_conductance_matrix(z_cond_vector, n_ele)
    if n_ele > 1:
        x_cond_mtx = build_x_cell_conductance_matrix(x_cond_vector, n_ele)
    else:
        x_cond_mtx = 0.0
    return x_cond_mtx + z_cond_mtx

def connect_cells(matrix, cell_ids, layer_ids, values, mtx_ids,
                  replace=False):
    if np.isscalar(values):
        values = np.full(len(cell_ids), values)
    if not len(cell_ids) == len(layer_ids):
        raise ValueError('cell and layer index lists must have equal length')
    for i in range(len(cell_ids)):
        mtx_id_0 = mtx_ids[cell_ids[i, 0]][:][layer_ids[i, 0]]
        mtx_id_1 = mtx_ids[cell_ids[i, 1]][:][layer_ids[i, 1]]
        if replace:
            matrix[mtx_id_0, mtx_id_1] = values[i]
            matrix[mtx_id_0, mtx_id_0] = -values[i]
            matrix[mtx_id_1, mtx_id_1] = -values[i]
            matrix[mtx_id_1, mtx_id_0] = values[i]
        else:
            matrix[mtx_id_0, mtx_id_1] += values[i]
            matrix[mtx_id_0, mtx_id_0] += -values[i]
            matrix[mtx_id_1, mtx_id_1] += -values[i]
            matrix[mtx_id_1, mtx_id_0] += values[i]

def create_index_lists(cells):
    n_cells = len(cells)
    index_list = []
    layer_ids = [[] for _ in range(cells[-1].n_layer)]
    for i, cell in enumerate(cells):
        index_array = \
            (cells[i-1].n_ele * cells[i-1].n_layer) * i \
            + cell.index_array
        index_list.append(index_array.tolist())

    for i in range(n_cells):
        for j in range(cells[i].n_layer):
            layer_ids[j].append(index_list[i][j])
    layer_index_list = []
    for sub_list in layer_ids:
        layer_index_list.append(np.hstack(sub_list))
    return index_list, layer_index_list

max_stars_repo_name: lufire/pemfc-core

max_issues_repo_name: lufire/pemfc-core

max_forks_repo_name: lufire/pemfc-core


text:
import pytest
import numpy as np
import torch

from experiments.online_fa import (
    run_all_fa_experiments,
    generate_and_sample_fa_model,
    generate_fa_model,
    sample_fa_observations,
    learn_fa_with_sklearn,
    learn_fa_with_online_gradients,
    learn_fa_with_online_em,
)

@pytest.mark.parametrize("n_experiments", [1, 3])
@pytest.mark.parametrize("n_trials", [1, 4])
@pytest.mark.parametrize("n_sample_sizes", [1, 4])
def test_fa_results_rows_and_columns(n_experiments, n_trials, n_sample_sizes):
    config = dict(
        observation_dim=10,
        latent_dim=5,
        spectrum_range=[0, 1],
        n_samples=list((np.arange(n_sample_sizes) + 1) * 10),
    )
    experiments_config = [config] * n_experiments
    results = run_all_fa_experiments(
        experiments_config,
        n_trials,
        gradient_optimiser='sgd',
        gradient_optimiser_kwargs=dict(lr=0.01),
        gradient_warm_up_time_steps=1,
        em_warm_up_time_steps=1,
    )
    expected_columns = [
        'observation_dim',
        'latent_dim',
        'spectrum_min',
        'spectrum_max',
        'n_samples',
        'covar_norm',
        'covar_distance_sklearn',
        'covar_distance_online_gradient',
        'covar_distance_online_em',
        'wasserstein_sklearn',
        'wasserstein_online_gradient',
        'wasserstein_online_em',
        'experiment',
        'trial'
    ]
    actual_columns = list(results.columns)
    assert len(actual_columns) == len(expected_columns)
    assert len(np.intersect1d(actual_columns, expected_columns)) == len(actual_columns)

    assert len(results) == n_experiments * n_trials * n_sample_sizes
    for i in range(n_experiments):
        assert (results['experiment'] == i + 1).sum() == n_trials * n_sample_sizes

@pytest.mark.parametrize("observation_dim", [10, 20])
@pytest.mark.parametrize("latent_dim", [5, 8])
@pytest.mark.parametrize("spectrum_range", [[0, 1], [0, 10]])
@pytest.mark.parametrize('n_samples', [10, 100])
def test_true_params_shape(observation_dim, latent_dim, spectrum_range, n_samples):
    mean, F, psi, covar, observations = generate_and_sample_fa_model(
        observation_dim, latent_dim, spectrum_range, n_samples, random_seed=0,
    )
    assert mean.shape == (observation_dim,)
    assert covar.shape == (observation_dim, observation_dim)

@pytest.mark.parametrize("observation_dim", [10, 20])
@pytest.mark.parametrize("latent_dim", [5, 8])
@pytest.mark.parametrize("spectrum_range", [[0, 1], [0, 10]])
def test_model_mean_matches_sample_mean(observation_dim, latent_dim, spectrum_range):
    n_samples = 100000
    c, F, psi, covar, observations = generate_and_sample_fa_model(
        observation_dim, latent_dim, spectrum_range, n_samples, random_seed=0,
    )
    sample_mean = observations.mean(dim=0)
    normalised_distance = torch.linalg.norm(c - sample_mean) / torch.linalg.norm(c)
    assert normalised_distance < 0.01

@pytest.mark.parametrize("observation_dim", [10, 20])
@pytest.mark.parametrize("latent_dim", [5, 8])
@pytest.mark.parametrize("spectrum_range", [[0, 1], [0, 10]])
def test_model_covariance_matches_sample_covariance(observation_dim, latent_dim, spectrum_range):
    n_samples = 100000
    c, F, psi, covar, observations = generate_and_sample_fa_model(
        observation_dim, latent_dim, spectrum_range, n_samples, random_seed=0,
    )
    sample_covar = torch.from_numpy(np.cov(observations.t().numpy())).float()
    normalised_distance = torch.linalg.norm(covar - sample_covar) / torch.linalg.norm(covar)
    assert normalised_distance < 0.1

@pytest.mark.parametrize("observation_dim", [10, 20])
@pytest.mark.parametrize("latent_dim", [5, 8])
@pytest.mark.parametrize("spectrum_range", [[0, 1], [0, 10]])
@pytest.mark.parametrize('n_samples', [10, 100])
def test_sampled_fa_observations_shape(observation_dim, latent_dim, spectrum_range, n_samples):
    c, F, psi = generate_fa_model(observation_dim, latent_dim, spectrum_range, random_seed=0)
    observations = sample_fa_observations(c, F, psi, n_samples, random_seed=0)
    assert observations.shape == (n_samples, observation_dim)

@pytest.mark.parametrize("observation_dim", [10, 20])
@pytest.mark.parametrize("latent_dim", [5, 8])
@pytest.mark.parametrize("spectrum_range", [[0, 1], [0, 10]])
@pytest.mark.parametrize('n_samples', [10, 100])
def test_sklearn_learned_params_shape(observation_dim, latent_dim, spectrum_range, n_samples):
    c, F, psi, covar, observations = generate_and_sample_fa_model(
        observation_dim, latent_dim, spectrum_range, n_samples, random_seed=0,
    )
    mean, covar = learn_fa_with_sklearn(observations, latent_dim, random_seed=0)
    assert mean.shape == (observation_dim,)
    assert covar.shape == (observation_dim, observation_dim)

@pytest.mark.parametrize("observation_dim", [10, 20])
@pytest.mark.parametrize("latent_dim", [5, 8])
@pytest.mark.parametrize("spectrum_range", [[0, 1], [0, 10]])
@pytest.mark.parametrize('n_samples', [10, 100])
def test_online_gradients_learned_params_shape(observation_dim, latent_dim, spectrum_range, n_samples):
    c, F, psi, covar, observations = generate_and_sample_fa_model(
        observation_dim, latent_dim, spectrum_range, n_samples, random_seed=0,
    )
    _, mean, covar = learn_fa_with_online_gradients(
        observations, latent_dim, optimiser_name='sgd', optimiser_kwargs=None, n_warm_up_time_steps=1, random_seed=0,
    )
    assert mean.shape == (observation_dim,)
    assert covar.shape == (observation_dim, observation_dim)

@pytest.mark.parametrize("observation_dim", [10, 20])
@pytest.mark.parametrize("latent_dim", [5, 8])
@pytest.mark.parametrize("spectrum_range", [[0, 1], [0, 10]])
@pytest.mark.parametrize('n_samples', [10, 100])
def test_online_em_learned_params_shape(observation_dim, latent_dim, spectrum_range, n_samples):
    c, F, psi, covar, observations = generate_and_sample_fa_model(
        observation_dim, latent_dim, spectrum_range, n_samples, random_seed=0,
    )
    _, mean, covar = learn_fa_with_online_em(
        observations, latent_dim, n_warm_up_time_steps=1, random_seed=0,
    )
    assert mean.shape == (observation_dim,)
    assert covar.shape == (observation_dim, observation_dim)

max_stars_repo_name: ayeright/swafa

max_issues_repo_name: ayeright/swafa

max_forks_repo_name: ayeright/swafa


text:
import gym

env = gym.make("MountainCar-v0")    # this creates an environment with a mountain, and a car
env.reset()

done = False

while not done:
    action = 2      # This environment has three actions, 0 = push car left, 1 = do nothing, 2 = push car right
    new_state, reward, done, _ = env.step(action)

    env.render()    # rendering the GUI

env.close()

DISCRETE_OS_SIZE = [20] * len(env.observation_space.high)

discrete_os_win_size = (env.observation_space.high-env.observation_space.low) / DISCRETE_OS_SIZE

# Initialising the Q-table :
import numpy as np
q_table = np.random.uniform(low = -2, high = 0, size = (DISCRETE_OS_SIZE + [env.action_space.n]))

print(q_table.shape)


text:
# vim: expandtab tabstop=4 shiftwidth=4

from pathlib import Path
from tempfile import TemporaryDirectory

import unittest

from numpy import array as nparray

import dspftw

class SaveLoadDataTests(unittest.TestCase):
    def try_save_load(self, data, data_type, endianness):
        with TemporaryDirectory() as tempdir:
            temp_path = Path(tempdir)
            file_path = temp_path.joinpath('data_file')

            num_written = dspftw.save_data(str(file_path), data, data_type, endianness)
            self.assertEqual(num_written, len(data))

            loaded_data = dspftw.load_data(str(file_path), data_type, endianness)

            for i in range(len(loaded_data)):
                self.assertEqual(loaded_data[i], data[i])

    def test_save_load(self):
        data = nparray(range(127))  # wrap-around causes values above 127 to fail for int8
        data_types = ('uint8', 'uint16', 'uint32', 'uint64',
                      'int8', 'int16', 'int32', 'int64',
                      'float32', 'float64', 'float128')
        endiannesses = ('little', 'big')

        for data_type in data_types:
            for endianness in endiannesses:
                self.try_save_load(data, data_type, endianness)

max_stars_repo_name: dspftw/dspftw

max_issues_repo_name: dspftw/dspftw

max_forks_repo_name: dspftw/dspftw

reason: from numpy


text:
# -*- coding: utf-8 -*-
# ProDy: A Python Package for Protein Dynamics Analysis
#
# Copyright (C) 2010-2012 Ahmet Bakan
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>

"""This module defines classes and functions for protein dynamics analysis.

Dynamics Models
===============

Following classes are designed for modeling and analysis of protein dynamics:

  * :class:`~.ANM` - Anisotropic network model, for coarse-grained NMA
  * :class:`~.GNM` - Gaussian network model, for coarse-grained dynamics
    analysis
  * :class:`~.PCA` - Principal component analysis of conformation ensembles
  * :class:`~.EDA` - Essential dynamics analysis of dynamics trajectories
  * :class:`~.NMA` - Normal mode analysis, for analyzing data from external
    programs

Usage of these classes are shown in :ref:`anm`, :ref:`gnm`, :ref:`pca`, and
:ref:`eda` examples.

Following classes are for analysis of individual modes or subsets of modes:

  * :class:`~.Mode` - analyze individual normal/principal/essential modes
  * :class:`~.ModeSet` - analyze subset of modes from a dynamics model
  * :class:`~.Vector` - analyze modified modes or deformation vectors

Customize ENMs
==============

Following classes allow for using structure or distance based, or other custom
force constants and cutoff distances in :class:`~.ANM` and :class:`~.GNM`
calculations:

  * :class:`~.Gamma` - base class for developing property custom force
    constant calculation methods
  * :class:`~.GammaStructureBased` - secondary structure based force constants
  * :class:`~.GammaVariableCutoff` - atom type based variable cutoff function

Function library
================

Dynamics of the functions described below accept a *modes* argument (may also
appear in different names), which may refer to one or more of the following:

  * a dynamics model, :class:`~.ANM`, :class:`~.GNM`, :class:`~.NMA`,
    :class:`~.PCA`, or :class:`~.EDA`
  * a :class:`~.Mode` obtained by indexing an NMA model, e.g. ``anm[0]``
  * a :class:`~.ModeSet` obtained by slicing an NMA model, e.g. ``anm[0:10]``

Some of these functions may also accept :class:`Vector` instances as *mode*
argument.  These are noted in function documentations.

Analyze models
==============

Following functions are for calculating atomic properties from normal modes:

  * :func:`~.calcCollectivity` - degree of collectivity of a mode
  * :func:`~.calcCovariance` - covariance matrix for given modes
  * :func:`~.calcCrossCorr` - cross-correlations of fluctuations
  * :func:`~.calcFractVariance` - fraction of variance explained by a mode
  * :func:`~.calcPerturbResponse` - response to perturbations in positions
  * :func:`~.calcProjection` - projection of conformations onto modes
  * :func:`~.calcSqFlucts` - square-fluctuations
  * :func:`~.calcTempFactors` - temperature factors fitted to exp. data

Compare models
==============

Following functions are for comparing normal modes or dynamics models:

  * :func:`~.calcOverlap` - overlap (correlation) between modes
  * :func:`~.calcCumulOverlap` - cumulative overlap between modes
  * :func:`~.calcSubspaceOverlap` - overlap between normal mode subspaces
  * :func:`~.calcCovOverlap` - covariance overlap between models
  * :func:`~.printOverlapTable` - formatted overlap table printed on screen

Generate conformers
===================

Following functions can be used to generate conformers along normal modes:

  * :func:`~.deformAtoms` - deform atoms along a mode
  * :func:`~.sampleModes` - deform along random combination of a set of modes
  * :func:`~.traverseMode` - traverse a mode along both directions

Editing models
==============

Following functions can be used to reduce, slice, or extrapolate models:

  * :func:`~.sliceMode` - take a slice of the normal mode
  * :func:`~.sliceModel` - take a slice of a model
  * :func:`~.sliceVector` - take a slice of a vector
  * :func:`~.reduceModel` - reduce a model to a subset of atoms
  * :func:`~.extrapolateModel` - extrapolate a coarse-grained model to all-atoms

Parse/write data
================

Following functions are parsing or writing normal mode data:

  * :func:`~.parseArray` - numeric arrays, e.g. coordinates, eigenvectors
  * :func:`~.parseModes` - normal modes
  * :func:`~.parseNMD` - normal mode, coordinate, and atomic data for NMWiz
  * :func:`~.parseSparseMatrix` - matrix data in sparse coordinate list format
  * :func:`~.writeArray` - numeric arrays, e.g. coordinates, eigenvectors
  * :func:`~.writeModes` - normal modes
  * :func:`~.writeNMD` - normal mode, coordinate, and atomic data
  * :func:`~.writeOverlapTable` - overlap between modes in a formatted table

Save/load models
================

Dynamics objects can be efficiently saved and loaded in later Python sessions
using the following functions:

  * :func:`~.loadModel`, :func:`~.saveModel` - load/save dynamics models
  * :func:`~.loadVector`, :func:`~.saveVector` - load/save modes or vectors

Short-hand functions
====================

Following allow for performing some dynamics calculations in one function call:

  * :func:`~.calcANM` - perform ANM calculations
  * :func:`~.calcGNM` - perform GNM calculations

Plotting functions
==================

Plotting functions are called by the name of the plotted data/property
and are prefixed with "show".  Function documentations refers to the
:mod:`matplotlib.pyplot` function utilized for actual plotting.
Arguments and keyword arguments are passed to the Matplotlib functions.

  * :func:`~.showMode` - mode shape
  * :func:`~.showOverlap` - overlap between modes
  * :func:`~.showSqFlucts` - square-fluctuations
  * :func:`~.showEllipsoid` - depict projection of a normal mode space on
    another
  * :func:`~.showContactMap` - contact map based on a Kirchhoff matrix
  * :func:`~.showProjection` - projection of conformations onto normal modes
  * :func:`~.showOverlapTable` - overlaps between two models
  * :func:`~.showScaledSqFlucts` - square-fluctuations fitted to experimental
    data
  * :func:`~.showNormedSqFlucts` - normalized square-fluctuations
  * :func:`~.showCrossProjection` - project conformations onto modes from
    different models
  * :func:`~.showCrossCorr` - cross-correlations between fluctuations
    in atomic positions
  * :func:`~.showCumulOverlap` - cumulative overlap of a mode with multiple
    modes from another model
  * :func:`~.showFractVars` - fraction of variances
  * :func:`~.showCumulFractVars` - cumulative fraction of variances
  * :func:`~.resetTicks` - change ticks in a plot

Visualize modes
===============

Finally, normal modes can be visualized and animated using VMD plugin
:ref:`nmwiz`. Following functions allow for running NMWiz from within Python:

  * :func:`~.viewNMDinVMD` - run VMD and load normal mode data
  * :func:`~.getVMDpath`, :func:`~.setVMDpath` - get/set path to VMD executable

Examples
========

Results from the example :ref:`pca-xray-calculations` will be used to
illustrate class methods and functions in the module.

>>> from prody import *
>>> import matplotlib.pyplot as plt
>>> import numpy as np

>>> p38_pca = loadModel('p38_xray.pca.npz')
>>> p38_anm = loadModel('1p38.anm.npz')
>>> p38_ensemble = loadEnsemble('p38_X-ray.ens.npz')
>>> p38_structure = parsePDB('p38_ref_chain.pdb')

.. plot::
   :nofigs:
   :context:

   import matplotlib.pyplot as plt

   plt.close('all')
"""

__author__ = 'Ahmet Bakan'
__copyright__ = 'Copyright (C) 2010-2012 Ahmet Bakan'

import prody
LOGGER = prody.LOGGER
SETTINGS = prody.SETTINGS

__all__ = []

import analysis
from analysis import *
__all__.extend(analysis.__all__)

import compare
from compare import *
__all__.extend(compare.__all__)

import editing
from editing import *
__all__.extend(editing.__all__)

import functions
from functions import *
__all__.extend(functions.__all__)

import plotting
from plotting import *
__all__.extend(plotting.__all__)

import sampling
from sampling import *
__all__.extend(sampling.__all__)

import nma
from nma import *
__all__.extend(nma.__all__)

import mode
from mode import *
__all__.extend(mode.__all__)

mode.calcCollectivity = calcCollectivity

import modeset
from modeset import *
__all__.extend(modeset.__all__)

modeset.NMA = NMA

import gamma
from gamma import *
__all__.extend(gamma.__all__)

import anm
from anm import *
__all__.extend(anm.__all__)

import gnm
from gnm import *
__all__.extend(gnm.__all__)

import pca
from pca import *
__all__.extend(pca.__all__)

max_stars_repo_name: gokceneraslan/ProDy

max_issues_repo_name: gokceneraslan/ProDy

max_forks_repo_name: gokceneraslan/ProDy