text:
import os
import sys
import string
import matplotlib
matplotlib.use('Agg')
from pylab import *
import numpy as np

class MatPlotter:
    fontsize = 15
    color = 0
    colors = ["r-", "b-", "k-", "c-", "m-", "y-"]
    colors += [x + "-" for x in colors]
    colors += ["g-", "g--"]
    curFigure = []
    legendNames = []
    fontsizeLegend = 14
    legendPlace = 'lower right'
    legendborderpad = None
    legendlabelsep = None

    def __init__(self, fontsize=15):
        # self.newFigure()
        self.fontsize = fontsize
        self.fontsizeLegend = fontsize - 1

    def formatLegend(
        self,
        newFontSize=14,
        newPlace='lower right',
        borderpad=None,
        labelsep=None
    ):
        self.fontsizeLegend = newFontSize
        self.legendPlace = newPlace
        self.legendborderpad = borderpad
        self.legendlabelsep = labelsep

    def newFigure(self, plotTitle="", fsize=rcParams['figure.figsize']):
        return self.newRPCFigure(plotTitle, fsize)

    def newRPCFigure(self, plotTitle="", fsize=rcParams['figure.figsize']):
        curFigure = figure(figsize=fsize)
        self.title = title(plotTitle, fontsize=self.fontsize)
        #subplots_adjust(left=0.085, right=0.975, top=0.975, bottom=0.085)
        subplots_adjust(right=0.975, top=0.975)

        #axis('equal')
        axis([0, 1, 0, 1])
        xticklocs, xticklabels = xticks(arange(0, 1.01, 0.1))
        setp(xticklabels, size=self.fontsize)
        yticklocs, yticklabels = yticks(arange(0, 1.01, 0.1))
        setp(yticklabels, size=self.fontsize)
        self.xlabel = xlabel("1-precision")
        self.xlabel.set_size(self.fontsize + 2)
        self.ylabel = ylabel("recall")
        self.ylabel.set_size(self.fontsize + 4)
        grid()
        hold(True)

    def newFPPIFigure(self, plotTitle="", fsize=rcParams['figure.figsize']):
        curFigure = figure(figsize=fsize)
        self.title = title(plotTitle, fontsize=self.fontsize)
        subplots_adjust(left=0.085, right=0.975, top=0.975, bottom=0.085)

        #axis('equal')
        axis([0, 100, 0, 1])
        xticklocs, xticklabels = xticks(arange(0, 100.01, 0.5))
        setp(xticklabels, size=self.fontsize)
        yticklocs, yticklabels = yticks(arange(0, 1.01, 0.1))
        setp(yticklabels, size=self.fontsize)
        self.xlabel = xlabel("false positives per image")
        self.xlabel.set_size(self.fontsize + 2)
        self.ylabel = ylabel("recall")
        self.ylabel.set_size(self.fontsize + 4)
        grid()
        hold(True)

    def newFreqFigure(
        self, plotTitle="", maxX=10, maxY=10, fsize=rcParams['figure.figsize']
    ):
        curFigure = figure(figsize=fsize)
        self.title = title(plotTitle, fontsize=self.fontsize)
        subplots_adjust(left=0.085, right=0.975, top=0.975, bottom=0.1)
        #axis('equal')

        axis([0, maxX, 0, maxY])
        xticklocs, xticklabels = xticks(
            arange(0, maxX + 0.01, maxX * 1.0 / 10)
        )
        setp(xticklabels, size=self.fontsize)
        yticklocs, yticklabels = yticks(
            arange(0, maxY + 0.01, maxY * 1.0 / 10)
        )
        setp(yticklabels, size=self.fontsize)
        self.xlabel = xlabel("False positive / ground truth rect")
        self.xlabel.set_size(self.fontsize + 2)
        self.ylabel = ylabel("True positives / ground truth rect")
        self.ylabel.set_size(self.fontsize + 4)
        grid()
        hold(True)

    def newFPPWFigure(self, plotTitle="", fsize=rcParams['figure.figsize']):
        curFigure = figure(figsize=fsize)
        self.title = title(plotTitle, fontsize=self.fontsize)
        subplots_adjust(left=0.085, right=0.975, top=0.975, bottom=0.085)

        self.xlabel = xlabel("false positive per windows (FPPW)")
        self.xlabel.set_size(self.fontsize + 2)
        self.ylabel = ylabel("miss rate")
        self.ylabel.set_size(self.fontsize + 4)

        grid()
        hold(True)

    def newLogFPPIFigure(self, plotTitle="", fsize=rcParams['figure.figsize']):
        curFigure = figure(figsize=fsize)
        self.title = title(plotTitle, fontsize=self.fontsize)
        subplots_adjust(left=0.085, right=0.975, top=0.975, bottom=0.1)

        #axis('equal')

        self.xlabel = xlabel("false positives per image")
        self.xlabel.set_size(self.fontsize + 2)
        self.ylabel = ylabel("miss rate")
        self.ylabel.set_size(self.fontsize + 4)
        grid()
        hold(True)

    def loadRPCData(self, fname):
        self.filename = fname
        self.prec = []
        self.rec = []
        self.score = []
        self.fppi = []
        file = open(fname)

        precScores = []
        for i in range(1, 10, 1):
            precScores.append(100 - i * 10)

        fppiScores = []
        for i in range(0, 500, 5):
            fppiScores.append(i * 1.0 / 100.0)

        precinfo = []
        fppiinfo = []
        eerinfo = []
        logAvInfo = []

        logAvMR = []
        self.lamr = 0
        self.eer = None
        firstLine = True
        leadingZeroCount = 0

        for line in file.readlines():
            vals = line.split()
            #vals=line.split(" ")
            #for val in vals:
            #       if val=="":
            #               vals.remove(val)
            self.prec.append(1 - float(vals[0]))
            self.rec.append(float(vals[1]))
            self.score.append(float(vals[2]))

            if (len(vals) > 3):
                self.fppi.append(float(vals[3]))
                if firstLine and not float(vals[3]) == 0:
                    firstLine = False

                    lamrcount = 1
                    self.lamr = 1 - float(vals[1])

                    lowest_fppi = math.ceil(
                        math.log(float(vals[3])) / math.log(10) * 10
                    )
                    print "lowest_fppi: ", lowest_fppi

                    # MA: temporarily commented out
                    # for i in range(lowest_fppi, 1, 1):
                    #       logAvMR.append(10** (i * 1.0 / 10))

                #self.score.append(float(vals[2][:-1]))
                #print 1-self.prec[-1], self.rec[-1], self.score[-1]
            if (len(self.prec) > 1):
                diff = (1 - self.prec[-1] - self.rec[-1]
                       ) * (1 - self.prec[-2] - self.rec[-2])
                if (diff < 0):
                    eerinfo.append(
                        "EER between: %.03f and %.03f\tScore:%f" %
                        (self.rec[-1], self.rec[-2], self.score[-1])
                    )
                    self.eer = (self.rec[-1] + self.rec[-2]) * 0.5
                if (diff == 0 and 1 - self.prec[-1] - self.rec[-1] == 0):
                    eerinfo.append(
                        "EER: %.03f\tScore:%f" %
                        (self.rec[-1], self.score[-1])
                    )
                    self.eer = self.rec[-1]

            #Remove already passed precision
            if (len(precScores) > 0 and (float(vals[0])) < precScores[0] / 100.0):
                precinfo.append(
                    "%d percent precision score: %f, recall: %.03f" %
                    (precScores[0], float(vals[2]), float(vals[1]))
                )
                while (
                    len(precScores) > 0 and
                    precScores[0] / 100.0 > float(vals[0])
                ):
                    precScores.pop(0)

            #Remove already passed precision
            if (len(vals) > 3):
                if (len(fppiScores) > 0 and (float(vals[3])) > fppiScores[0]):
                    fppiinfo.append(
                        "%f fppi score: %f, recall: %.03f" %
                        (fppiScores[0], float(vals[2]), float(vals[1]))
                    )
                    while (
                        len(fppiScores) > 0 and fppiScores[0] < float(vals[3])
                    ):
                        fppiScores.pop(0)

                if (len(logAvMR) > 0 and (float(vals[3])) > logAvMR[0]):
                    while (len(logAvMR) > 0 and logAvMR[0] < float(vals[3])):
                        logAvInfo.append(
                            "%f fppi, miss rate: %.03f, score: %f" %
                            (logAvMR[0], 1 - float(vals[1]), float(vals[2]))
                        )
                        self.lamr += 1 - float(vals[1])
                        lamrcount += 1
                        logAvMR.pop(0)

                lastMR = 1 - float(vals[1])

        if (len(vals) > 3):
            for i in logAvMR:
                logAvInfo.append(
                    "%f fppi, miss rate: %.03f, extended" % (i, lastMR)
                )
                self.lamr += lastMR
                lamrcount += 1

        for i in precinfo:
            print i
        print
        for i in fppiinfo:
            print i
        print
        for i in eerinfo:
            print i
        print
        print "Recall at first false positive: %.03f" % self.rec[0]
        if (len(vals) > 3):
            print
            for i in logAvInfo:
                print i
            self.lamr = self.lamr * 1.0 / lamrcount
            print "Log average miss rate in [10^%.01f, 10^0]: %.03f" % (
                lowest_fppi / 10.0, self.lamr
            )

        print
        print
        file.close()

    def loadFreqData(self, fname):
        self.filename = fname
        self.prec = []
        self.rec = []
        self.score = []
        file = open(fname)

        for line in file.readlines():
            vals = line.split()

            self.prec.append(float(vals[0]))
            self.rec.append(float(vals[1]))
            self.score.append(float(vals[2]))

        file.close()

    def loadFPPWData(self, fname):
        self.loadFreqData(fname)

    def finishPlot(self, axlimits=[0, 1.0, 0, 1.0]):
        # MA:
        #self.legend = legend(self.legendNames, self.legendPlace, pad = self.legendborderpad, labelsep = self.legendlabelsep)
        self.legend = legend(self.legendNames, self.legendPlace)

        lstrings = self.legend.get_texts()
        setp(lstrings, fontsize=self.fontsizeLegend)
        #line= plot( [1 - axlimits[0], 0], [axlimits[3], 1 - axlimits[3] ] , 'k')
        line = plot([1, 0], [0, 1], 'k')

    def finishFreqPlot(self):
        self.legend = legend(
            self.legendNames,
            self.legendPlace,
            pad=self.legendborderpad,
            labelsep=self.legendlabelsep
        )
        lstrings = self.legend.get_texts()
        setp(lstrings, fontsize=self.fontsizeLegend)

    def show(self, plotEER=True, axlimits=[0, 1.0, 0, 1.0]):
        if (plotEER):
            self.finishPlot(axlimits)
            axis(axlimits)
        else:
            self.finishFreqPlot()

        show()

    def saveCurrentFigure(self, plotEER, filename, axlimits=[0, 1.0, 0, 1.0]):
        if (plotEER):
            self.finishPlot(axlimits)
            axis(axlimits)
        else:
            self.finishFreqPlot()

        print "Saving: " + filename
        savefig(filename)

    def plotRFP(self, numImages, fname, line="r-"):
        print 'NOT YET IMPLEMENTED'

    def plotRPC(
        self,
        fname,
        descr="line",
        style="-1",
        axlimits=[0, 1.0, 0, 1.0],
        linewidth=2,
        dashstyle=[],
        addEER=False
    ):
        self.loadRPCData(fname)

        #axis(axlimits);
        if (style == "-1"):
            if dashstyle != []:
                line = plot(
                    self.prec,
                    self.rec,
                    self.colors[self.color],
                    dashes=dashstyle
                )
            else:
                line = plot(self.prec, self.rec, self.colors[self.color])
            self.color = self.color + 1
            self.color = self.color % len(self.colors)
        else:
            if dashstyle != []:
                line = plot(self.prec, self.rec, style, dashes=dashstyle)
            else:
                line = plot(self.prec, self.rec, style)

        axis(axlimits)

        if addEER and self.eer != None:
            descr += " (%.01f%%)" % (self.eer * 100)

        setp(line, 'linewidth', linewidth)
        self.legendNames = self.legendNames + [descr]

    def plotFPPI(
        self,
        fname,
        descr="line",
        style="-1",
        axlimits=[0, 2, 0, 1],
        linewidth=2,
        dashstyle=[]
    ):
        self.loadRPCData(fname)

        if (style == "-1"):
            if dashstyle != []:
                line = plot(
                    self.fppi,
                    self.rec,
                    self.colors[self.color],
                    dashes=dashstyle
                )
            else:
                line = plot(self.fppi, self.rec, self.colors[self.color])
            self.color = self.color + 1
            self.color = self.color % len(self.colors)
        else:
            if dashstyle != []:
                line = plot(self.fppi, self.rec, style, dashes=dashstyle)
            else:
                line = plot(self.fppi, self.rec, style)

        axis(axlimits)

        setp(line, 'linewidth', linewidth)
        self.legendNames = self.legendNames + [descr]

    def plotFreq(
        self, fname, descr="line", style="-1", linewidth=2, dashstyle=[]
    ):
        self.loadFreqData(fname)
        if (style == "-1"):
            if dashstyle != []:
                line = plot(
                    self.prec,
                    self.rec,
                    self.colors[self.color],
                    dashes=dashstyle
                )
            else:
                line = plot(self.prec, self.rec, self.colors[self.color])
            self.color = self.color + 1
            self.color = self.color % len(self.colors)
        else:
            if dashstyle != []:
                line = plot(self.prec, self.rec, style, dashes=dashstyle)
            else:
                line = plot(self.prec, self.rec, style)

        setp(line, 'linewidth', linewidth)
        self.legendNames = self.legendNames + [descr]

    def plotFPPW(
        self,
        fname,
        descr="line",
        style="-1",
        axlimits=[5e-6, 1e0, 1e-2, 0.5],
        linewidth=2,
        dashstyle=[]
    ):
        self.loadFPPWData(fname)
        if (style == "-1"):
            if dashstyle != []:
                line = loglog(
                    self.prec,
                    self.rec,
                    self.colors[self.color],
                    dashes=dashstyle
                )
            else:
                line = loglog(self.prec, self.rec, self.colors[self.color])
            self.color = self.color + 1
            self.color = self.color % len(self.colors)
        else:
            if dashstyle != []:
                line = loglog(self.prec, self.rec, style, dashes=dashstyle)
            else:
                line = loglog(self.prec, self.rec, style)

        xticklocs, xticklabels = xticks([1e-5, 1e-4, 1e-3, 1e-2, 1e-1, 1e0])
        setp(xticklabels, size=self.fontsize)
        yticklocs, yticklabels = yticks(
            array(
                [
                    0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.1,
                    0.2, 0.3, 0.4, 0.5
                ]
            ), (
                "0.01", "0.02", "0.03", "0.04", "0.05", "0.06", "0.07", "0.08",
                "0.09", "0.1", "0.2", "0.3", "0.4", "0.5"
            )
        )
        setp(yticklabels, size=self.fontsize)

        axis(axlimits)

        gca().yaxis.grid(True, 'minor')
        setp(line, 'linewidth', linewidth)

        self.legendNames = self.legendNames + [descr]

    def plotLogFPPI(
        self,
        fname,
        descr="line",
        style="-1",
        axlimits=[5e-3, 1e1, 1e-1, 1],
        linewidth=2,
        dashstyle=[],
        addlamr=False
    ):
        self.loadRPCData(fname)
        if (style == "-1"):
            if dashstyle != []:
                line = loglog(
                    self.fppi, [1 - x for x in self.rec],
                    self.colors[self.color],
                    dashes=dashstyle
                )
            else:
                line = loglog(
                    self.fppi, [1 - x
                                for x in self.rec], self.colors[self.color]
                )

            self.color = (self.color + 1) % len(self.colors)
        else:
            if dashstyle != []:
                line = loglog(
                    self.fppi, [1 - x for x in self.rec],
                    style,
                    dashes=dashstyle
                )
            else:
                line = loglog(self.fppi, [1 - x for x in self.rec], style)

        gca().yaxis.grid(True, 'minor')

        m = min(self.fppi)
        lax = axlimits[0]
        for i in self.fppi:
            if (i != m):
                lax = math.floor(log(i) / math.log(10))
                leftlabel = math.pow(10, lax)
                break

        m = max(self.fppi)
        rightlabel = math.pow(10, math.ceil(log(m) / math.log(10))) + 0.01

        k = leftlabel
        ticks = [k]
        while k < rightlabel:
            k = k * 10
            ticks.append(k)

        xticklocs, xticklabels = xticks(ticks)
        setp(xticklabels, size=self.fontsize)
        yticklocs, yticklabels = yticks(
            arange(0.1, 1.01, 0.1), (
                "0.1", "0.2", "0.3", "0.4", "0.5", "0.6", "0.7", "0.8", "0.9",
                "1.0"
            )
        )
        setp(yticklabels, size=self.fontsize)

        axlimits[0] = lax
        axis(axlimits)

        setp(line, 'linewidth', linewidth)

        if addlamr:
            descr += " (%.01f%%)" % (self.lamr * 100)

        self.legendNames = self.legendNames + [descr]

max_stars_repo_name: mdcatapult/deepfigures-open

max_issues_repo_name: mdcatapult/deepfigures-open

max_forks_repo_name: mdcatapult/deepfigures-open


text:
# Import all dependencies.
# Note: Please import matplotlib before importing any psyneulink dependencies.
import numpy as np

from matplotlib import pyplot as plt
# from scipy.special import erfinv # need to import this to make us of the UniformToNormalDist function.
import psyneulink as pnl

# --------------------------------- Global Variables ----------------------------------------
# Now, we set the global variables, weights and initial values as in the paper.
# WATCH OUT !!! In the paper the weight "Mutual inhibition among response units" is not defined, but needs to be set to
# 0 in order to reproduce the paper.

SD = 0.15       # noise determined by standard deviation (SD)
a = 0.50        # Parameter describing shape of the FitzHugh–Nagumo cubic nullcline for the fast excitation variable v
d = 0.5         # Uncorrelated Activity
k = 1.5         # Scaling factor for transforming NE release (u ) to gain (g ) on potentiated units
G = 0.5         # Base level of gain applied to decision and response units
dt = 0.02       # time step size
C = 0.90        # LC coherence (see Gilzenrat et al. (2002) on more details on LC coherence

initial_hv = 0.07                     # Initial value for h(v)
initial_w = 0.14                      # initial value u
initial_v = (initial_hv - (1 - C) * d) / C  # get initial v from initial h(v)

# Weights:
inpwt = 1.5       # inpwt (Input to decision layer)
crswt = 1 / 3       # crswt (Crosstalk input to decision layer)
inhwt = 1.0       # inhwt (Mutual inhibition among decision units)
respinhwt = 0     # respinhwt (Mutual inhibition among response units)  !!! WATCH OUT: this parameter is not mentioned
# in the original paper, most likely since it was set 0
decwt = 3.5       # decwt (Target decision unit to response unit)
selfdwt = 2.5     # selfdwt (Self recurrent conn. for each decision unit)
selfrwt = 2.0     # selfrwt (Self recurrent conn. for response unit)
lcwt = 0.3        # lcwt (Target decision unit to LC)
decbias = 1.75    # decbias (Bias input to decision units)
respbias = 1.75   # respbias (Bias input to response units)
tau_v = 0.05    # Time constant for fast LC excitation variable v
tau_u = 5.00    # Time constant for slow LC recovery variable (‘NE release’) u
trials = 1100   # number of trials to reproduce Figure 3 from Nieuwenhuis et al. (2005)

# Setting seed (if noise is applied)
# np.random.seed(22)    # Set noise with seed generator to 22 to use the UniformToNormalDist function to have a python
# seed compatible with MATLAB random seed generator (rsg=22)
# Please see https://github.com/jonasrauber/randn-matlab-python for further documentation

# Create mechanisms ---------------------------------------------------------------------------------------------------

# Input Layer --- [ Target 1, Target 2, Distractor ]

# First, we create the 3 layers of the behavioral network, i.e. INPUT LAYER, DECISION LAYER, and RESPONSE LAYER.
input_layer = pnl.TransferMechanism(
    size=3,                      # Number of units in input layer
    initial_value=[[0.0, 0.0, 0.0]],     # Initial input values
    name='INPUT LAYER'                  # Define the name of the layer; this is optional,
)                 # but will help you to overview your model later on

# Create Decision Layer  --- [ Target 1, Target 2, Distractor ]
decision_layer = pnl.LCA(
    size=3,                            # Number of units in input layer
    initial_value=[[0.0, 0.0, 0.0]],    # Initial input values
    time_step_size=dt,                 # Integration step size
    leak=-1.0,                         # Sets off diagonals to negative values
    self_excitation=selfdwt,           # Set diagonals to self excitate
    competition=inhwt,                 # Set off diagonals to inhibit
    function=pnl.Logistic(bias=decbias),   # Set the Logistic function with bias = decbias
    # noise=pnl.UniformToNormalDist(standard_dev = SD).function, # The UniformToNormalDist function will
    integrator_mode=True,               # set the noise with a seed generator that is compatible with
    name='DECISION LAYER'               # MATLAB random seed generator 22 (rsg=22)
)

# decision_layer.set_log_conditions('RESULT')  # Log RESULT of the decision layer
decision_layer.set_log_conditions('value')  # Log value of the decision layer

for output_state in decision_layer.output_states:
    output_state.value *= 0.0                                       # Set initial output values for decision layer to 0

# Create Response Layer  --- [ Target1, Target2 ]
response_layer = pnl.LCA(
    size=2,                                        # Number of units in input layer
    initial_value=[[0.0, 0.0]],                    # Initial input values
    time_step_size=dt,                             # Integration step size
    leak=-1.0,                                     # Sets off diagonals to negative values
    self_excitation=selfrwt,                       # Set diagonals to self excitate
    competition=respinhwt,                         # Set off diagonals to inhibit
    function=pnl.Logistic(bias=respbias),          # Set the Logistic function with bias = decbias
    # noise=pnl.UniformToNormalDist(standard_dev = SD).function,
    integrator_mode=True,
    name='RESPONSE LAYER'
)

response_layer.set_log_conditions('RESULT')     # Log RESULT of the response layer
for output_state in response_layer.output_states:
    output_state.value *= 0.0                   # Set initial output values for response layer to 0

# Connect mechanisms --------------------------------------------------------------------------------------------------
# Weight matrix from Input Layer --> Decision Layer
input_weights = np.array([
    [inpwt, crswt, crswt],                    # Input weights are diagonals, cross weights are off diagonals
    [crswt, inpwt, crswt],
    [crswt, crswt, inpwt]
])

# Weight matrix from Decision Layer --> Response Layer
output_weights = np.array([
    [decwt, 0.0],  # Weight from T1 and T2 but not distractor unit (row 3 set to all zeros) to response layer
    [0.0, decwt],  # Need a 3 by 2 matrix, to project from decision layer with 3 units to response layer with 2 units
    [0.0, 0.0]
])

decision_process = pnl.Process(
    pathway=[
        input_layer,
        input_weights,
        decision_layer,
        output_weights,
        response_layer
    ],
    name='DECISION PROCESS'
)

# Abstracted LC to modulate gain --------------------------------------------------------------------

# This LCControlMechanism modulates gain.
LC = pnl.LCControlMechanism(
    integration_method="EULER",       # We set the integration method to Euler like in the paper
    threshold_FHN=a,                  # Here we use the Euler method for integration and we want to set the parameters,
    uncorrelated_activity_FHN=d,      # for the FitzHugh–Nagumo system.
    time_step_size_FHN=dt,
    mode_FHN=C,
    time_constant_v_FHN=tau_v,
    time_constant_w_FHN=tau_u,
    a_v_FHN=-1.0,
    b_v_FHN=1.0,
    c_v_FHN=1.0,
    d_v_FHN=0.0,
    e_v_FHN=-1.0,
    f_v_FHN=1.0,
    a_w_FHN=1.0,
    b_w_FHN=-1.0,
    c_w_FHN=0.0,
    t_0_FHN=0.0,
    base_level_gain=G,                # Additionally, we set the parameters k and G to compute the gain equation
    scaling_factor_gain=k,
    initial_v_FHN=initial_v,          # Initialize v
    initial_w_FHN=initial_w,          # Initialize w
    objective_mechanism=pnl.ObjectiveMechanism(
        function=pnl.Linear,
        monitored_output_states=[(
            decision_layer,  # Project output of T1 and T2 but not distractor from decision layer to LC
            np.array([[lcwt], [lcwt], [0.0]])
        )],
        name='Combine values'
    ),
    modulated_mechanisms=[decision_layer, response_layer],  # Modulate gain of decision & response layers
    name='LC'
)

# Log value of LC
LC.set_log_conditions('value')

# Set initial gain to G + k*initial_w, when the System runs the very first time,
# since the decison layer executes before the LC and hence needs one initial gain value to start with.
for output_state in LC.output_states:
    output_state.value *= G + k * initial_w

# Now, we specify the processes of the System, which in this case is just the decision_process
task = pnl.System(processes=[decision_process])

# Create Stimulus -----------------------------------------------------------------------------------------------------

# In the paper, each period has 100 time steps, so we will create 11 time periods.
# As described in the paper in figure 3, during the first 3 time periods input to distractor units is fixed to 1.
# Then T1 gets turned on during time period 4 with an input of 1.
# T2 gets turns on with some lag from T1 onset on, in this example we turn T2 on with Lag 2 and an input of 1
# Between T1 and T2 and after T2 the distractor unit is on.
# We create one array with 3 numbers, one for each input unit and repeat this array 100 times for one time period
# We do this 11 times. T1 is on for time4, T2 is on for time7 to model Lag3
num_time_steps = 100  # Each stimulus is presented for two units of time which is equivalent to 100 time steps
stimulus_T1 = np.repeat(np.array([[0, 0, 1]]), num_time_steps, axis=0)
stimulus_T2 = np.repeat(np.array([[0, 0, 1]]), num_time_steps, axis=0)
stimulus_T3 = np.repeat(np.array([[0, 0, 1]]), num_time_steps, axis=0)
stimulus_T4 = np.repeat(np.array([[1, 0, 0]]), num_time_steps, axis=0)    # Turn T1 on
stimulus_T5 = np.repeat(np.array([[0, 0, 1]]), num_time_steps, axis=0)
stimulus_T6 = np.repeat(np.array([[0, 1, 0]]), num_time_steps, axis=0)    # Turn T2 on --> example for Lag 2
stimulus_T7 = np.repeat(np.array([[0, 0, 1]]), num_time_steps, axis=0)
stimulus_T8 = np.repeat(np.array([[0, 0, 1]]), num_time_steps, axis=0)
stimulus_T9 = np.repeat(np.array([[0, 0, 1]]), num_time_steps, axis=0)
stimulus_T10 = np.repeat(np.array([[0, 0, 1]]), num_time_steps, axis=0)
stimulus_T11 = np.repeat(np.array([[0, 0, 1]]), num_time_steps, axis=0)

# Concatenate the 11 arrays to one array with 1100 rows and 3 columns.
time = np.concatenate((stimulus_T1, stimulus_T2, stimulus_T3, stimulus_T4, stimulus_T5, stimulus_T6,
                       stimulus_T7, stimulus_T8, stimulus_T9, stimulus_T10, stimulus_T11), axis=0)

# assign inputs to input_layer (Origin Mechanism) for each trial
stim_list_dict = {input_layer: time}

# show the system
task.show()
# run the system
task.run(stim_list_dict, num_trials=trials)

# This displays a diagram of the System
# task.show_graph()

LC_results = LC.log.nparray()        # get logged results
LC_results_w = np.zeros([trials])    # get LC_results_w
for i in range(trials):
    LC_results_w[i] = LC_results[3][i+1][3][0]
LC_results_v = np.zeros([trials])    # get LC_results_v
for i in range(trials):
    LC_results_v[i] = LC_results[3][i+1][2][0]

def h_v(v, C, d):                   # Compute h(v)
    return C * v + (1 - C) * d
LC_results_hv = np.zeros([trials])    # get LC_results_hv
for i in range(trials):
    LC_results_hv[i] = h_v(LC_results_v[i],C,d)

#Plot the Figure 3 from the paper
t = np.linspace(0,trials,trials)            # Create array for x axis with same length then LC_results_v
fig = plt.figure()                          # Instantiate figure
ax = plt.gca()                              # Get current axis for plotting
ax2 = ax.twinx()                            # Create twin axis with a different y-axis on the right side of the figure
ax.plot(t, LC_results_hv, label="h(v)")      # Plot h(v)
ax2.plot(t, LC_results_w, label="w", color = 'red') # Plot w
h1, l1 = ax.get_legend_handles_labels()
h2, l2 = ax2.get_legend_handles_labels()
ax.legend(h1 + h2, l1 + l2, loc=2)          # Create legend on one side
ax.set_xlabel('Time (ms)')                  # Set x axis lable
ax.set_ylabel('LC Activity')                # Set left y axis label
ax2.set_ylabel('NE Output')                 # Set right y axis label
plt.title('Nieuwenhaus 2005 PsyNeuLink Lag 2 without noise', fontweight='bold') # Set title
ax.set_ylim((-0.2,1.0))                     # Set left y axis limits
ax2.set_ylim((0.0, 0.4))                    # Set right y axis limits
plt.show()

ext: py

max_stars_repo_name: mihaic/psyneulink

max_issues_repo_name: mihaic/psyneulink

max_forks_repo_name: mihaic/psyneulink

reason: import numpy,from scipy


text:
""" Flag data that lies outside the zone of Marie's paper """

import numpy as np
import sys
sys.path.append('..')
from mass_age_functions import *

id_all = np.load("id_all.npz")['arr_0']
nobj = len(id_all)

flags = np.ones(nobj) # start with all flagged

# [M/H] > -0.8
feh = np.load("feh_all.npz")['arr_0']
good_feh = feh > -0.8

# 4000 < Teff < 5000
teff = np.load("teff_all.npz")['arr_0']
good_teff = np.logical_and(teff>4000, teff<5000)

# 1.8 < logg < 3.3
logg = np.load("logg_all.npz")['arr_0']
good_logg = np.logical_and(logg>1.8, logg<3.3)

# -0.25 < [C/M] < 0.15
cm = np.load("cm_all.npz")['arr_0']
good_cm = np.logical_and(cm > -0.25, cm < 0.15)

# -0.1 < [N/M] < 0.45
nm = np.load("nm_all.npz")['arr_0']
good_nm = np.logical_and(nm > -0.1, nm < 0.45)

# -0.1 < [(C+N)/M] < 0.15
cplusm = calc_sum(feh, cm, nm)
good_cplusm = np.logical_and(cplusm > -0.1, cplusm < 0.15)

# -0.6 < [C/N] < 0.2
cn = calc_cn(cm, nm)
good_cn = np.logical_and(cn > -0.6, cn < 0.2)

keep = np.logical_and.reduce(
        (good_feh, good_teff, good_logg,
        good_cm, good_nm, good_cplusm, good_cn))

flags[keep] = 0

np.savez("flags.npz", flags)

max_stars_repo_name: aceilers/TheCannon

max_issues_repo_name: aceilers/TheCannon

max_forks_repo_name: aceilers/TheCannon

reason: import numpy


text:
import torch
import numpy as np
from scipy.spatial.transform import Rotation as ROT
import torch.utils.data as data
import os
import csv
import random
from PIL import Image
class KITTIDataLoader(data.Dataset):
    def __init__(self):
        super(KITTIDataLoader, self).__init__()

    def initialize(self, opt):
        self.opt = opt
        self.dataroot = './datasets/dataset_kitti'

        self.opt.bound = 5
        with open(os.path.join(self.dataroot, 'id_train.txt'), 'r') as fp:
            self.ids_train = [s.strip() for s in fp.readlines() if s]

        self.ids = self.ids_train
        self.dataset_size = int(len(self.ids))// (opt.bound*2)

        self.pose_dict = {}
        pose_path = os.path.join(self.dataroot, 'poses.txt')
        with open(pose_path) as csv_file:
            csv_reader = csv.reader(csv_file, delimiter=' ')
            for row in csv_reader:
                id = row[0]
                self.pose_dict[id] = []
                for col in row[1:-1]:
                    self.pose_dict[id].append(float(col))
                self.pose_dict[id] = np.array(self.pose_dict[id])

    def __getitem__(self, index):
        id = self.ids[index]
        id_num = int(id.split('_')[-1])
        while True:
            delta = random.choice([x for x in range(-self.opt.bound, self.opt.bound+1) if x != 0] )
            id_target = id.split('_')[0] +'_' + str(id_num + delta).zfill(len(id.split('_')[-1]))
            if id_target in self.pose_dict.keys(): break

        B = self.load_image(id) / 255. * 2 - 1
        B = torch.from_numpy(B.astype(np.float32)).permute((2,0,1))
        A = self.load_image(id_target) / 255. * 2 - 1
        A = torch.from_numpy(A.astype(np.float32)).permute((2,0,1))

        poseB = self.pose_dict[id]
        poseA = self.pose_dict[id_target]
        TB = poseB[3:].reshape(3, 1)
        RB = ROT.from_euler('xyz',poseB[0:3]).as_dcm()
        TA = poseA[3:].reshape(3, 1)
        RA = ROT.from_euler('xyz',poseA[0:3]).as_dcm()
        T = RA.T.dot(TB-TA)/50.

        mat = np.block(
            [ [RA.T@RB, T],
              [np.zeros((1,3)), 1] ] )

        return {'A': A, 'B': B, 'RT': mat.astype(np.float32)}

    def load_image(self, id):
        image_path = os.path.join(self.dataroot, 'images', id + '.png')
        image = np.asarray(Image.open(image_path).convert('RGB'))
        return image

    def __len__(self):
        return self.dataset_size

    def name(self):
        return 'KITTIDataLoader'


text:
#!/usr/bin/env python

# Trying to generate samples from the conditional model.

# NOTE: My convention of indexing arrays from 0 does not agree with
# the convention in Matt's notes.

# Daniel Klein, 9/14/2011

import sys
import numpy as np
from scipy.io import savemat
from random import random

from utility import log_weighted_sample, window_permutations

# Parameters
params = {'N': 5,
          'T': 500000,
          'L': 2,
          'Delta': 5,
          'theta_method': ('sparse_unique', {'p': 1.0, 'scale': 3.0}),
          # 'theta_method': ('cascade_2', {'strength': 8.0, 'decay': 0.8}),
          # 'S_method': ('random_uniform', {'p_min': 0.05, 'p_max': 0.2}),
          'S_method': ('random_periodic', {'baseline': 0.005, 'scale': 0.002})}
if not params['T'] % params['Delta'] == 0:
    print 'Error: T must be a multiple of Delta'
    sys.exit()
params['M'] = params['T'] / params['Delta']

# Generate thetas
theta = np.zeros((params['N'], params['N'], params['L']))
method_name, method_params = params['theta_method']
if method_name == 'sparse_unique':
    for i in range(params['N']):
        if random() < method_params['p']:
            j = np.random.randint(0, params['N'])
            l = np.random.randint(0, params['L'])
            theta[i,j,l] = np.random.normal(0, method_params['scale'])
if method_name == 'cascade_2':
    for i in range(params['N']-1):
        theta[i,(i+1),1] = method_params['strength']*(method_params['decay']**i)
    theta[params['N']-1,0,0] = 0.5 * method_params['strength']

# Generate S (calling it "windows" in code)
print 'Generating window permutations'
windows = []
method_name, method_params = params['S_method']
if method_name in ['random_uniform', 'random_periodic']:
    if method_name == 'random_periodic':
        periods = params['M'] * np.random.normal(1, 1, params['N'])
        phases = np.random.uniform(0.0, 2.0*np.pi, params['N'])
    for k in range(params['M']):
        window = []
        if method_name == 'random_uniform':
            p = np.random.uniform(method_params['p_min'],method_params['p_max'])
            w_raw = np.random.binomial(1, p, (params['N'],params['Delta']))
        if method_name == 'random_periodic':
            w_raw = np.empty((params['N'], params['Delta']))
            for i in range(params['N']):
                p = (method_params['baseline'] +
                     method_params['scale'] * np.sin(periods[i]*k + phases[i]))
                w_raw[i,:] = np.random.binomial(1, p, (1, params['Delta']))
        w = np.array(w_raw, dtype='uint8')
        windows.append(window_permutations(w))
n_w = map(len, windows)

# Tabulate log-potential functions: h_1, h_2, ..., h_M
print 'Tabulating log-potential functions'
h = [np.empty(n_w[0])]
s_padded = np.zeros((params['N'],2*params['Delta']), dtype='uint8')
hits = np.zeros((params['N'],params['N'],params['L']), dtype='int32')
for w, s in enumerate(windows[0]):
    s_padded[:,params['Delta']:(2*params['Delta'])] = s
    hits[:,:,:] = 0
    for l in range(params['L']):
        tmin, tmax = params['Delta'] - (l+1), 2*params['Delta'] - (l+1)
        s_lagged = s_padded[:,tmin:tmax]
        hits[:,:,l] = np.tensordot(s_lagged, s, axes = (1,1))
    h[0][w] = np.sum(theta * hits)
for k in range(1, params['M']):
    h.append(np.empty((n_w[k-1], n_w[k])))
    for w_prev, s_prev in enumerate(windows[k-1]):
        s_padded[:,0:params['Delta']] = s_prev
        for w, s in enumerate(windows[k]):
            s_padded[:,params['Delta']:(2*params['Delta'])] = s
            hits[:,:,:] = 0
            for l in range(params['L']):
                tmin, tmax = params['Delta'] - (l+1), 2*params['Delta'] - (l+1)
                s_lagged = s_padded[:,tmin:tmax]
                hits[:,:,l] = np.tensordot(s_lagged, s, axes = (1,1))
            h[k][w_prev,w] = np.sum(theta * hits)

# Run the backward algorithm
print 'Running backward algorithm'
b = [np.zeros(n_w[params['M']-1])]
for k in range(params['M']-1, 0, -1):
    b = [np.empty(n_w[k-1])] + b
    for w_prev in range(n_w[k-1]):
        for w in range(n_w[k]):
            b[0][w_prev] = np.logaddexp.reduce(h[k][w_prev,:] + b[1])
b = [None] + b

# Sample (binned) spike trains by sampling permutations of jitter windows
print 'Sampling'
x = np.empty((params['N'], params['T']), dtype='uint8')
w_samp = log_weighted_sample(h[0] + b[1])
x[:,0:params['Delta']] = windows[0][w_samp]
for k in range(1, params['M']):
    w_samp = log_weighted_sample(h[k][w_samp,:] + b[k+1])
    x[:,(k*params['Delta']):((k+1)*params['Delta'])] = windows[k][w_samp]

# Write sample to file
savemat('sample.mat', {'theta': theta, 'sample': np.array(x, dtype='float32')},
        oned_as = 'column')

# Output
print 'Parameters'
for param in params:
    print '%s: %s' % (param, str(params[param]))
print

print 'Theta'
print theta
print

print 'x'
print x
print

print 'log_kappa'
print np.logaddexp.reduce(h[0] + b[1])


text:
import numpy as np
import tensorflow as tf

from .globals import GLOBAL_TYPE

tf.keras.backend.set_floatx(GLOBAL_TYPE)

def class_sums_from_generator(n_classes, training_generator, generator_length=-1):
    """Return an array of class counts in a dataset

    training_generator: when iterated for num_training_batches
        (or completion, if non-repeating), returns the whole
        set in batches of (x, labs)
    generator_length: break the generator iteration, if generator
        repeats (e.g. is infinite)

    Returns:
        class_sums (np.array): integer sums of each class present
            in the generator
    """
    class_sums = np.zeros((n_classes,), dtype=np.int64)
    for b, (_, label) in enumerate(training_generator):
        classes, _,  counts = tf.unique_with_counts(
            tf.reshape(tf.argmax(label, axis=-1), [-1])
        )
        for i in range(len(classes)):
            class_sums[classes[i]] += counts[i]
        if b == generator_length:
            break
    return class_sums

# NOTE: perhaps instead of for loop we could be gathering
# E.g. adjacency[label, class(rolled)] = counts(rolled)
def adjacency_from_generator(n_classes, training_generator, generator_length=-1):
    """Calculate class adjencies of dataset

    generator dim 4 req, (batch, x, y, class)

    Returns:
        adjacency: np.arrary shape (n_classes, n_classes)
            A matrix of adjacencies, per class

    """

    # TODO axis and shift aren't generalisable: 1 or -1
    def update_adjacency(
        adj_matrix, adj_to_class,
        flattened_labels,
        axis, shift
    ):
        """
        Count everything adjacent (defined by axis, shift)
        to instances of the class "adj_to_class"
        in a given batch of flattened_labels

        Add to the adj_matrix[adj_to_class] += counts
        """
        def crop_out(arr, shift, axis):
            """Crop out the edge of an arr

            Crop an edge in a given dimension by slicing array
            from that axis. Has effect of shifting elements by
            a given shift in given axis, with no roll-around.
            """
            if len(arr.shape) != 3:
                raise NotImplementedError(
                    f"Dim must be 3 not {arr.shape}"
                )
            if axis == 1:
                if shift > 0:
                    cropped = arr[:, :-shift, :]
                else:
                    cropped = arr[:, -shift:, :]
            elif axis == 2:
                if shift > 0:
                    cropped = arr[:, :, :-shift]
                else:
                    cropped = arr[:, :, -shift:]
            return cropped
        adjacent_classes = tf.where(
            # Redefine labels having removed axis that will
            # be rolled-around by the adjacency-finding shift
            crop_out(flattened_labels, -shift, axis) == adj_to_class,
            # Shift labels to put adjacencies in-place of the
            # original label cell
            crop_out(flattened_labels, shift, axis),
            # Blank-out the cells that we are not counting from
            # (e.g. that are not adjacent to adj_to_class)
            -1,
        )
        classes, _, counts = tf.unique_with_counts(
            tf.reshape(adjacent_classes, [-1])
        )
        # TODO is this a gather?
        for i in range(len(classes)):
            # -1 used to indicate non-adjacent cells
            if classes[i] >= 0:
                adj_matrix[c, classes[i]] += counts[i]

    adjacency = np.zeros((n_classes, n_classes))

    for b, (_, label) in enumerate(training_generator):
        assert len(label.shape) == 4
        flat_label = tf.argmax(label, axis=-1)
        for c in range(n_classes):
            for counting, ax, shf in (
                ("above", 1, 1), ("below", 1, -1),
                ("left", 2, 1), ("right", 2, -1),
            ):
                update_adjacency(
                    adjacency, c,
                    flat_label,
                    ax, shf,
                )

        if b == generator_length:
            break

    return adjacency

max_stars_repo_name: j-bernardi/bayesian-label-smoothing

max_issues_repo_name: j-bernardi/bayesian-label-smoothing

max_forks_repo_name: j-bernardi/bayesian-label-smoothing


text:
from collections import defaultdict

import numpy as np
from .recommendation_algorithm import (RecommendationAlgorithm,
                                       register_algorithm)

@register_algorithm('xals')
class MangakiXALS(RecommendationAlgorithm):
    def __init__(self, nb_components=10, nb_iterations=10, lambda_=0.1):
        super().__init__()
        self.M = None
        self.U = None
        self.VT = None
        self.nb_components = nb_components
        self.nb_iterations = nb_iterations
        self.lambda_ = lambda_

    @property
    def is_serializable(self):
        return True

    def make_matrix(self, X, y):
        matrix = defaultdict(dict)
        means = np.zeros((self.nb_users,))
        for (user, work), rating in zip(X, y):
            matrix[user][work] = rating
            means[user] += rating
        for user in matrix:
            means[user] /= len(matrix[user])
        for (user, work) in X:
            matrix[user][work] -= means[user]
        return matrix, means

    def fit_user(self, user, matrix):
        Ru = np.array(list(matrix[user].values()), ndmin=2).T
        Vu = self.VT[:, list(matrix[user].keys())]
        Gu = self.lambda_ * len(matrix[user]) * np.eye(self.nb_components +
                                                       self.nb_tags)
        self.U[[user], :] = np.linalg.solve(Vu.dot(Vu.T) + Gu, Vu.dot(Ru)).T

    def fit_work(self, work, matrixT):
        Ri = np.array(list(matrixT[work].values()), ndmin=2).T
        Ui = self.U[list(matrixT[work].keys()), :].T
        Gi = self.lambda_ * len(matrixT[work]) * np.eye(self.nb_components +
                                                        self.nb_tags)
        Pi = Ui[-self.nb_tags:, :]
        Ti = self.T[work]
        newV = np.linalg.solve(Ui.dot(Ui.T) + Gi, Ui.dot(Ri - Ti.dot(Pi)))
        self.VT[:self.nb_components, work] = newV[:self.nb_components, 0]

    def factorize(self, matrix, random_state):
        # Preprocessings
        matrixT = defaultdict(dict)
        for user in matrix:
            for work in matrix[user]:
                matrixT[work][user] = matrix[user][work]
        # Init
        self.U = np.random.rand(self.nb_users, self.nb_components +
                                self.nb_tags)
        self.VT = np.concatenate((np.random.rand(self.nb_components,
                                                 self.nb_works), self.T.T))
        # ALS
        for i in range(self.nb_iterations):
            self.M = self.U.dot(self.VT)
            # print('shape X_train', self.X_train.shape)
            y_pred = self.predict(self.X_train)
            print('Step {}: {}'.format(i, self.compute_rmse(
                self.predict(self.X_train), self.y_train)))
            for user in matrix:
                self.fit_user(user, matrix)
            for work in matrixT:
                self.fit_work(work, matrixT)

    def fit(self, X, y):
        self.X_train = X
        self.y_train = y
        if self.verbose_level:
            print("Computing M: (%i × %i)" % (self.nb_users, self.nb_works))
        matrix, self.means = self.make_matrix(X, y)

        self.chrono.save('fill and center matrix')

        self.factorize(matrix, random_state=42)
        if self.verbose_level:
            print('Shapes', self.U.shape, self.VT.shape)
        self.M = self.U.dot(self.VT)

        # self.save('backup.pickle')

        self.chrono.save('factor matrix')

    def predict(self, X):
        return (self.M[X[:, 0].astype(np.int64), X[:, 1].astype(np.int64)] +
                self.means[X[:, 0].astype(np.int64)])

    def get_shortname(self):
        return 'xals-%d' % self.nb_components

max_stars_repo_name: Akulen/mangaki-zero

max_issues_repo_name: Akulen/mangaki-zero

max_forks_repo_name: Akulen/mangaki-zero


text:
import os
import sys
import numpy
import pylab
import scipy.ndimage
import myutils as utils
import scipy.stats
import mahotas.polygon
import toolbox
import types
import math

LOOP_MAX_WIDTH = 400*1E-3 #400 microns
PIXELS_PER_MM_HOR = 320

"""
def rgb2gray(rgb):
    r, g, b = rgb[:,:,0], rgb[:,:,1], rgb[:,:,2]
    gray = 0.2989 * r + 0.5870 * g + 0.1140 * b
    return gray
"""
def img2float(fn):
    im = scipy.misc.imread(fn, flatten=True)
    return im #return rgb2gray(im)

def show_img(img,show=True, save=False, xy=None):
    #if show:
    #  pylab.imshow(img)
    #  pylab.gray()
    #  pylab.show()
    filename = os.tempnam()+".png"
    pylab.imsave(filename, img, cmap=pylab.cm.Greys_r)
    if xy is not None:
      x0 = xy[0]-5
      y0 = xy[1]-5
      x1 = x0+5
      y1 = y0+5
      os.system("convert %s -fill red -draw 'rectangle %d,%d,%d,%d' %s_" % (filename, x0, y0, x1, y1, filename))
      os.system("mv %s_ %s" % (filename, filename))
    if show:
      os.system("display %s" % filename)
    if not save:
      os.unlink(filename)
    else:
      print filename

def find_loop(img, debug=False, pixels_per_mm_horizontal=PIXELS_PER_MM_HOR, chi_angle=0):
    if type(img) == types.StringType:
        raw_img = img2float(img)
    else:
        # already a numpy array
        raw_img = img

    im = scipy.ndimage.gaussian_filter(raw_img, 1)

    if debug:
        show_img(im)

    im = im[1:-1,1:-2]
    im = utils.expand(im,2,mode="mirror")

    mask = new_mask(im.shape, .7, .7, 5)

    i1f = numpy.fft.fft2(im)
    i2f = numpy.multiply(i1f,mask)
    res = numpy.fft.ifft2(i2f)

    imgnumpy = im - res.real
    if debug:
        show_img(imgnumpy)

    dx = scipy.ndimage.sobel(imgnumpy,axis=0,mode='constant')  # horizontal derivative
    dy = scipy.ndimage.sobel(imgnumpy,axis=1,mode='constant')  # vertical derivative
    mag = numpy.hypot(dx,dy)  # magnitude
    mag *= 255.0/numpy.max(mag)
    if debug:
        show_img(mag)
    imgnumpy = numpy.uint8(mag)

    for offset in [20,15,10,5,0]:
        T = toolbox.get_robust_background_threshold(imgnumpy)
        thresholded = imgnumpy > (T+offset)
        if debug:
            show_img(thresholded)

        binimg = scipy.ndimage.morphology.binary_opening(thresholded, iterations=2)

        if debug:
            show_img(binimg)

        if numpy.sum(binimg) > 50:
          break
    else:
        return ("No loop detected", -1, -1)

    min1, max1, min2, max2 = mahotas.bbox(binimg)
    if debug:
        show_img(binimg)
    binimg2 = numpy.zeros_like(binimg)
    loop_max_width_pixels = LOOP_MAX_WIDTH*pixels_per_mm_horizontal
    if chi_angle == 0:
      # remove 200 microns from right
      binimg2[min1:max1, max2-loop_max_width_pixels:max2]=binimg[min1:max1, max2-loop_max_width_pixels:max2]
    else:
      # remove 200 micros from top
      binimg2[max1-loop_max_width_pixels:max1,  min2:max2]=binimg[max1-loop_max_width_pixels:max1,  min2:max2]
    if debug:
        show_img(binimg2)

    min1, max1, min2, max2 = mahotas.bbox(binimg2)
    bounded_img = binimg2[min1:max1, min2:max2]
    hull = mahotas.polygon.fill_convexhull(bounded_img)
    cy, cx = map(int, scipy.ndimage.measurements.center_of_mass(hull))
    x = int(min2)+cx
    y = int(min1)+cy
    if debug:
        show_img(binimg2, xy=(x,y))

    return ("Coord",x,y)

def new_mask(shape,sigma,sigma2,mulsigma):
    h = numpy.zeros(shape[0])
    w = numpy.zeros(shape[1])

    h[(shape[0]//2)-(mulsigma*sigma):(shape[0]//2)+(mulsigma*sigma)]=1
    w[(shape[1]//2)-(mulsigma*sigma2):(shape[1]//2)+(mulsigma*sigma2)]=1

    b1 = scipy.ndimage.filters.gaussian_filter(h,sigma)
    b2 = scipy.ndimage.filters.gaussian_filter(w,sigma2)

    h0 = numpy.empty_like(b1)
    h1 = numpy.empty_like(b2)

    h0[:shape[0] // 2] = b1[shape[0] - shape[0] // 2:]
    h0[shape[0] // 2:] = b1[:shape[0] - shape[0] // 2]
    h1[:shape[1] // 2] = b2[shape[1] - shape[1] // 2:]
    h1[shape[1] // 2:] = b2[:shape[1] - shape[1] // 2]

    g = numpy.outer(h0,h1)
    return g

if __name__ == '__main__':
    for filename in sys.argv[1:]:
        print find_loop(filename, debug=False) #True)

max_stars_repo_name: mxcube/lucid

max_issues_repo_name: mxcube/lucid

max_forks_repo_name: mxcube/lucid

reason: import numpy,import scipy


text:
"""Contains image metrics."""

import numpy as np

def ssim(x: np.ndarray, y: np.ndarray) -> float:
    """Calculate the structural similarity metric (SSIM) between images.

    Note:
        This function requires skimage.
    """
    from skimage.metrics import structural_similarity

    assert x.dtype == y.dtype

    return structural_similarity(x, y, channel_axis=2, multichannel=True)

def psnr(x: np.ndarray, y: np.ndarray) -> float:
    """Takes two images, returns the peak signal-to-noise ratio."""
    assert x.dtype == y.dtype

    mserr = mse(x, y)
    max_i = 255 if x.dtype == np.uint8 else 1

    return 10 * np.log10((max_i * max_i) / mserr)

def mse(x: np.ndarray, y: np.ndarray) -> float:
    """Takes two images, returns the mean-squared-error."""
    assert x.dtype == y.dtype

    return np.mean(np.square(x - y))

max_stars_repo_name: tilleyd/compy

max_issues_repo_name: tilleyd/compy

max_forks_repo_name: tilleyd/compy


text:
import argparse
from operator import imod
import os
import sys
import numpy as np
from numpy.lib.nanfunctions import _nancumprod_dispatcher

from warmup_scheduler import scheduler

PROJECT_ROOT = '/Users/dmitry/Documents/Personal/study/kaggle/tantum'
sys.path.append(PROJECT_ROOT)

import glob
import albumentations
import pandas as pd
import tantum
import torch
import torch.nn as nn
from efficientnet_pytorch import EfficientNet
from sklearn import metrics, model_selection, preprocessing
from torchvision.datasets import STL10

from torch.optim.lr_scheduler import CosineAnnealingWarmRestarts, CosineAnnealingLR, ReduceLROnPlateau

from tantum.callbacks import EarlyStopping, OutOfFold, WeightAndBiasesCallback, Checkpoint
from tantum.datasets.datasets import ImageDataset
from tantum.trainer.v1.base_fitter import Model
from torch.nn import functional as F
from tantum.model.simple_net import Net
from tantum.datasets.datasets import create_folds, MnistDataset

from tantum.datasets.augmentations import cutmix, mixup, fmix
from tantum.model.attention.multihead import MultiHeadAttention

from tantum.model import get_model
from tantum.optimizer import get_optimizer
from tantum.loss import get_criterion
from tantum.scheduler import get_scheduler

INPUT_PATH = "../../../input"
MODEL_PATH = "../models/"
MODEL_NAME = os.path.basename(__file__)[:-3]
TRAIN_BATCH_SIZE = 32
VALID_BATCH_SIZE = 32
IMAGE_SIZE = 192
EPOCHS = 20

class Config:
    # target_size=10
    model_type='Net'
    optim='AdamW'
    criterion='CrossEntropyLoss'
    scheduler='CosineAnnealingWarmRestarts'

    scheduler_params={'T_0':10, "last_epoch":-1, 'eta_min':1e-4, "T_mult":1}
    model_params = {'target_size': 10}
    use_cutmix=False

class ConvPart(nn.Module):
    def __init__(self):
        super(ConvPart, self).__init__()
        self.conv1 = nn.Conv2d(1, 10, kernel_size=5)
        self.conv2 = nn.Conv2d(10, 20, kernel_size=5)
        self.conv2_drop = nn.Dropout2d()

    def forward(self, x):
        x = x.view(-1,1,28,28)
        x = F.relu(F.max_pool2d(self.conv1(x), 2))
        x = F.relu(F.max_pool2d(self.conv2_drop(self.conv2(x)), 2))
        return x

class NetMultiheadAttention(nn.Module):
    def __init__(self):
        super().__init__()
        self.conv = ConvPart()
        self.attn1 = MultiHeadAttention(n_head=4, d_model=20, d_k=8, d_v=8)
        self.optim = torch.optim.Adam(self.parameters(), lr=1e-4)
        self.fc1 = nn.Linear(320, 50)
        self.fc2 = nn.Linear(50, 10)
        self.final = nn.Linear(20,10)

    def forward(self, x):
        z = self.conv(x)
        q = torch.reshape(z, (z.size(0), -1 , z.size(1)))
        q, w = self.attn1(q, q, q)
        q = torch.reshape(q, (z.size(0), z.size(1), z.size(2), z.size(3)))
        z = q.mean(3).mean(2)
        z = q.view(-1, 320)
        z = self.fc1(z)
        z = self.fc2(z)
        return F.log_softmax(z, -1), q

class MnistModel(Model):
    def __init__(self):
        super().__init__()
        # self.model = get_model(Config)()
        self.model = NetMultiheadAttention()
        self.opt = get_optimizer(Config)(self.parameters(), lr=1e-4)
        self.loss_fn = get_criterion(Config)()

    def monitor_metrics(self, outputs, targets):
        outputs = torch.argmax(outputs, dim=1).cpu().detach().numpy()
        targets = targets.cpu().detach().numpy()
        accuracy = metrics.accuracy_score(targets, outputs)
        return {"accuracy": accuracy}

    def fetch_scheduler(self, train_loader):
        scheduler = get_scheduler(Config)(self.opt, **Config.scheduler_params)
        return scheduler

    def fetch_optimizer(self):
        return self.opt

    def forward(self, image, targets=None):
        mix_decision = np.random.rand()
        image = image.float()
        targets = targets.long()
        if Config.use_cutmix and mix_decision < 0.25:
            image, image_labels = cutmix(image, targets, 1.)
        # elif mix_decision >=0.25 and mix_decision < 0.5:
            # image, image_labels = fmix(image, targets, alpha=1., decay_power=5., shape=(28,28))

        outputs, att = self.model(image.float())

        if targets is not None:

            if Config.use_cutmix and mix_decision < 0.25:
                loss = self.loss_fn(outputs, image_labels[0]) * image_labels[2] + self.loss_fn(outputs, image_labels[1]) * (1. - image_labels[2])
                metrics = self.monitor_metrics(outputs, image_labels[0])
                return outputs, loss, metrics
            else:
                loss = self.loss_fn(outputs, targets)
                metrics = self.monitor_metrics(outputs, targets)

            return outputs, loss, metrics
        return outputs, 0, {}
# dataset = STL10("stl10", split='train', download=True)
if __name__ == "__main__":

    ## Reading Data
    df = pd.read_csv(os.path.join(PROJECT_ROOT, 'train.csv') )
    ## Create Folds
    folds = create_folds(df, 5)

    ## Callbacks

    oof_callback = OutOfFold(output_path='./')

    es_callback = EarlyStopping(
        monitor="valid_loss",
        model_path=os.path.join('./', 'mnist' + ".pt"),
        patience=3,
        mode="min",
    )

    chc_callback = Checkpoint(save_path='./', file_name='mnist_checkpoint')

    # wb_callback = WeightAndBiasesCallback(project='mnist')

    for fold in range(5):

        model = MnistModel()

        trn_idx = folds[folds['fold'] != fold].index
        val_idx = folds[folds['fold'] == fold].index

        train_folds = folds.loc[trn_idx].reset_index(drop=True)
        valid_folds = folds.loc[val_idx].reset_index(drop=True)

        y = train_folds['label'].values
        X = train_folds.drop(['label','fold'], axis=1).values

        y_val = valid_folds['label'].values
        X_val = valid_folds.drop(['label','fold'], axis=1).values

        valid_labels = valid_folds['label'].values

        train_dataset = MnistDataset(
            X = X,
            y = y,
        )

        valid_dataset = MnistDataset(
            X = X_val,
            y = y_val,
        )

        model.fit(
            train_dataset,
            valid_dataset=valid_dataset,
            train_bs=TRAIN_BATCH_SIZE,
            valid_bs=VALID_BATCH_SIZE,
            device="cpu",
            epochs=EPOCHS,
            callbacks=[],
            fp16=False,
            n_jobs=0,
            valid_labels=valid_labels,
            valid_folds=valid_folds,
            target_size=10,
            fold=fold
        )

max_stars_repo_name: dmitryshendryk/tantum

max_issues_repo_name: dmitryshendryk/tantum

max_forks_repo_name: dmitryshendryk/tantum

reason: import numpy,from numpy


text:
import scipy as _sp

def cylinder(target, throat_diameter='throat.diameter'):
    r"""
    Calcuate the throat perimeter assuming a circular cross-section

    Parameters
    ----------
    target : OpenPNM Object
        The object which this model is associated with. This controls the
        length of the calculated array, and also provides access to other
        necessary properties.

    throat_diameter : string
        The dictionary key of the array containing the throat diameter values
    """
    return target[throat_diameter]*_sp.pi

def cuboid(target, throat_diameter='throat.diameter'):
    r"""
    Calcuate the throat perimeter assuming a square cross-section

    Parameters
    ----------
    target : OpenPNM Object
        The object which this model is associated with. This controls the
        length of the calculated array, and also provides access to other
        necessary properties.

    throat_diameter : string
        The dictionary key of the array containing the throat diameter values
    """
    return target[throat_diameter]*4

max_stars_repo_name: bryanwweber/OpenPNM

max_issues_repo_name: kvmkrao/OpenPNM

max_forks_repo_name: kvmkrao/OpenPNM

reason: import scipy


text:
# This Python 3 environment comes with many helpful analytics libraries installed

# First Let us import all the required modules.

# For linear algebra
import numpy as np

# For data processing, CSV file I/O (e.g. pd.read_csv)
import pandas as pd

# For plotting the Curves
import matplotlib.pyplot as pt

# For Statistcal Graphics
import seaborn as sns

# For Optimiztion and Integration
from sklearn import preprocessing
from scipy.integrate import odeint
import time
from datetime import datetime
from scipy import integrate, optimize
import warnings
warnings.filterwarnings('ignore')

#Now,Let us import some ML libraries
import lightgbm as lgb
import xgboost as xgb
from xgboost import plot_importance, plot_tree
from sklearn.model_selection import RandomizedSearchCV, GridSearchCV
from sklearn import linear_model
from sklearn.metrics import mean_squared_error

# Taking the input file
train = pd.read_csv("../input/traincsv/train.csv")

# Creating the dataframes from the input file.
population = float(10617423)
country_df = pd.DataFrame()
country_df['ConfirmedCases'] = train.loc[train['Country_Region']=='Georgia'].ConfirmedCases.diff().fillna(0)
country_df = country_df[41:]
country_df['day_count'] = list(range(1,len(country_df)+1))

# Storing the data in required entities
ydata = [i for i in country_df.ConfirmedCases]
xdata = country_df.day_count
ydata = np.array(ydata, dtype=float)
xdata = np.array(xdata, dtype=float)

# Initialising the Data
N = population
inf0 = ydata[0]
sus0 = N - inf0
rec0 = 0.0

# Defining the basic ODE's of SIR Model
def sir_model(y, x, beta, gamma):
    sus = -beta * y[0] * y[1] / N
    rec = gamma * y[1]
    inf = -(sus + rec)
    return sus, inf, rec

# Defining the Function to integrate ODE's
def fit_odeint(x, beta, gamma):
    return integrate.odeint(sir_model, (sus0, inf0, rec0), x, args=(beta, gamma))[:,1]

# Storing the values of the OPtimized parameters.
Opt, pcov = optimize.curve_fit(fit_odeint, xdata, ydata)
fitted = fit_odeint(xdata, *Opt)

# Plotting the fitting of curve into data points.
print("Fitting a SECOND order equation into Infectives data yields:-")
pt.plot(xdata, ydata, 'o')
pt.plot(xdata, fitted)
pt.title("Fit of SIR model for Georgia infected cases")
pt.ylabel("Population infected")
pt.xlabel("Days")
pt.show()

# Total population
N = population
# Initialization of infected and recovered individuals, I0 and R0.
I0, R0 = 2, 0
# S0 be the susceptible to infection initially.
S0 = N - I0 - R0
# Contact rate(beta) and mean recovery rate(gamma).
beta, gamma = Opt[0],Opt[1]
# A grid of time points (in days)
t = np.linspace(0, 200, 200)

# The SIR model definition in differential equations.
def der(y, t, N, beta, gamma):
    S, I, R = y
    dSdt = -beta * S * I / N
    dIdt = beta * S * I / N - gamma * I
    dRdt = gamma * I
    return dSdt, dIdt, dRdt

# Initialization
y0 = S0, I0, R0
# Integrating the SIR equations over the time grid, t.
ret = odeint(der, y0, t, args=(N, beta, gamma))
S, I, R = ret.T

# Plotting the data on three separate curves for S(t), I(t) and R(t)
fig = pt.figure(facecolor='y')
ax = fig.add_subplot(111, facecolor='#dddddd', axisbelow=True)
ax.plot(t, S/(1.25*N), 'b', label='Susceptible')
ax.plot(t, I/5000, 'r', label='Infected')
ax.plot(t, R/400000, 'g', label='Recovered ')
ax.set_xlabel('Time /days')
ax.set_ylabel('Number (1000s)')
ax.set_ylim(0,1.2)
ax.yaxis.set_tick_params(length=0)
ax.xaxis.set_tick_params(length=0)
ax.grid(b=True, which='major', c='w', lw=2, ls='-')
legend = ax.legend()
legend.get_frame().set_alpha(0.5)
for spine in ('top', 'right', 'bottom', 'left'):
    ax.spines[spine].set_visible(False)
pt.show()

print("1.0. The Model Parameters are:-")
r=Opt[0]/N
a=Opt[1]
print("1.1. Transmissibility(r)          =",Opt[0]/N)
print("1.2. Recovery Rate(a)             =",Opt[1])
print("1.3. Relative Removal Rate (Rho)  =",a/r)
print("1.4. Contact rate(Sigma)          =",r/a)
print("2.0. Basic Reproduction Rate(Ro)  =",(r*(10617421))/a)
print("3.0. Maximum Number of Infectives =",max(I))

#Since Max No of infectives occur ar S~=Rho
temp=[]
for i in range(len(S)):
        if (S[i]-(a/r))>0:
            temp.append([(S[i]-(a/r)),i])
#Since min(temp) is at 85th day,It has max no of Infectives
print("4.0. The time location of maximum infectives is at 85th day from the day of initial infection which means 26-05-2020.")

temp2=[]
for i in range(len(I)):
    if I[i]<1:
        x=i
        break

#Since at 179th day from the initial day of infection, Infectives become less than 1, It is duration of epidemic.
print("5.0. The total number of Recovered people is around",R[i])
print("6.0. The Duration of the epidemic is 179 days.So it is expected that epidemic ends on 18-8-2020 in Georgia,US.")

print("Hence,The SIR Model is fitted into the COVID 19 data of US state,Georgia. All the model parameters and other required quantities are estimated")

max_stars_repo_path: Machine Learning Projects/Covid-19 Projects/SIR Model/SIR_Model.py


text:
import sys
import os
import json
import networkx as nx
from ast import literal_eval
sys.path.append(os.path.join(os.path.dirname(__file__), '../..'))
import multi_layer_network.src.event_baseline2 as eb2

def run(entity_edgelist, entity_heads, event_heads, outputs_prefix):
    # entity_heads will be modified in eb2.get_resolved_entity: ('entity2.json')
    entity2cluster = eb2.get_resolved_entity(entity_edgelist, entity_heads)
    print(len(entity2cluster))
    path_to_output = outputs_prefix + "event.jl"
    event_jl = eb2.event_baseline_linking(event_heads, entity2cluster)
    with open(path_to_output, 'w') as f:
        for line in event_jl:
            json.dump(line, f)
            f.write('\n')
    return event_jl

def load_edgelist_from_file(edgelist_file):
    G = nx.Graph()
    with open(edgelist_file, "r") as edges:
        G.add_nodes_from(literal_eval(edges.readline()))
        for edge in edges:
            edge_nodes = literal_eval(edge)
            G.add_edge(edge_nodes[0], edge_nodes[1])
    return G

def run_with_file_io(edgelist_file, entity_head_file, event_head_file, outputs_prefix):
    '''

    :param edgelist_file: file path of the entity edgelist
    :param entity_head_file: file path of the entity json head
    :param event_head_file: file path of the event json head
    :param outputs_prefix: path of the output folder of jl file
    :return: event jl
    '''
    edgelist = load_edgelist_from_file(edgelist_file)
    with open(entity_head_file) as f:
        entity_head = json.load(f)
    with open(event_head_file) as f:
        event_head = json.load(f)
    event_jl = run(edgelist, entity_head, event_head, outputs_prefix)
    return event_jl

max_stars_repo_name: shinyichen/gaia-clustering

max_issues_repo_name: shinyichen/gaia-clustering

max_forks_repo_name: shinyichen/gaia-clustering

reason: import networkx


text:
# Author: Matthew Wicker
# Companion Code for paper: Analysis of 3D Deep Learning in an Adversarial Setting
# CVPR 2019

"""
This file impliments the PointNet model (Qi et. al. 2017) in keras
It is aware of weights that are saved in the Models directory of this
repository. So if you would like to modify/retrain this model, then
please ensure the weights and architecture are changed accordingly.
"""

import h5py
import numpy as np
import numpy as np
import os
import tensorflow as tf
from keras import optimizers
from keras.layers import Input
from keras.models import Model
from keras.layers import Dense, Flatten, Reshape, Dropout
from keras.layers import Convolution1D, MaxPooling1D, BatchNormalization
from keras.layers import Lambda
from keras.utils import np_utils
from keras import backend as K
import copy
def mat_mul(A, B):
    return tf.matmul(A, B)

"""
This function declares the PointNet architecture.
@ Param classes - Integer, defining the number of classes that the
		 model will be predicting.
@ Param load_weights - Boolean, if classes is set for 10 or 40 we
		will load pretrained weights into the model.
"""
def PointNet(classes=40, load_weights=True, num_points =2048):
	num_points = num_points
	input_points = Input(shape=(num_points, 3))
	x = Convolution1D(64, 1, activation='relu',
		          input_shape=(num_points, 3))(input_points)
	x = BatchNormalization()(x)
	x = Convolution1D(128, 1, activation='relu')(x)
	x = BatchNormalization()(x)
	x = Convolution1D(1024, 1, activation='relu')(x)
	x = BatchNormalization()(x)
	x = MaxPooling1D(pool_size=num_points)(x)
	x = Dense(512, activation='relu')(x)
	x = BatchNormalization()(x)
	x = Dense(256, activation='relu')(x)
	x = BatchNormalization()(x)
	x = Dense(9, weights=[np.zeros([256, 9]), np.array([1, 0, 0, 0, 1, 0, 0, 0, 1]).astype(np.float32)])(x)
	input_T = Reshape((3, 3))(x)

	# forward net
	g = Lambda(mat_mul, arguments={'B': input_T})(input_points)
	g = Convolution1D(64, 1, input_shape=(num_points, 3), activation='relu')(g)
	g = BatchNormalization()(g)
	g = Convolution1D(64, 1, input_shape=(num_points, 3), activation='relu')(g)
	g = BatchNormalization()(g)

	# feature transform net
	f = Convolution1D(64, 1, activation='relu')(g)
	f = BatchNormalization()(f)
	f = Convolution1D(128, 1, activation='relu')(f)
	f = BatchNormalization()(f)
	f = Convolution1D(1024, 1, activation='relu')(f)
	f = BatchNormalization()(f)
	f = MaxPooling1D(pool_size=num_points)(f)
	f = Dense(512, activation='relu')(f)
	f = BatchNormalization()(f)
	f = Dense(256, activation='relu')(f)
	f = BatchNormalization()(f)
	f = Dense(64 * 64, weights=[np.zeros([256, 64 * 64]), np.eye(64).flatten().astype(np.float32)])(f)
	feature_T = Reshape((64, 64))(f)

	# forward net
	g = Lambda(mat_mul, arguments={'B': feature_T})(g)
	g = Convolution1D(64, 1, activation='relu')(g)
	g = BatchNormalization()(g)
	g = Convolution1D(128, 1, activation='relu')(g)
	g = BatchNormalization()(g)
	g = Convolution1D(1024, 1, activation='relu')(g)
	g = BatchNormalization()(g)

	# global_feature
	global_feature = MaxPooling1D(pool_size=num_points)(g)

	# point_net_cls
	c = Dense(512, activation='relu')(global_feature)
	c = BatchNormalization()(c)
	c = Dropout(rate=0.7)(c)
	c = Dense(256, activation='relu')(c)
	c = BatchNormalization()(c)
	c = Dropout(rate=0.7)(c)
	c = Dense(classes, activation='softmax')(c)
	prediction = Flatten()(c)

	model = Model(inputs=input_points, outputs=prediction)

	model.compile(optimizer='adam',
              loss='categorical_crossentropy',
              metrics=['accuracy'])

	if(classes == 40 and load_weights):
    		model.load_weights('Models/PointNet-ModelNet40.h5')
	if(classes == 10 and load_weights):
    		model.load_weights('Models/PointNet-ModelNet10.h5')
	if(classes == 2 and load_weights):
    		model.load_weights('Models/PointNet-KITTI.h5')
	if(classes == 3 and load_weights):
    		model.load_weights('Models/PointNet-KITTI3.h5')
	print model.summary()

	return model

def predict(x_in, model):
    x_in = np.squeeze(x_in)
    val = model.predict(np.asarray([x_in]))
    val = np.squeeze(val)
    cl = np.argmax(val)
    return val[cl], cl

"""
This method returns the activations of the max pooling layer for
the specified inputs. Expects the following input
@Param - model, the Keras model that is outputted from the PointNet() function
@Param - point_cloud, the point cloud input that we want the maxpooling layer for
"""
def get_max_pool(model, point_cloud):
	layer_name = 'max_pooling1d_3'
	intermediate_layer_model = Model(inputs=model.input,
		                     outputs=model.get_layer(layer_name).output)
	inp = np.asarray(point_cloud)
	activations = intermediate_layer_model.predict(inp)
	return activations

def get_latent_activations(model, point_cloud):
	layer_name = 'batch_normalization_15'
	intermediate_layer_model = Model(inputs=model.input,
                                     outputs=model.get_layer(layer_name).output)
	value_test = np.asarray([point_cloud])
	intermediate_output = intermediate_layer_model.predict(value_test)
	return intermediate_output

def get_critical_set(model, point_cloud):
	latent = get_latent_activations(model, point_cloud)[0]
	critical_set = np.argmax(latent, axis=0)
	critical_set = set(critical_set)
	return critical_set

def get_critical_set_bb(model, point_cloud):
    critical_set = []
    values = []
    #v_init, c_init = predict(point_cloud, model)
    pc = copy.deepcopy(point_cloud)
    for i in range(len(pc)):
        val = copy.deepcopy(pc[i])
        pc[i] = [0.0,0.0,0.0]
        v, c = predict(pc, model)
        values.append(v)
        pc[i] = val
    unique = np.unique(values,return_index=True)
    return unique[1]


text:
# Remember for seqlogo you need to run it using python3
import argparse
import numpy as np
import os
import pandas as pd

parser = argparse.ArgumentParser()
parser.add_argument("-r","--folderToRead",type=str,help="Input a folder in which you want to read files from")
parser.add_argument("-w","--folderToWrite",type=str,help="Input a folder in which you want to write files into")
parser.add_argument("-i","--IDfile",type=str,help="Input an ID file to identify the motifs and match them to the RBP")

args = parser.parse_args()

TMPfolder=args.folderToRead
WriteFolder=args.folderToWrite
idFile=args.IDfile

# Format a float number to have 3 decimal places
def format(value):
    return "%.3f" % value

#Get IDs that are relevant for Homo sapiens
with open(idFile) as f:
    idlines=[line.strip().split("\t") for line in f][1:]

relevantIDs = [i[0] for i in idlines if i[3]=="Homo_sapiens"]

# Go through every PFM file in the folder
for filename in os.listdir(TMPfolder):
    if (".txt" in filename) and (filename.split("_")[0] in relevantIDs):
        data=pd.read_csv(TMPfolder+filename,sep="\t",header=0,index_col=0)
        data_matrix = data.values

        PWM_cluster = np.log2(data_matrix/0.25)
        writefilename="_".join(filename.split(".txt")[0].split("_")[0:2])+"_pwm.txt"
        #print(PWM_cluster)
        with open(WriteFolder+"PWMs_RBPs_Compendium/"+writefilename,"w") as fw:
            for line in PWM_cluster:
                for i in line:
                    fw.write(str(format(i)))
                    fw.write("\t")
                fw.write("\n")


text:
## dual axis for the x-axes (relation between photon energy and wavelength)
def top_tick_function(p): return 1e6*(np.pi*2)/p
ax3 = ax.twiny()
ax3.set_xscale('log') ## LOG!
#ax2.axis['top'].major_ticklabels.set_visible(False)
ax3.set_xlim(np.array(ax.get_xlim()))
ax3.set_xlabel('characteristic feature size (μm)')

## If we wish nice round numbers on the secondary axis ticks...
## ... we must give them correct positions. To do so, we need to numerically invert the tick values:
top_ax_limits = ([top_tick_function(lim) for lim in ax.get_xlim()])
xticks2 = 10**(matplotlib.ticker.MaxNLocator(nbins=8, steps=[1]).tick_values(np.log10(top_ax_limits[1]), np.log10(top_ax_limits[0])))
from scipy.optimize import brentq
def top_tick_function_inv(r, target_value=0):
    return (brentq(lambda r,q: top_tick_function(r) - target_value, ax.get_xlim()[1], ax.get_xlim()[0], 1e-6))
valid_xtick2loc, valid_xtick2val = [], []
for xtick2 in xticks2:
    try:
        valid_xtick2loc.append(top_tick_function_inv(0, xtick2))
        valid_xtick2val.append(xtick2)
    except ValueError:      ## (skip tick if the ticker.MaxNLocator gave invalid target_value to brentq optimization)
        pass
## Finally, we set the positions and tick values
print("valid_xtick2loc",valid_xtick2loc,"valid_xtick2val", valid_xtick2val)
ax3.set_xticks(valid_xtick2loc)
from matplotlib.ticker import FixedFormatter, NullFormatter
ax3.xaxis.set_major_formatter(FixedFormatter(["%g" % toptick for toptick in valid_xtick2val]))
ax3.set_minor_formatter(NullFormatter())

max_stars_repo_name: FilipDominec/gayplot

max_issues_repo_name: FilipDominec/gayplot

max_forks_repo_name: FilipDominec/gayplot

reason: from scipy


text:
# -*- coding: utf-8 -*-
'''
@time: 2019/9/12 15:16

@ author: javis
'''
import torch
import numpy as np
import time, os
from sklearn.metrics import f1_score, roc_auc_score
from torch import nn

def mkdirs(path):
    if not os.path.exists(path):
        os.makedirs(path)

# 计算F1score,每一类单独计算最后计算平均的auc
def calc_auc(y_true, y_pre, threshold=0.5):
    labels = y_true.cpu().detach().numpy().astype(np.int)
    outputs = y_pre.cpu().detach().numpy()
    return roc_auc_score(labels, outputs, 'micro')

# 打印时间
def print_time_cost(since):
    time_elapsed = time.time() - since
    return '{:.0f}m{:.0f}s\n'.format(time_elapsed // 60, time_elapsed % 60)

# 调整学习率
def adjust_learning_rate(optimizer, lr):
    for param_group in optimizer.param_groups:
        param_group['lr'] = lr
    return lr

# 多标签使用类别权重
class WeightedMultilabel(nn.Module):
    def __init__(self, weights: torch.Tensor):
        super(WeightedMultilabel, self).__init__()
        self.cerition = nn.BCEWithLogitsLoss(reduction='none')
        self.weights = weights

    def forward(self, outputs, targets):
        loss = self.cerition(outputs, targets)
        return (loss * self.weights).mean()


text:
import networkx as nx
import utils

def comp(G):
    for node in G.node:
        G.node[node]['seq']=utils.rc(G.node[node]['seq'])

    genome2length=dict()
    #relabel the offsets, determine the length of all genomes in the graph, then l-pos
    for sample in G.graph['paths']:
        maxp=0
        for node,data in G.nodes(data=True):
            if sample in data['offsets']:
                if data['offsets'][sample]+len(data['seq'])>maxp:
                    maxp=data['offsets'][sample]+len(data['seq'])
        genome2length[sample]=maxp

    for sample in G.graph['paths']:
        for node,data in G.nodes(data=True):
            if sample in data['offsets']:
                G.node[node]['offsets'][sample]=genome2length[sample]-(G.node[node]['offsets'][sample]+len(data['seq']))

    G.reverse(copy=False)
    return G

def comp_cmd(args):
    g=nx.DiGraph()
    g.graph['paths']=[]
    utils.read_gfa(args.graph[0],None,None,g,targetsample=None)
    g=comp(g)
    utils.write_gfa(g,"",outputfile=args.graph[0].replace('.gfa','.rc.gfa'), nometa=False)

max_stars_repo_name: ChriKub/reveal

max_issues_repo_name: ChriKub/reveal

max_forks_repo_name: ChriKub/reveal


text:
# -*- coding: utf-8 -*-
#
from __future__ import division

import math
import numpy
import scipy.special
import sympy

def untangle(data):
    weights, points = zip(*data)
    return (
        numpy.concatenate(points),
        numpy.repeat(weights, [len(grp) for grp in points]),
    )

def n_outer(a):
    """Given a list (tuple, array) of arrays, this method computes their outer
    product. If the dimension of the input arrays is larger than one, the
    product is formed across the first dimension; all other dimensions must
    coincide in size.

    Examples:
    n_outer([np.ones(4), np.ones(5)]).shape == (4, 5)
    n_outer([np.ones(4), np.ones(5), np.ones(6)]).shape == (4, 5, 6)
    n_outer([np.ones(4, 3, 7), np.ones(5, 3, 7)]).shape == (4, 5, 3, 7)
    """
    # <https://stackoverflow.com/a/45376730/353337>
    d = len(a)

    # If the elements are more than one-dimensional, assert that the extra
    # dimensions are all equal.
    s0 = a[0].shape
    for arr in a:
        assert s0[1:] == arr.shape[1:]

    out = a[0]
    for k in range(1, d):
        # Basically outer products. Checkout `numpy.outer`'s implementation for
        # comparison.
        out = numpy.multiply(
            # Insert a newaxis after k `:`
            out[(slice(None),) * k + (numpy.newaxis,)],
            # Insert a newaxis at the beginning
            a[k][numpy.newaxis],
        )
    return out

def compute_dobrodeev(n, I0, I2, I22, I4, pm_type, i, j, k, symbolic=False):
    """Compute some helper quantities used in

    L.N. Dobrodeev,
    Cubature rules with equal coefficients for integrating functions with
    respect to symmetric domains,
    USSR Computational Mathematics and Mathematical Physics,
    Volume 18, Issue 4, 1978, Pages 27-34,
    <https://doi.org/10.1016/0041-5553(78)90064-2>.
    """
    t = 1 if pm_type == "I" else -1

    binomial = sympy.binomial if symbolic else scipy.special.binom
    fact = sympy.factorial if symbolic else math.factorial
    sqrt = sympy.sqrt if symbolic else numpy.sqrt

    L = binomial(n, i) * 2 ** i
    M = fact(n) // (fact(j) * fact(k) * fact(n - j - k)) * 2 ** (j + k)
    N = L + M
    F = I22 / I0 - I2 ** 2 / I0 ** 2 + (I4 / I0 - I22 / I0) / n
    R = (
        -(j + k - i) / i * I2 ** 2 / I0 ** 2
        + (j + k - 1) / n * I4 / I0
        - (n - 1) / n * I22 / I0
    )
    H = (
        1
        / i
        * (
            (j + k - i) * I2 ** 2 / I0 ** 2
            + (j + k) / n * ((i - 1) * I4 / I0 - (n - 1) * I22 / I0)
        )
    )
    Q = L / M * R + H - t * 2 * I2 / I0 * (j + k - i) / i * sqrt(L / M * F)

    G = 1 / N
    a = sqrt(n / i * (I2 / I0 + t * sqrt(M / L * F)))
    b = sqrt(n / (j + k) * (I2 / I0 - t * sqrt(L / M * F) + t * sqrt(k / j * Q)))
    c = sqrt(n / (j + k) * (I2 / I0 - t * sqrt(L / M * F) - t * sqrt(j / k * Q)))
    return G, a, b, c

max_stars_repo_name: gdmcbain/quadpy

max_issues_repo_name: gdmcbain/quadpy

max_forks_repo_name: gdmcbain/quadpy

reason: import numpy,import scipy,import sympy