text:
#================================PlotFuncs.py==================================#
# Created by Ciaran O'Hare 2020

# Description:
# This file has many functions which are used throughout the project, but are
# all focused around the bullshit that goes into making the plots

#==============================================================================#

from numpy import *
from numpy.random import *
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import matplotlib.gridspec as gridspec
from matplotlib.colors import ListedColormap
from matplotlib import colors
from scipy.spatial import Delaunay
from scipy.spatial import ConvexHull
from scipy.stats import zscore,chi2,multivariate_normal
from scipy.special import erfinv
from scipy.stats import gaussian_kde
import matplotlib.ticker as mticker
from mpl_toolkits.axes_grid1 import make_axes_locatable
from scipy.stats import norm
from Params import *

#==============================================================================#
# My custom figure save
def MySaveFig(fig,pltname,pngsave=True):
    fig.savefig(pltdir+pltname+'.pdf',bbox_inches='tight')
    if pngsave:
        fig.savefig(pltdir_png+pltname+'.png',bbox_inches='tight')
#==============================================================================#

#==============================================================================#
# My preferred style of colorbar
def cbar(mappable,extend='neither',\
                label='',lfs=35,labelpad=40,labelrotation=-90):
    ax = mappable.axes
    fig = ax.figure
    divider = make_axes_locatable(ax)
    cax = divider.append_axes("right", size="5%", pad=0.2)
    cbar = fig.colorbar(mappable, cax=cax,extend=extend)
    cbar.ax.tick_params(which='minor',length=8,width=2)
    cbar.ax.tick_params(which='major',length=10,width=3)
    cbar.set_label(label,fontsize=lfs,rotation=labelrotation,labelpad=labelpad)
    return cbar
#==============================================================================#

#==============================================================================#
def MySquarePlot(xlab='',ylab='',\
                 lw=2.5,lfs=45,tfs=25,size_x=13,size_y=12,Grid=False):
    plt.rcParams['axes.linewidth'] = lw
    plt.rc('text', usetex=True)
    plt.rc('font', family='serif',size=tfs)
    mpl.rcParams['text.latex.preamble'] = [r'\usepackage{mathpazo}']

    fig = plt.figure(figsize=(size_x,size_y))
    ax = fig.add_subplot(111)

    ax.set_xlabel(xlab,fontsize=lfs)
    ax.set_ylabel(ylab,fontsize=lfs)

    ax.tick_params(which='major',direction='in',width=2,length=13,right=True,top=True,pad=7)
    ax.tick_params(which='minor',direction='in',width=1,length=10,right=True,top=True)
    if Grid:
        ax.grid()
    return fig,ax

def MyDoublePlot(xlab1='',ylab1='',xlab2='',ylab2='',\
                 wspace=0.25,lw=2.5,lfs=45,tfs=25,size_x=20,size_y=11,Grid=False):
    plt.rcParams['axes.linewidth'] = lw
    plt.rc('text', usetex=True)
    plt.rc('font', family='serif',size=tfs)
    mpl.rcParams['text.latex.preamble'] = [r'\usepackage{mathpazo}']

    fig, axarr = plt.subplots(1, 2,figsize=(size_x,size_y))
    gs = gridspec.GridSpec(1, 2)
    gs.update(wspace=wspace)
    ax1 = plt.subplot(gs[0])
    ax2 = plt.subplot(gs[1])
    ax1.tick_params(which='major',direction='in',width=2,length=13,right=True,top=True,pad=7)
    ax1.tick_params(which='minor',direction='in',width=1,length=10,right=True,top=True)
    ax2.tick_params(which='major',direction='in',width=2,length=13,right=True,top=True,pad=7)
    ax2.tick_params(which='minor',direction='in',width=1,length=10,right=True,top=True)

    ax1.set_xlabel(xlab1,fontsize=lfs)
    ax1.set_ylabel(ylab1,fontsize=lfs)

    ax2.set_xlabel(xlab2,fontsize=lfs)
    ax2.set_ylabel(ylab2,fontsize=lfs)

    if Grid:
        ax1.grid()
        ax2.grid()
    return fig,ax1,ax2

def MyDoublePlot_Vertical(xlab1='',ylab1='',xlab2='',ylab2='',\
                 hspace=0.25,lw=2.5,lfs=45,tfs=35,size_x=20,size_y=10,Grid=False):
    plt.rcParams['axes.linewidth'] = lw
    plt.rc('text', usetex=True)
    plt.rc('font', family='serif',size=tfs)
    mpl.rcParams['text.latex.preamble'] = [r'\usepackage{mathpazo}']

    fig, axarr = plt.subplots(2,1,figsize=(size_x,size_y))
    gs = gridspec.GridSpec(2, 1)
    gs.update(hspace=hspace)
    ax1 = plt.subplot(gs[0])
    ax2 = plt.subplot(gs[1])

    ax1.tick_params(which='major',direction='in',width=2,length=13,right=True,top=True,pad=7)
    ax1.tick_params(which='minor',direction='in',width=1,length=10,right=True,top=True)

    ax2.tick_params(which='major',direction='in',width=2,length=13,right=True,top=True,pad=7)
    ax2.tick_params(which='minor',direction='in',width=1,length=10,right=True,top=True)

    ax1.set_xlabel(xlab1,fontsize=lfs)
    ax1.set_ylabel(ylab1,fontsize=lfs)

    ax2.set_xlabel(xlab2,fontsize=lfs)
    ax2.set_ylabel(ylab2,fontsize=lfs)

    if Grid:
        ax1.grid()
        ax2.grid()
    return fig,ax1,ax2

def MyTriplePlot(xlab1='',ylab1='',xlab2='',ylab2='',xlab3='',ylab3='',\
                 wspace=0.25,lw=2.5,lfs=45,tfs=25,size_x=20,size_y=7,Grid=False):
    plt.rcParams['axes.linewidth'] = lw
    plt.rc('text', usetex=True)
    plt.rc('font', family='serif',size=tfs)
    mpl.rcParams['text.latex.preamble'] = [r'\usepackage{mathpazo}']

    fig, axarr = plt.subplots(1, 3,figsize=(size_x,size_y))
    gs = gridspec.GridSpec(1, 3)
    gs.update(wspace=wspace)
    ax1 = plt.subplot(gs[0])
    ax2 = plt.subplot(gs[1])
    ax3 = plt.subplot(gs[2])

    ax1.tick_params(which='major',direction='in',width=2,length=13,right=True,top=True,pad=7)
    ax1.tick_params(which='minor',direction='in',width=1,length=10,right=True,top=True)

    ax2.tick_params(which='major',direction='in',width=2,length=13,right=True,top=True,pad=7)
    ax2.tick_params(which='minor',direction='in',width=1,length=10,right=True,top=True)

    ax3.tick_params(which='major',direction='in',width=2,length=13,right=True,top=True,pad=7)
    ax3.tick_params(which='minor',direction='in',width=1,length=10,right=True,top=True)

    ax1.set_xlabel(xlab1,fontsize=lfs)
    ax1.set_ylabel(ylab1,fontsize=lfs)

    ax2.set_xlabel(xlab2,fontsize=lfs)
    ax2.set_ylabel(ylab2,fontsize=lfs)

    ax3.set_xlabel(xlab3,fontsize=lfs)
    ax3.set_ylabel(ylab3,fontsize=lfs)

    if Grid:
        ax1.grid()
        ax2.grid()
        ax3.grid()
    return fig,ax1,ax2,ax3

def MyTriplePlot_Vertical(xlab1='',ylab1='',xlab2='',ylab2='',xlab3='',ylab3='',\
                 hspace=0.25,lw=2.5,lfs=45,tfs=35,size_x=20,size_y=15,Grid=False):
    plt.rcParams['axes.linewidth'] = lw
    plt.rc('text', usetex=True)
    plt.rc('font', family='serif',size=tfs)
    mpl.rcParams['text.latex.preamble'] = [r'\usepackage{mathpazo}']

    fig, axarr = plt.subplots(3,1,figsize=(size_x,size_y))
    gs = gridspec.GridSpec(3, 1)
    gs.update(hspace=hspace)
    ax1 = plt.subplot(gs[0])
    ax2 = plt.subplot(gs[1])
    ax3 = plt.subplot(gs[2])

    ax1.tick_params(which='major',direction='in',width=2,length=13,right=True,top=True,pad=7)
    ax1.tick_params(which='minor',direction='in',width=1,length=10,right=True,top=True)

    ax2.tick_params(which='major',direction='in',width=2,length=13,right=True,top=True,pad=7)
    ax2.tick_params(which='minor',direction='in',width=1,length=10,right=True,top=True)

    ax3.tick_params(which='major',direction='in',width=2,length=13,right=True,top=True,pad=7)
    ax3.tick_params(which='minor',direction='in',width=1,length=10,right=True,top=True)

    ax1.set_xlabel(xlab1,fontsize=lfs)
    ax1.set_ylabel(ylab1,fontsize=lfs)

    ax2.set_xlabel(xlab2,fontsize=lfs)
    ax2.set_ylabel(ylab2,fontsize=lfs)

    ax3.set_xlabel(xlab3,fontsize=lfs)
    ax3.set_ylabel(ylab3,fontsize=lfs)

    if Grid:
        ax1.grid()
        ax2.grid()
        ax3.grid()
    return fig,ax1,ax2,ax3
#==============================================================================#

#==============================================================================#
# Colormap stuff:

from copy import copy
def cmap_setunderwhite(cmap):
    cmap_c = copy(plt.get_cmap(cmap))
    cmap_c.set_under('white', 1.0)
    return cmap_c

def reverse_colourmap(cmap, name = 'my_cmap_r'):
    reverse = []
    k = []

    for key in cmap._segmentdata:
        k.append(key)
        channel = cmap._segmentdata[key]
        data = []

        for t in channel:
            data.append((1-t[0],t[2],t[1]))
        reverse.append(sorted(data))

    LinearL = dict(zip(k,reverse))
    my_cmap_r = mpl.colors.LinearSegmentedColormap(name, LinearL)
    return my_cmap_r

def col_alpha(col,alpha=0.1):
    rgb = colors.colorConverter.to_rgb(col)
    bg_rgb = [1,1,1]
    return [alpha * c1 + (1 - alpha) * c2
            for (c1, c2) in zip(rgb, bg_rgb)]
#==============================================================================#

#==============================================================================#
# Set up axion plot as desired:
def AxionPlotSetup(CAST_text_on=True,Haloscopes_text_on=True,QCD_text_on=True,\
                        StellarBounds_text_on=True):
    fig,ax = AxionPhoton.FigSetup(Shape='Square',\
                    ylab='$|g_{a\gamma}|$ [GeV$^{-1}$]',mathpazo=True,\
                    g_min=1e-12,g_max=1e-9,m_min=1e-6,m_max=1e0,tfs=35)

    AxionPhoton.QCDAxion(ax,text_on=False,thick_lines=True)
    AxionPhoton.StellarBounds(ax,text_on=False)
    AxionPhoton.Helioscopes(ax,projection=False,text_on=False,line_alpha=1)

    m_min=1e-6
    m_max=1e0
    plt.plot([m_min,m_max],[2e-10*0.1*m_min,2e-10*0.1*m_max],'-',\
                        color='brown',lw=1.5,zorder=0,alpha=0.1)
    plt.plot([m_min,m_max],[2e-10*10*m_min,2e-10*10*m_max],'-',\
                        color='brown',lw=1.5,zorder=0,alpha=0.1)

    fs = 25
    AxionPhoton.NeutronStars(ax,col=[0.0, 0.66, 0.42],fs=15,\
                            RescaleByMass=False,text_on=False)
    AxionPhoton.ADMX(ax,col='darkred',fs=fs,text_on=False)
    AxionPhoton.RBF_UF(ax,col='darkred',fs=fs-2,text_on=False)
    AxionPhoton.HAYSTAC(ax,col='darkred',text_on=False)
    AxionPhoton.CAPP(ax,col='darkred',fs=fs-4,text_on=False)
    AxionPhoton.ORGAN(ax,col='darkred',text_on=False)

    if CAST_text_on:
        plt.gcf().text(0.3,0.75,r'{\bf CAST}',color='k',fontsize=50)
        plt.gcf().text(0.295,0.75,r'{\bf CAST}',color='w',fontsize=50)

    if StellarBounds_text_on:
        plt.gcf().text(0.89,0.59,r'{\bf Stellar bounds}',color='k',\
                            fontsize=30,horizontalalignment='right')
        plt.gcf().text(0.888,0.59,r'{\bf Stellar bounds}',color='w',\
                            fontsize=30,horizontalalignment='right')

    if Haloscopes_text_on:
        plt.gcf().text(0.188,0.135,r'{\bf Haloscopes}',color='k',\
                        rotation=90,fontsize=32,rotation_mode='anchor',zorder=0)
        plt.gcf().text(0.185,0.1351,r'{\bf Haloscopes}',color='w',\
                        rotation=90,fontsize=32,rotation_mode='anchor',zorder=0)

    if QCD_text_on:
        cols = cm.get_cmap('YlOrBr')
        plt.gcf().text(0.561,0.13,r'{\bf KSVZ}',rotation=63,\
                                fontsize=32,rotation_mode='anchor')
        plt.gcf().text(0.611,0.13,r'{\bf DFSZ II}',rotation=63,\
                                fontsize=32,rotation_mode='anchor')
        plt.gcf().text(0.56,0.13,r'{\bf KSVZ}',rotation=63,fontsize=32,\
                                    rotation_mode='anchor',color=cols(0.9))
        plt.gcf().text(0.61,0.13,r'{\bf DFSZ II}',rotation=63,fontsize=32,\
                                    rotation_mode='anchor',color=cols(0.9))

    plt.title('space',color='w',fontsize=15)

    return fig,ax

#==============================================================================#
class AxionPhoton():
    def FigSetup(xlab=r'$m_a$ [eV]',ylab='',\
                     g_min = 1.0e-19,g_max = 1.0e-6,\
                     m_min = 1.0e-12,m_max = 1.0e7,\
                     lw=2.5,lfs=45,tfs=25,tickdir='out',\
                     Grid=False,Shape='Rectangular',mathpazo=False,TopAndRightTicks=False):

            plt.rcParams['axes.linewidth'] = lw
            plt.rc('text', usetex=True)
            plt.rc('font', family='serif',size=tfs)

            if mathpazo:
                mpl.rcParams['text.latex.preamble'] = [r'\usepackage{mathpazo}']

            if Shape=='Wide':
                fig = plt.figure(figsize=(16.5,5))
            elif Shape=='Rectangular':
                fig = plt.figure(figsize=(16.5,11))
            elif Shape=='Square':
                fig = plt.figure(figsize=(16,16))

            ax = fig.add_subplot(111)

            ax.set_xlabel(xlab,fontsize=lfs)
            ax.set_ylabel(ylab,fontsize=lfs)

            ax.tick_params(which='major',direction=tickdir,width=2.5,length=13,right=TopAndRightTicks,top=TopAndRightTicks,pad=7)
            ax.tick_params(which='minor',direction=tickdir,width=1,length=10,right=TopAndRightTicks,top=TopAndRightTicks)

            ax.set_yscale('log')
            ax.set_xscale('log')
            ax.set_xlim([m_min,m_max])
            ax.set_ylim([g_min,g_max])

            locmaj = mpl.ticker.LogLocator(base=10.0, subs=(1.0, ), numticks=50)
            locmin = mpl.ticker.LogLocator(base=10.0, subs=arange(2, 10)*.1,numticks=100)
            ax.xaxis.set_major_locator(locmaj)
            ax.xaxis.set_minor_locator(locmin)
            ax.xaxis.set_minor_formatter(mpl.ticker.NullFormatter())

            locmaj = mpl.ticker.LogLocator(base=10.0, subs=(1.0, ), numticks=100)
            locmin = mpl.ticker.LogLocator(base=10.0, subs=arange(2, 10)*.1,numticks=100)
            ax.yaxis.set_major_locator(locmaj)
            ax.yaxis.set_minor_locator(locmin)
            ax.yaxis.set_minor_formatter(mpl.ticker.NullFormatter())

            if Shape=='Rectangular':
                plt.xticks(rotation=20)

            if Grid:
                ax.grid(zorder=0)
            return fig,ax

    def QCDAxion(ax,coupling='Photon',
                      C_logwidth=10,KSVZ_on=True,DFSZ_on=True,
                      cmap='YlOrBr',fs=18,RescaleByMass=False,text_on=True,thick_lines=False):
        if RescaleByMass:
            rs1 = 1.0
            rs2 = 0.0
        else:
            rs1 = 0.0
            rs2 = 1.0
        ## QCD Axion band:
        g_min,g_max = ax.get_ylim()
        m_min,m_max = ax.get_xlim()

        # Mass-coupling relation
        def g_x(C_ag,m_a):
            return 2e-10*C_ag*m_a
        KSVZ = 1.92
        DFSZ = 0.75

        if rs1==0:
            # Plot Band
            n = 200
            g = logspace(log10(g_min),log10(g_max),n)
            m = logspace(log10(m_min),log10(m_max),n)
            QCD = zeros(shape=(n,n))
            for i in range(0,n):
                QCD[:,i] = norm.pdf(log10(g)-log10(g_x(1.0,m[i])),0.0,0.8)
            cols = cm.get_cmap(cmap)

            cols.set_under('w') # Set lowest color to white
            vmin = amax(QCD)/(C_logwidth/4.6)
            plt.contourf(m, g, QCD, 50,cmap=cols,vmin=vmin,vmax=0.9,zorder=0)
            plt.contourf(m, g, QCD, 50,cmap=cols,vmin=vmin,vmax=0.9,zorder=0)
            plt.contourf(m, g, QCD, 50,cmap=cols,vmin=vmin,vmax=0.9,zorder=0)

            # QCD Axion models
            rot = 45.0
            trans_angle = plt.gca().transData.transform_angles(array((rot,)),array([[0, 0]]))[0]
            m2 = array([1e-9,5e-8])
            if KSVZ_on:
                if thick_lines:
                    plt.plot(m,g_x(KSVZ,m),'-',linewidth=5,color='k',zorder=0)
                    plt.plot(m,g_x(KSVZ,m),'-',linewidth=3,color=cols(0.7),zorder=0)
                else:
                    plt.plot(m,g_x(KSVZ,m),'-',linewidth=2,color=cols(1.0),zorder=0)
                if text_on:
                    plt.text(1e-8,g_x(KSVZ,1e-8)*1.05,r'{\bf KSVZ}',fontsize=fs,rotation=trans_angle,color=cols(1.0),ha='left',va='bottom',rotation_mode='anchor')
            if DFSZ_on:
                if thick_lines:
                    plt.plot(m,g_x(DFSZ,m),'-',linewidth=5,color='k',zorder=0)
                    plt.plot(m,g_x(DFSZ,m),'-',linewidth=3,color=cols(0.7),zorder=0)
                else:
                    plt.plot(m,g_x(DFSZ,m),'-',linewidth=2,color=cols(1.0),zorder=0)
                if text_on:
                    plt.text(5e-8,g_x(DFSZ,5e-8)/1.5,r'{\bf DFSZ II}',fontsize=fs,rotation=trans_angle,color=cols(1.0),ha='left',va='top',rotation_mode='anchor')
        else:
            C_min,C_max = ax.get_ylim()
            n = 200
            C = logspace(log10(C_min),log10(C_max),n)
            m = logspace(log10(m_min),log10(m_max),n)
            QCD = zeros(shape=(n,n))
            for i in range(0,n):
                QCD[:,i] = norm.pdf(log10(C),0.0,0.8)
            cols = cm.get_cmap(cmap)
            cols.set_under('w') # Set lowest color to white
            vmin = amax(QCD)/(C_logwidth/2)
            plt.contourf(m, C, QCD, 50,cmap=cols,vmin=vmin,vmax=0.9,zorder=0)
            plt.contourf(m, C, QCD, 50,cmap=cols,vmin=vmin,vmax=0.9,zorder=0)
            plt.contourf(m, C, QCD, 50,cmap=cols,vmin=vmin,vmax=0.9,zorder=0)
            if thick_lines:
                plt.plot([1e-9,1e0],[0.75,0.75],'-',lw=5,color='k')
                plt.plot([1e-9,1e0],[0.75,0.75],'-',lw=3,color='k')
            else:
                plt.plot([1e-9,1e0],[0.75,0.75],'-',lw=2,color='k')
            if text_on:
                plt.text(1e-2,0.75/3,r'{\bf DFSZ II}',fontsize=fs,color='k')
        return

    def ADMX(ax,col=[0.8, 0.0, 0.0],projection=False,fs=15,RescaleByMass=False,text_on=True):
        if RescaleByMass:
            rs1 = 1.0
            rs2 = 0.0
        else:
            rs1 = 0.0
            rs2 = 1.0
        # 2018: arXiv[1804.05750]
        # 2019: arXiv[1910.08638]
        y2 = ax.get_ylim()[1]
        dat = loadtxt(limit_dir+"/AxionPhoton/ADMX.txt")
        plt.fill_between(dat[:,0],dat[:,1]/(rs1*2e-10*dat[:,0]+rs2),y2=y2,edgecolor=None,facecolor=col,zorder=0.1)
        dat = loadtxt(limit_dir+"AxionPhoton/ADMX2018.txt")
        plt.fill_between(dat[:,0],dat[:,1]/(rs1*2e-10*dat[:,0]+rs2),y2=y2,edgecolor=None,facecolor=col,zorder=0.1)
        dat = loadtxt(limit_dir+"AxionPhoton/ADMX2019_1.txt")
        plt.fill_between(dat[:,0],dat[:,1]/(rs1*2e-10*dat[:,0]+rs2),y2=y2,edgecolor=None,facecolor=col,zorder=0.1)
        dat = loadtxt(limit_dir+"/AxionPhoton/ADMX2019_2.txt")
        plt.fill_between(dat[:,0],dat[:,1]/(rs1*2e-10*dat[:,0]+rs2),y2=y2,edgecolor=None,facecolor=col,zorder=0.1)
        dat = loadtxt(limit_dir+"/AxionPhoton/ADMX_Sidecar.txt")
        plt.fill_between(dat[:,0],dat[:,1]/(rs1*2e-10*dat[:,0]+rs2),y2=y2,edgecolor=None,facecolor=col,zorder=0.1)

        if projection:
            # ADMX arXiv[1804.05750]
            dat = loadtxt(limit_dir+"/AxionPhoton/Projections/ADMX_Projected.txt")
            plt.plot(dat[:,0],dat[:,1]/(rs1*2e-10*dat[:,0]+rs2),'-',linewidth=1.5,color=col,zorder=0)
            plt.fill_between(dat[:,0],dat[:,1]/(rs1*2e-10*dat[:,0]+rs2),y2=y2,edgecolor=None,facecolor=col,zorder=0,alpha=0.1)
            if rs1==0:
                if text_on:
                    plt.text(4e-5,9e-16,r'{\bf ADMX}',fontsize=20,color=col,rotation=0,ha='left',va='top')
                plt.plot([4e-5,3e-5],[9e-16,2.1e-15],'k-',lw=1.5)
            else:
                plt.text(0.9e-6,0.15,r'{\bf ADMX}',fontsize=fs,color=col,rotation=0,ha='left',va='top')
        else:
            if rs1==0:
                if text_on:
                    plt.text(0.7e-6,1e-13,r'{\bf ADMX}',fontsize=fs,color=col,rotation=90,ha='left',va='top')

        return

    def NeutronStars(ax,col=[0.1, 0.5, 0.2],fs=15,RescaleByMass=False,text_on=True):
        # Neutron stars arXiv:[2004.00011]
        y2 = ax.get_ylim()[1]
        if RescaleByMass:
            rs1 = 1.0
            rs2 = 0.0
        else:
            rs1 = 0.0
            rs2 = 1.0
        dat = loadtxt(limit_dir+'/AxionPhoton/NeutronStars.txt')
        plt.fill_between(dat[:,0],dat[:,1]/(rs1*2e-10*dat[:,0]+rs2),y2=y2,edgecolor=None,facecolor=col,zorder=0.11)
        if rs1==0:
            if text_on:
                plt.text(5e-6,1e-12,r'{\bf Neutron stars}',fontsize=fs,color='w',ha='left')
        else:
            if text_on:
                plt.text(1e-7,4e3,r'{\bf Neutron}',fontsize=fs,color=col,ha='center')
                plt.text(1e-7,1e3,r'{\bf stars}',fontsize=fs,color=col,ha='center')
            plt.plot([3.5e-7,7e-6],[6e3,2e4],lw=1.5,color=col)

    def RBF_UF(ax,col ='darkred',fs=13,RescaleByMass=False,text_on=True):
        # UF: Phys. Rev. D42, 1297 (1990).
        # RBF: Phys. Rev. Lett. 59, 839 (1987).
        if RescaleByMass:
            rs1 = 1.0
            rs2 = 0.0
        else:
            rs1 = 0.0
            rs2 = 1.0
        y2 = ax.get_ylim()[1]
        dat = loadtxt(limit_dir+"/AxionPhoton/RBF_UF_Haloscopes.txt")
        plt.fill_between(dat[:,0],dat[:,1]/(rs1*2e-10*dat[:,0]+rs2),y2=y2,edgecolor=None,facecolor=col,zorder=0.1)
        if rs1==0:
            if text_on:
                plt.text(0.4e-5,3e-11,r'{\bf RBF+UF}',fontsize=fs,color='w',rotation=-90,ha='left',va='top')
        else:
            if text_on:
                plt.text(0.7e-5,4e3,r'{\bf RBF}',fontsize=fs,color='w',rotation=0,ha='center',va='top')
                plt.text(0.7e-5,1e3,r'{\bf UF}',fontsize=fs,color='w',rotation=0,ha='center',va='top')

        return

    def HAYSTAC(ax,col=[0.88, 0.07, 0.37],fs=13,RescaleByMass=False,projection=True,text_on=True):
        # HAYSTAC arXiv:[1803.03690]
        if RescaleByMass:
            rs1 = 1.0
            rs2 = 0.0
            zo = 3
        else:
            rs1 = 0.0
            rs2 = 1.0
            zo = 0
        y2 = ax.get_ylim()[1]
        dat = loadtxt(limit_dir+"/AxionPhoton/HAYSTAC.txt")
        if rs1==0:
            plt.plot([dat[0,0],dat[0,0]],[dat[0,1]/(rs1*2e-10*dat[0,0]+rs2),y2/(rs1*2e-10*dat[0,0]+rs2)],color=col,zorder=zo,lw=3)
            if text_on:
                if projection:
                    plt.text(2.4e-5,4e-12,r'{\bf HAYSTAC}',fontsize=fs,color=col,rotation=-90,ha='left',va='top')
                else:
                    plt.text(2.4e-5,5e-13,r'{\bf HAYSTAC}',fontsize=fs,color=col,rotation=-90,ha='left',va='top')
        else:
            plt.plot([dat[0,0],dat[0,0]],[dat[0,1]/(rs1*2e-10*dat[0,0]+rs2),y2/(rs1*2e-10*dat[0,0]+rs2)],color='k',zorder=zo,lw=4)
            plt.plot([dat[0,0],dat[0,0]],[dat[0,1]/(rs1*2e-10*dat[0,0]+rs2),y2/(rs1*2e-10*dat[0,0]+rs2)],color=col,zorder=zo,lw=3)
            if text_on:
                plt.text(dat[0,0],y2*1.2,r'{\bf HAYSTAC}',fontsize=fs,color=col,rotation=40,ha='left',rotation_mode='anchor')
            plt.plot(dat[0,0],dat[0,1]/(rs1*2e-10*dat[0,0]+rs2),'.',markersize=15,color=col,markeredgecolor='k',zorder=zo)
        return

    def CAPP(ax,col=[1, 0.1, 0.37],fs=15,RescaleByMass=False,text_on=True):
        # CAPP arXiv:[2001.05102]
        y2 = ax.get_ylim()[1]
        if RescaleByMass:
            rs1 = 1.0
            rs2 = 0.0
            zo = 3
        else:
            rs1 = 0.0
            rs2 = 1.0
            zo = 0
        dat = loadtxt(limit_dir+"/AxionPhoton/CAPP-8TB.txt")
        if rs1==0:
            plt.plot([dat[0,0],dat[0,0]],[dat[0,1]/(rs1*2e-10*dat[0,0]+rs2),y2/(rs1*2e-10*dat[0,0]+rs2)],color=col,zorder=zo,lw=3)
            if text_on:
                plt.text(1e-5,1e-13,r'{\bf CAPP}',fontsize=fs,color=col,rotation=-90,ha='center',va='top')
        else:
            plt.plot([dat[0,0],dat[0,0]],[dat[0,1]/(rs1*2e-10*dat[0,0]+rs2),y2/(rs1*2e-10*dat[0,0]+rs2)],color='k',zorder=zo,lw=4)
            plt.plot([dat[0,0],dat[0,0]],[dat[0,1]/(rs1*2e-10*dat[0,0]+rs2),y2/(rs1*2e-10*dat[0,0]+rs2)],color=col,zorder=zo,lw=3)
            if text_on:
                plt.text(dat[0,0],y2*1.8,r'{\bf CAPP}',fontsize=fs,color=col,rotation=40,ha='left',va='top',rotation_mode='anchor')
            plt.plot(dat[0,0],dat[0,1]/(rs1*2e-10*dat[0,0]+rs2),'.',markersize=15,color=col,markeredgecolor='k',zorder=zo)
        return

    def QUAX(ax,col='crimson',fs=15,RescaleByMass=False,text_on=True):
        # QUAX arXiv:[1903.06547]
        if RescaleByMass:
            rs1 = 1.0
            rs2 = 0.0
            zo = 3
        else:
            rs1 = 0.0
            rs2 = 1.0
            zo = 0
        y2 = ax.get_ylim()[1]
        dat = loadtxt(limit_dir+"/AxionPhoton/QUAX.txt")

        if rs1==0:
            plt.plot([dat[0,0],dat[0,0]],[dat[0,1]/(rs1*2e-10*dat[0,0]+rs2),y2/(rs1*2e-10*dat[0,0]+rs2)],color=col,lw=2,zorder=zo)
            if text_on:
                plt.text(5.2e-5,4e-11,r'{\bf QUAX}',fontsize=fs,color=col,rotation=-90,ha='center',va='top')
        else:
            plt.plot([dat[0,0],dat[0,0]],[dat[0,1]/(rs1*2e-10*dat[0,0]+rs2),y2/(rs1*2e-10*dat[0,0]+rs2)],color='k',lw=4,zorder=zo)
            plt.plot([dat[0,0],dat[0,0]],[dat[0,1]/(rs1*2e-10*dat[0,0]+rs2),y2/(rs1*2e-10*dat[0,0]+rs2)],color=col,lw=3,zorder=zo)
            if text_on:
                plt.text(dat[0,0]*1.3,y2*1.2,r'{\bf QUAX}',fontsize=fs,color=col,rotation=40,ha='left',rotation_mode='anchor')
            plt.plot(dat[0,0],dat[0,1]/(rs1*2e-10*dat[0,0]+rs2),'.',markersize=15,color=col,markeredgecolor='k',zorder=zo)
        return

    def ABRACADABRA(ax,col=[0.83, 0.07, 0.37],fs=15,projection=False,RescaleByMass=False,text_on=True):
        # ABRACADABRA arXiv:[1810.12257]
        if RescaleByMass:
            rs1 = 1.0
            rs2 = 0.0
        else:
            rs1 = 0.0
            rs2 = 1.0
        y2 = ax.get_ylim()[1]
        dat = loadtxt(limit_dir+"/AxionPhoton/ABRACADABRA.txt")
        n = shape(dat)[0]
        plt.fill_between(dat[:,0],dat[:,1]/(rs1*2e-10*dat[:,0]+rs2),y2=y2,edgecolor=None,facecolor=col,zorder=2)
        x = dat[arange(0,n,20),0]
        y = dat[arange(0,n,20),1]
        y[-1] = y2
        plt.plot(x,y/(rs1*2e-10*x+rs2),'k-',lw=1,zorder=10,alpha=0.5)
        if text_on:
            if rs1==0:
                plt.text(1.5e-9,3e-8,r'{\bf ABRA}',fontsize=fs,color='w',rotation=0,ha='center',va='top',zorder=10)
                plt.text(1.5e-9,1e-8,r'10 cm',fontsize=fs,color='w',rotation=0,ha='center',va='top',zorder=10)

        if projection:
            dat = loadtxt(limit_dir+"/AxionPhoton/Projections/ABRACADABRA.txt")
            plt.plot(dat[:,0],dat[:,1]/(rs1*2e-10*dat[:,0]+rs2),'-',linewidth=1.5,color=col,zorder=0)
            plt.fill_between(dat[:,0],dat[:,1]/(rs1*2e-10*dat[:,0]+rs2),y2=y2,edgecolor=None,facecolor=col,zorder=0,alpha=0.1)
            if rs1==0:
                if text_on:
                    plt.text(1e-12,2.5e-18,r'{\bf ABRACADABRA}',fontsize=fs-1,color=col,rotation=13,ha='left',va='top')
            else:
                if text_on:
                    plt.text(1.3e-9,1.0e2,r'{\bf ABRACADABRA}',fontsize=fs-1,color=col,rotation=0,ha='left',va='top')
                plt.plot([dat[-1,0],dat[-1,0]],[dat[-1,1]/(rs1*2e-10*dat[-1,0]+rs2),1e6],lw=1.5,color=col,zorder=0)
        return

    def ORGAN(ax,col='crimson',projection=False,fs=15,RescaleByMass=False,text_on=True):
        # ORGAN arXiv[1706.00209]
        y2 = ax.get_ylim()[1]
        if RescaleByMass:
            rs1 = 1.0
            rs2 = 0.0
        else:
            rs1 = 0.0
            rs2 = 1.0
        col = [0.8, 0.0, 0.0]
        dat = loadtxt(limit_dir+"/AxionPhoton/ORGAN.txt")
        plt.fill_between(dat[:,0],dat[:,1]/(rs1*2e-10*dat[:,0]+rs2),y2=y2,edgecolor=col,facecolor=col,zorder=0.1,lw=2)

        if projection:
            dat = loadtxt(limit_dir+"/AxionPhoton/Projections/ORGAN_Projected.txt")
            plt.plot(dat[:,0],dat[:,1]/(rs1*2e-10*dat[:,0]+rs2),'-',linewidth=1.5,color=col,zorder=0)
            plt.fill_between(dat[:,0],dat[:,1]/(rs1*2e-10*dat[:,0]+rs2),y2=y2,edgecolor=None,facecolor=col,zorder=0,alpha=0.1)
            if rs1==0:
                if text_on:
                    plt.text(5e-4,1.15e-14,r'{\bf ORGAN}',fontsize=18,color=col,rotation=0,ha='left',va='top')
                plt.plot([5e-4,1.5e-4],[1.3e-14,6e-13],'k-',lw=1.5)
            else:
                plt.text(1.2e-4,1e3,r'{\bf ORGAN}',fontsize=18,color=col,rotation=-90,ha='left',va='top')

        else:
            if rs1==0:
                if text_on:
                    plt.text(110e-6,6e-11,r'{\bf ORGAN}',fontsize=fs,color=col,rotation=-90,ha='left',va='top')
        return

    def MADMAX(ax,col=[0.6, 0.1, 0.1],fs=18,RescaleByMass=False,text_on=True):
        # MADMAX arXiv[2003.10894]
        y2 = ax.get_ylim()[1]
        if RescaleByMass:
            rs1 = 1.0
            rs2 = 0.0
        else:
            rs1 = 0.0
            rs2 = 1.0
        dat = loadtxt(limit_dir+"/AxionPhoton/Projections/MADMAX.txt")
        plt.plot(dat[:,0],dat[:,1]/(rs1*2e-10*dat[:,0]+rs2),'-',linewidth=1.5,color=col,zorder=0)
        plt.fill_between(dat[:,0],dat[:,1]/(rs1*2e-10*dat[:,0]+rs2),y2=y2,edgecolor=None,facecolor=col,zorder=0,alpha=0.2)
        if rs1==0:
            plt.text(1.5e-4,3.5e-15,r'{\bf MADMAX}',fontsize=18,color=col,rotation=0,ha='left',va='top')
            plt.plot([3e-4,1.3e-4],[4.5e-15,2.6e-14],'k-',lw=1.5)
        else:
            if text_on:
                plt.text(4e-5,3.5e-1,r'{\bf MADMAX}',fontsize=fs,color=col,rotation=0,ha='left',va='top')
        return

    def KLASH(ax,col=[0.6, 0.1, 0.2],fs=15,RescaleByMass=False,text_on=True):
        # KLASH arXiv:[1707.06010]
        y2 = ax.get_ylim()[1]
        if RescaleByMass:
            rs1 = 1.0
            rs2 = 0.0
        else:
            rs1 = 0.0
            rs2 = 1.0
        dat = loadtxt(limit_dir+"/AxionPhoton/Projections/KLASH.txt")
        plt.plot(dat[:,0],dat[:,1]/(rs1*2e-10*dat[:,0]+rs2),'-',linewidth=1.5,color=col,zorder=0)
        plt.fill_between(dat[:,0],dat[:,1]/(rs1*2e-10*dat[:,0]+rs2),y2=y2,facecolor=col,zorder=0,alpha=0.3)
        if text_on:
            if rs1==0:
                plt.text(1e-7,1e-12,r'{\bf KLASH}',rotation=90,fontsize=fs,color=col,ha='left',va='top')
            else:
                plt.text(2.5e-7,1.3e0,r'{\bf KLASH}',rotation=90,fontsize=fs,color=col,ha='left',va='top',rotation_mode='anchor')

    def BRASS(ax,col=[0.5, 0.1, 0.2],fs=15,RescaleByMass=False,text_on=True):
        # BRASS http://www.iexp.uni-hamburg.de/groups/astroparticle/brass/brassweb.htm
        y2 = ax.get_ylim()[1]
        if RescaleByMass:
            rs1 = 1.0
            rs2 = 0.0
        else:
            rs1 = 0.0
            rs2 = 1.0
        dat = loadtxt(limit_dir+"/AxionPhoton/Projections/BRASS.txt")
        plt.plot(dat[:,0],dat[:,1]/(rs1*2e-10*dat[:,0]+rs2),'-',linewidth=1.5,color=col,zorder=0)
        plt.fill_between(dat[:,0],dat[:,1]/(rs1*2e-10*dat[:,0]+rs2),y2=y2,facecolor=col,zorder=0,alpha=0.1)
        if text_on:
            if rs1==0:
                plt.text(2.3e-3,0.6e-10,r'{\bf BRASS}',rotation=56,fontsize=fs,color=col,ha='left',va='top')
            else:
                plt.text(1e-3,0.12e3,r'{\bf BRASS}',rotation=15,fontsize=fs,color=col,ha='left',va='top')

    def TOORAD(ax,col=[0.8, 0.1, 0.2],fs=15,RescaleByMass=False,text_on=True):
        # TOORAD arXiv[1807.08810]
        y2 = ax.get_ylim()[1]
        if RescaleByMass:
            rs1 = 1.0
            rs2 = 0.0
        else:
            rs1 = 0.0
            rs2 = 1.0
        dat = loadtxt(limit_dir+"/AxionPhoton/Projections/TOORAD.txt")
        plt.plot(dat[:,0],dat[:,1]/(rs1*2e-10*dat[:,0]+rs2),'-',linewidth=1.5,color=col,zorder=0)
        plt.fill_between(dat[:,0],dat[:,1]/(rs1*2e-10*dat[:,0]+rs2),y2=y2,facecolor=col,zorder=0,alpha=0.1)
        if rs1==0:
            if text_on:
                plt.text(0.5e-2,5e-14,r'{\bf TOORAD}',rotation=0,fontsize=18,color=col,ha='left',va='top')
            plt.plot([0.5e-2,0.21e-2],[5e-14,1e-13],'k-',lw=1.5)
        else:
            if text_on:
                plt.text(0.6e-3,4e-1,r'{\bf TOORAD}',rotation=-25,fontsize=fs,color=col,ha='left',va='top')

    def LAMPOST(ax,col=[0.8, 0.1, 0.2],fs=15,RescaleByMass=False,text_on=True):
        # LAMPOST arXiv[1803.11455]
        y2 = ax.get_ylim()[1]
        if RescaleByMass:
            rs1 = 1.0
            rs2 = 0.0
        else:
            rs1 = 0.0
            rs2 = 1.0
        dat = loadtxt(limit_dir+"/AxionPhoton/Projections/LAMPOST.txt",delimiter=',')
        plt.plot(dat[:,0],dat[:,1]/(rs1*2e-10*dat[:,0]+rs2),'-',linewidth=1.5,color=col,zorder=0)
        plt.fill_between(dat[:,0],dat[:,1]/(rs1*2e-10*dat[:,0]+rs2),y2=y2,facecolor=col,zorder=0,alpha=0.1)

        if rs1==0:
            if text_on:
                plt.text(0.8e-1,5e-12,r'{\bf LAMPOST}',rotation=-90,fontsize=fs,color=col,ha='left',va='top')
        else:
            if text_on:
                plt.text(0.9e-1,1.9e-1,r'{\bf LAMPOST}',rotation=0,fontsize=fs,color=col,ha='left',va='top')

    # Low mass ALP haloscopes
    def DANCE(ax,col=[0.8, 0.1, 0.2],fs=15,text_on=True):
        # DANCE arXiv[1911.05196]
        y2 = ax.get_ylim()[1]
        dat = loadtxt(limit_dir+"/AxionPhoton/Projections/DANCE.txt")
        plt.plot(dat[:,0],dat[:,1],'-',linewidth=1.5,color=col,zorder=0)
        plt.fill_between(dat[:,0],dat[:,1],y2=y2,facecolor=col,zorder=0,alpha=0.1)
        if text_on:
            plt.text(1.7e-12,2e-13,r'{\bf DANCE}',rotation=50,fontsize=fs,color=col,ha='left',va='top')

    def aLIGO(ax,col=[0.8, 0.1, 0.2],fs=15,text_on=True):
        # aLIGO arXiv[1903.02017]
        y2 = ax.get_ylim()[1]
        dat = loadtxt(limit_dir+"/AxionPhoton/Projections/aLIGO.txt")
        plt.plot(dat[:,0],dat[:,1],'-',linewidth=1.5,color=col,zorder=0)
        plt.fill_between(dat[:,0],dat[:,1],y2=y2,facecolor=col,zorder=0,alpha=0.1)
        if text_on:
            plt.text(0.2e-9,0.35e-13,r'{\bf aLIGO}',rotation=0,fontsize=fs,color=col,ha='left',va='top')

    def ADBC(ax,col=[0.8, 0.1, 0.2],fs=15,text_on=True):
        # ADBC arXiv[1809.01656]
        y2 = ax.get_ylim()[1]
        dat = loadtxt(limit_dir+"/AxionPhoton/Projections/ADBC.txt")
        plt.plot(dat[:,0],dat[:,1],'-',linewidth=1.5,color=col,zorder=0)
        plt.fill_between(dat[:,0],dat[:,1],y2=y2,facecolor=col,zorder=0,alpha=0.1)
        if text_on:
            plt.text(1e-11,1.3e-12,r'{\bf ADBC}',rotation=26,fontsize=fs,color=col,ha='left',va='top')

    def SHAFT(ax,col='red',fs=16,text_on=True):
        # SHAFT arXiv:[2003.03348]
        y2 = ax.get_ylim()[1]
        dat = loadtxt(limit_dir+"/AxionPhoton/SHAFT.txt")
        n = shape(dat)[0]
        x = dat[arange(0,n,2),0]
        y = dat[arange(0,n,2),1]
        y[-1] = y2
        plt.plot(x,y,'k-',lw=1,zorder=10,alpha=0.5)
        plt.fill_between(dat[:,0],dat[:,1],y2=y2,edgecolor=None,facecolor=col,zorder=1.8)
        if text_on:
            plt.text(0.8e-10,3e-10,r'{\bf SHAFT}',fontsize=fs,color='w',rotation=0,ha='center',va='top',zorder=9)
        return

    def ALPS(ax,projection=True,col=[0.8, 0.25, 0.33],fs=15,RescaleByMass=False,text_on=True):
        # ALPS-I arXiv:[1004.1313]
        if RescaleByMass:
            rs1 = 1.0
            rs2 = 0.0
        else:
            rs1 = 0.0
            rs2 = 1.0

        y2 = ax.get_ylim()[1]
        dat = loadtxt(limit_dir+"/AxionPhoton/ALPS.txt")
        plt.plot(dat[:,0],dat[:,1],'k-',lw=2.5,zorder=1.53,alpha=0.5)
        plt.fill_between(dat[:,0],dat[:,1]/(rs1*2e-10*dat[:,0]+rs2),y2=y2,edgecolor=None,facecolor=col,zorder=1.53,lw=0.01)
        if text_on:
            plt.text(1e-5,7e-8,r'{\bf ALPS-I}',fontsize=20,color='w')
        if projection:
            dat = loadtxt(limit_dir+"/AxionPhoton/Projections/ALPS-II.txt")
            plt.plot(dat[:,0],dat[:,1]/(rs1*2e-10*dat[:,0]+rs2),'-',lw=1.5,zorder=1.5,color='k',alpha=0.5)
            if text_on:
                if RescaleByMass:
                    plt.text(9e-4,2.5e3,r'{\bf ALPS-II}',fontsize=20,color='k',rotation=20,alpha=0.5)
                else:
                    plt.text(1.5e-3,3e-9,r'{\bf ALPS-II}',rotation=60,fontsize=18,color='w',zorder=10)
        return

    def OSQAR(ax,col=[0.6, 0.2, 0.25],fs=15,text_on=True):
        # OSQAR arXiv:[]
        y2 = ax.get_ylim()[1]
        dat = loadtxt(limit_dir+"/AxionPhoton/OSQAR.txt")
        plt.plot(dat[:,0],dat[:,1],'k-',lw=2.5,zorder=1.52,alpha=0.5)
        plt.fill_between(dat[:,0],dat[:,1],y2=y2,edgecolor=None,facecolor=col,zorder=1.52,lw=0.01)
        if text_on:
            plt.text(1e-5,1.5e-8,r'{\bf OSQAR}',fontsize=17,color='w')
        return

    def PVLAS(ax,col=[0.4, 0.2, 0.2],fs=15,text_on=True):
        # PVLAS arXiv:[]
        y2 = ax.get_ylim()[1]
        dat = loadtxt(limit_dir+"/AxionPhoton/PVLAS.txt")
        plt.plot(dat[:,0],dat[:,1],'k-',lw=2.5,zorder=1.51,alpha=0.4)
        plt.fill_between(dat[:,0],dat[:,1],y2=y2,edgecolor=None,facecolor=col,zorder=1.51,lw=0.01)
        if text_on:
            plt.text(2e-3,9e-8,r'{\bf PVLAS}',fontsize=17,color='w',rotation=45)
        return

    def CROWS(ax,col=[0.7, 0.2, 0.2],fs=15,text_on=True):
        # CROWS arXiv:[1310.8098]
        y2 = ax.get_ylim()[1]
        dat = loadtxt(limit_dir+"/AxionPhoton/CROWS.txt")
        plt.plot(dat[:,0],dat[:,1],'k-',lw=2.5,zorder=1.54,alpha=0.4)
        plt.fill_between(dat[:,0],dat[:,1],y2=y2,edgecolor=None,facecolor=col,zorder=1.54,lw=0.01)
        if text_on:
            plt.text(1e-7,1.5e-7,r'{\bf CROWS}',fontsize=17,color='w',rotation=0)
        return

    ####################################################
    def Helioscopes(ax,col=[0.5, 0.0, 0.13],fs=25,projection=True,RescaleByMass=False,text_on=True,line_alpha=0.5):
        # CAST arXiv:[1705.02290]
        y2 = ax.get_ylim()[1]
        if RescaleByMass:
            rs1 = 1.0
            rs2 = 0.0
        else:
            rs1 = 0.0
            rs2 = 1.0
        dat = loadtxt(limit_dir+"/AxionPhoton/CAST_highm.txt")
        plt.plot(dat[:,0],dat[:,1]/(rs1*2e-10*dat[:,0]+rs2),'k-',lw=2,zorder=1.49,alpha=line_alpha)
        plt.fill_between(dat[:,0],dat[:,1]/(rs1*2e-10*dat[:,0]+rs2),y2=y2,edgecolor='k',facecolor=col,zorder=1.49,lw=0.1)
        dat = loadtxt(limit_dir+"/AxionPhoton/CAST.txt")
        plt.plot(dat[:,0],dat[:,1]/(rs1*2e-10*dat[:,0]+rs2),'k-',lw=2,zorder=1.5,alpha=line_alpha)
        plt.fill_between(dat[:,0],dat[:,1]/(rs1*2e-10*dat[:,0]+rs2),y2=y2,edgecolor='k',facecolor=col,zorder=1.5,lw=0.1)
        if text_on:
            if rs1==0:
                plt.text(1e-1,1.5e-9,r'{\bf CAST}',fontsize=fs+4,color='w',rotation=0,ha='center',va='top')
            else:
                plt.text(4e-2,5e3,r'{\bf CAST}',fontsize=fs+4,color='w',rotation=0,ha='center',va='top')

        if projection:
            # IAXO arXiv[1212.4633]
            IAXO_col = 'purple'
            IAXO = loadtxt(limit_dir+"/AxionPhoton/Projections/IAXO.txt")
            plt.plot(IAXO[:,0],IAXO[:,1]/(rs1*2e-10*IAXO[:,0]+rs2),'--',linewidth=2.5,color=IAXO_col,zorder=0.5)
            plt.fill_between(IAXO[:,0],IAXO[:,1]/(rs1*2e-10*IAXO[:,0]+rs2),y2=y2,edgecolor=None,facecolor=IAXO_col,zorder=0,alpha=0.3)
            if text_on:
                if rs1==0:
                    plt.text(0.35e-1,0.2e-11,r'{\bf IAXO}',fontsize=fs,color=IAXO_col,rotation=45)
                else:
                    plt.text(0.7e-2,0.15e1,r'{\bf IAXO}',fontsize=fs,color=IAXO_col,rotation=-20)
        return

    def Haloscopes(ax,projection=True,fs=20,text_on=True):
        AxionPhoton.ADMX(ax,projection=projection,fs=fs,text_on=text_on)
        AxionPhoton.RBF_UF(ax,fs=fs-2,text_on=text_on)
        AxionPhoton.HAYSTAC(ax,projection=projection,text_on=text_on)
        AxionPhoton.ABRACADABRA(ax,fs=fs,projection=projection,text_on=text_on)
        AxionPhoton.SHAFT(ax,text_on=text_on)
        AxionPhoton.CAPP(ax,fs=fs-4,text_on=text_on)
        AxionPhoton.ORGAN(ax,projection=projection,text_on=text_on)

        if projection:
            AxionPhoton.MADMAX(ax,text_on=text_on)
            AxionPhoton.KLASH(ax,text_on=text_on)
            AxionPhoton.TOORAD(ax,text_on=text_on)
            AxionPhoton.BRASS(ax,text_on=text_on)
            AxionPhoton.ADBC(ax,text_on=text_on)
            AxionPhoton.DANCE(ax,text_on=text_on)
            AxionPhoton.aLIGO(ax,text_on=text_on)
        else:
            AxionPhoton.QUAX(ax,text_on=text_on)
        return

    def LSW(ax,projection=True,text_on=True):
        AxionPhoton.PVLAS(ax,text_on=text_on)
        AxionPhoton.ALPS(ax,projection=projection,text_on=text_on)
        AxionPhoton.OSQAR(ax,text_on=text_on)
        AxionPhoton.CROWS(ax,text_on=text_on)
        return

    def AstroBounds(ax,projection=True,fs=15,text_on=True):
        y2 = ax.get_ylim()[1]
        ### Astrophysical constraints

        # SN-gamma rays arXiv:[1410.3747]
        SNgamma_col = [0.05, 0.5, 0.06]
        SNgamma = loadtxt(limit_dir+"/AxionPhoton/SN-gamma.txt")
        plt.plot(SNgamma[:,0],SNgamma[:,1],'k-',alpha=0.6,zorder=0.21,lw=2)
        plt.fill_between(SNgamma[:,0],SNgamma[:,1],y2=y2,edgecolor=None,facecolor=SNgamma_col,zorder=0.21)

        # M87 Limits from arXiv:[1703.07354]
        M87_col = [0.0, 0.66, 0.42]
        M87 = loadtxt(limit_dir+"/AxionPhoton/M87.txt")
        plt.plot(M87[:,0],M87[:,1],'k-',lw=2,alpha=0.8,zorder=0.2)
        plt.fill_between(M87[:,0],M87[:,1],y2=y2,edgecolor=None,facecolor=M87_col,zorder=0.2)

        # HYDRA-A arXiv:[1304.0989]
        HYDRA_col = [0.24, 0.71, 0.54]
        HYDRA = loadtxt(limit_dir+"/AxionPhoton/HYDRA_A.txt")
        plt.plot(HYDRA[:,0],HYDRA[:,1],'k-',alpha=0.6,zorder=0.23,lw=2)
        plt.fill_between(HYDRA[:,0],HYDRA[:,1],y2=y2,edgecolor=None,facecolor=HYDRA_col,zorder=0.23)

        # HESS arXiv:[1304.0700]
        HESS_col = [0.0, 0.55, 0.3]
        HESS = loadtxt(limit_dir+"/AxionPhoton/HESS.txt")
        plt.plot(HESS[:,0],HESS[:,1],'k-',alpha=0.6,zorder=0.2,lw=2)
        plt.fill_between(HESS[:,0],HESS[:,1],y2=y2,edgecolor=None,facecolor=HESS_col,zorder=0.2)

        # Fermi NGC1275 arXiv:[1603.06978]
        Fermi_col = [0.0, 0.42, 0.24]
        Fermi1 = loadtxt(limit_dir+"/AxionPhoton/Fermi1.txt")
        Fermi2 = loadtxt(limit_dir+"/AxionPhoton/Fermi2.txt")
        plt.fill_between(Fermi1[:,0],Fermi1[:,1],y2=y2,edgecolor=Fermi_col,facecolor=Fermi_col,zorder=0.24,lw=3)
        plt.fill(Fermi2[:,0],1.01*Fermi2[:,1],edgecolor=Fermi_col,facecolor=Fermi_col,lw=3,zorder=0.24)
        Fermi1 = loadtxt(limit_dir+"/AxionPhoton/Fermi_bound.txt")
        Fermi2 = loadtxt(limit_dir+"/AxionPhoton/Fermi_hole.txt")
        plt.plot(Fermi1[:,0],Fermi1[:,1],'k-',alpha=0.5,lw=1.5,zorder=0.24)
        plt.plot(Fermi2[:,0],Fermi2[:,1],'k-',alpha=0.5,lw=1.5,zorder=0.24)

        # Optical telescope [astro-ph/0611502]
        Telescopes_col = [0.09, 0.45, 0.27]
        Telescopes = loadtxt(limit_dir+"/AxionPhoton/Telescopes.txt")
        plt.fill_between(Telescopes[:,0],Telescopes[:,1],y2=y2,edgecolor=None,facecolor=Telescopes_col,zorder=0.2)
        if text_on:
            plt.text(3.3,4e-12,r'{\bf Telescopes}',fontsize=fs,color=Telescopes_col,rotation=-90,ha='left',va='top')
            plt.text(4.8e-10,1.2e-11,r'{\bf Fermi}',fontsize=fs,color='w',ha='left',va='top')
            plt.text(2e-8,1.6e-11,r'{\bf HESS}',fontsize=fs+1,color=HESS_col,ha='left',va='top')
            plt.text(1.5e-12,4e-11,r'{\bf Hydra}',fontsize=fs-2,color='w',ha='left',va='top')
            plt.text(3e-12,2e-11,r'\quad {\bf A}',fontsize=fs-2,color='w',ha='left',va='top')
            plt.text(1.4e-12,4e-12,r'\quad {\bf M87}',fontsize=fs,color='w',ha='left',va='top')
            plt.text(3e-11,2e-11,r'{\bf SN}-$\gamma$',fontsize=fs,color='w',ha='left',va='top')

        # Chandra arXiv:[1907.05475]
        Chandra_col = [0.0, 0.3, 0.24]
        Chandra = loadtxt(limit_dir+'/AxionPhoton/Chandra.txt')
        plt.plot(Chandra[:,0],Chandra[:,1],'k-',alpha=0.8,lw=2,zorder=0.1)
        plt.fill_between(Chandra[:,0],Chandra[:,1],y2=y2,edgecolor=None,facecolor=Chandra_col,zorder=0.1)
        if text_on:
            if projection==False:
                plt.text(1.1e-11,2e-12,r'{\bf Chandra}',fontsize=fs,color=Chandra_col,rotation=0,ha='left',va='top')

        if projection==True:
            # Fermi nearby SN prospects arXiv:[1609.02350]
            FermiSN = loadtxt(limit_dir+"/AxionPhoton/Projections/FermiSN.txt")
            plt.fill_between(FermiSN[:,0],FermiSN[:,1],y2=y2,edgecolor=Fermi_col,linewidth=1.5,facecolor=Fermi_col,zorder=0.1,alpha=0.2)
            if text_on:
                plt.text(1e-9,4e-12,r'{\bf Fermi SN}',fontsize=fs,color=Fermi_col,rotation=43,ha='left',va='top')

        return

    def Cosmology(ax,fs=30,text_on=True):
        y2 = ax.get_ylim()[1]
        ## Cosmology constraints see arXiv:[1210.3196] for summary
        # Xray Background
        XRAY_col = [0.03, 0.57, 0.82]
        XRAY = loadtxt(limit_dir+"/AxionPhoton/XRAY.txt")
        plt.plot(XRAY[:,0],XRAY[:,1],color='k',alpha=0.5,zorder=0.3,lw=2)
        plt.fill_between(XRAY[:,0],XRAY[:,1],y2=1e-11,edgecolor=None,facecolor=XRAY_col,zorder=0.3)

        # Extragalactic background light
        EBL_col =  [0.0, 0.2, 0.6]
        EBL = loadtxt(limit_dir+"/AxionPhoton/EBL.txt")
        EBL2 = loadtxt(limit_dir+"/AxionPhoton/EBL2.txt")
        plt.plot(EBL[:,0],EBL[:,1],'k',lw=2.5,zorder=0.4,alpha=0.8)
        plt.fill_between(EBL[:,0],EBL[:,1],y2=y2,edgecolor=None,facecolor=EBL_col,zorder=0.5)
        plt.fill_between(EBL2[:,0],EBL2[:,1],y2=y2,edgecolor=None,facecolor=EBL_col,zorder=0.5)

        # Ionisation fraction
        x_ion_col = [0.27, 0.51, 0.71]
        x_ion = loadtxt(limit_dir+"/AxionPhoton/x_ion.txt")
        plt.plot(x_ion[:,0],x_ion[:,1],'k',lw=2.5,zorder=0.4,alpha=0.8)
        plt.fill_between(x_ion[:,0],x_ion[:,1],y2=y2,edgecolor=None,facecolor=x_ion_col,zorder=0.5)

        # BBN+N_eff arXiv:[2002.08370]
        BBN_col = [0.27, 0.51, 0.71]
        BBN = loadtxt(limit_dir+"/AxionPhoton/BBN_Neff.txt")
        plt.plot(BBN[:,0],BBN[:,1],'k',lw=2,zorder=0.4,alpha=0.5)
        plt.fill_between(BBN[:,0],BBN[:,1],y2=y2,edgecolor=None,facecolor=BBN_col,zorder=0.4)

        if text_on:
            plt.text(3e3,0.8e-16,r'{\bf X-rays}',fontsize=fs,color='w',rotation=-50,ha='left',va='top')
            plt.text(1e4,5e-14,r'{\bf EBL}',fontsize=fs+5,color='w',rotation=-55,ha='left',va='top')
            plt.text(100.5744,5.1720e-11,r'{\bf Ionisation}',fontsize=fs-6,color='w',rotation=-90,ha='left',va='top')
            plt.text(40,4.1720e-11,r'{\bf fraction}',fontsize=fs-6,color='w',rotation=-90,ha='left',va='top')
            plt.text(0.05e6,8e-10,r'{\bf BBN}+$N_{\rm eff}$',fontsize=fs,color='w',rotation=-55,ha='left',va='top')

    def StellarBounds(ax,fs=30,text_on=True):
        y2 = ax.get_ylim()[1]
        # Stellar physics constraints

        # Globular clusters arXiv:[1406.6053]
        HB_col = [0.0, 0.66, 0.42]
        HB = loadtxt(limit_dir+"/AxionPhoton/HorizontalBranch.txt")
        plt.plot(HB[:,0],HB[:,1],color='k',alpha=0.5,zorder=1,lw=2)
        plt.fill_between(HB[:,0],HB[:,1],y2=y2,edgecolor=None,facecolor=HB_col,zorder=1)

        # Solar neutrino B8 bound arXiv:[1501.01639]
        SolarNu_col = [0.01, 0.75, 0.24]
        SolarNu = loadtxt(limit_dir+"/AxionPhoton/SolarNu.txt")
        plt.plot(SolarNu[:,0],SolarNu[:,1],color='k',lw=2,alpha=0.5,zorder=1)
        plt.fill_between(SolarNu[:,0],SolarNu[:,1],y2=y2,edgecolor=None,facecolor=SolarNu_col,zorder=1)

        # SN1987A-neutrinos updated arXiv:[1808.10136]
        SN = loadtxt(limit_dir+"/AxionPhoton/SN1987A_2019.txt")
        plt.fill_between(SN[:,0],SN[:,1],y2=y2,edgecolor=None,facecolor='ForestGreen',zorder=0.1)
        # SN1987A-decay arXiv:[1702.02964]
        SN = loadtxt(limit_dir+"/AxionPhoton/SN1987A_decay.txt")
        plt.fill_between(SN[:,0],SN[:,1],y2=y2,edgecolor=None,facecolor='ForestGreen',zorder=0.1)

        if text_on:
            plt.text(0.4e6,6e-7,r'{\bf SN1987A}',fontsize=fs-9,color='w',rotation=-60,ha='left',va='top')
            plt.text(1e1,3e-9,r'{\bf Solar} $\nu$',fontsize=fs+3,color='w')
            plt.text(1.4e0,1.5e-10,r'{\bf Horizontal branch}',fontsize=fs-7,color='w')

#==============================================================================#

from matplotlib import patches
from matplotlib import text as mtext
import numpy as np
import math

class CurvedText(mtext.Text):
    """
    A text object that follows an arbitrary curve.
    """
    def __init__(self, x, y, text, axes, **kwargs):
        super(CurvedText, self).__init__(x[0],y[0],' ', **kwargs)

        axes.add_artist(self)

        ##saving the curve:
        self.__x = x
        self.__y = y
        self.__zorder = self.get_zorder()

        ##creating the text objects
        self.__Characters = []
        for c in text:
            if c == ' ':
                ##make this an invisible 'a':
                t = mtext.Text(0,0,'a')
                t.set_alpha(0.0)
            else:
                t = mtext.Text(0,0,c, **kwargs)

            #resetting unnecessary arguments
            t.set_ha('center')
            t.set_rotation(0)
            t.set_zorder(self.__zorder +1)

            self.__Characters.append((c,t))
            axes.add_artist(t)

    ##overloading some member functions, to assure correct functionality
    ##on update
    def set_zorder(self, zorder):
        super(CurvedText, self).set_zorder(zorder)
        self.__zorder = self.get_zorder()
        for c,t in self.__Characters:
            t.set_zorder(self.__zorder+1)

    def draw(self, renderer, *args, **kwargs):
        """
        Overload of the Text.draw() function. Do not do
        do any drawing, but update the positions and rotation
        angles of self.__Characters.
        """
        self.update_positions(renderer)

    def update_positions(self,renderer):
        """
        Update positions and rotations of the individual text elements.
        """

        #preparations

        ##determining the aspect ratio:
        ##from https://stackoverflow.com/a/42014041/2454357

        ##data limits
        xlim = self.axes.get_xlim()
        ylim = self.axes.get_ylim()
        ## Axis size on figure
        figW, figH = self.axes.get_figure().get_size_inches()
        ## Ratio of display units
        _, _, w, h = self.axes.get_position().bounds
        ##final aspect ratio
        aspect = ((figW * w)/(figH * h))*(ylim[1]-ylim[0])/(xlim[1]-xlim[0])

        #points of the curve in figure coordinates:
        x_fig,y_fig = (
            np.array(l) for l in zip(*self.axes.transData.transform([
            (i,j) for i,j in zip(self.__x,self.__y)
            ]))
        )

        #point distances in figure coordinates
        x_fig_dist = (x_fig[1:]-x_fig[:-1])
        y_fig_dist = (y_fig[1:]-y_fig[:-1])
        r_fig_dist = np.sqrt(x_fig_dist**2+y_fig_dist**2)

        #arc length in figure coordinates
        l_fig = np.insert(np.cumsum(r_fig_dist),0,0)

        #angles in figure coordinates
        rads = np.arctan2((y_fig[1:] - y_fig[:-1]),(x_fig[1:] - x_fig[:-1]))
        degs = np.rad2deg(rads)

        rel_pos = 10
        for c,t in self.__Characters:
            #finding the width of c:
            t.set_rotation(0)
            t.set_va('center')
            bbox1  = t.get_window_extent(renderer=renderer)
            w = bbox1.width
            h = bbox1.height

            #ignore all letters that don't fit:
            if rel_pos+w/2 > l_fig[-1]:
                t.set_alpha(0.0)
                rel_pos += w
                continue

            elif c != ' ':
                t.set_alpha(1.0)

            #finding the two data points between which the horizontal
            #center point of the character will be situated
            #left and right indices:
            il = np.where(rel_pos+w/2 >= l_fig)[0][-1]
            ir = np.where(rel_pos+w/2 <= l_fig)[0][0]

            #if we exactly hit a data point:
            if ir == il:
                ir += 1

            #how much of the letter width was needed to find il:
            used = l_fig[il]-rel_pos
            rel_pos = l_fig[il]

            #relative distance between il and ir where the center
            #of the character will be
            fraction = (w/2-used)/r_fig_dist[il]

            ##setting the character position in data coordinates:
            ##interpolate between the two points:
            x = self.__x[il]+fraction*(self.__x[ir]-self.__x[il])
            y = self.__y[il]+fraction*(self.__y[ir]-self.__y[il])

            #getting the offset when setting correct vertical alignment
            #in data coordinates
            t.set_va(self.get_va())
            bbox2  = t.get_window_extent(renderer=renderer)

            bbox1d = self.axes.transData.inverted().transform(bbox1)
            bbox2d = self.axes.transData.inverted().transform(bbox2)
            dr = np.array(bbox2d[0]-bbox1d[0])

            #the rotation/stretch matrix
            rad = rads[il]
            rot_mat = np.array([
                [math.cos(rad), math.sin(rad)*aspect],
                [-math.sin(rad)/aspect, math.cos(rad)]
            ])

            ##computing the offset vector of the rotated character
            drp = np.dot(dr,rot_mat)

            #setting final position and rotation:
            t.set_position(np.array([x,y])+drp)
            t.set_rotation(degs[il])

            t.set_va('center')
            t.set_ha('center')

            #updating rel_pos to right edge of character
            rel_pos += w-used

max_stars_repo_name: cajohare/solax

max_issues_repo_name: cajohare/solax

max_forks_repo_name: cajohare/solax

reason: import numpy,from numpy,from scipy


text:
from collections import defaultdict
import sys
import os
import json

import numpy as np
# from gensim.models.fasttext import load_facebook_model
import fasttext
from string import punctuation
from ruwordnet.ruwordnet_reader import RuWordnet

class FasttextVectorizer:
    def __init__(self, model_path):
        # self.model = load_facebook_model(model_path)
        self.model = fasttext.load_model(model_path)
        self.vector_size = self.model.get_dimension()
        print('Model loaded')

    # -------------------------------------------------------------
    # vectorize ruwordnet
    # -------------------------------------------------------------

    def vectorize_ruwordnet(self, synsets, output_path):
        ids, vectors = self.__get_ruwordnet_vectors(synsets)
        self.save_as_w2v(ids, vectors, output_path)

    def __get_ruwordnet_vectors(self, synsets):
        ids = []
        vectors = np.zeros((len(synsets), self.vector_size))
        for i, (_id, texts) in enumerate(synsets.items()):
            ids.append(_id)
            # vectors[i, :] = self.__get_avg_vector(texts)
            vectors[i, :] = self.model.get_sentence_vector((" ".join(texts)).lower())
            if i in range(2):
                print("number",i, (" ".join(texts)).lower(), vectors[i, :])
        return ids, vectors

    def __get_avg_vector(self, texts):
        sum_vector = np.zeros(self.vector_size)
        for text in texts:
            words = [i.strip(punctuation) for i in text.split()]
            sum_vector += np.sum(self.__get_data_vectors(words), axis=0)/len(words)
        return sum_vector/len(texts)

    # -------------------------------------------------------------
    # vectorize data
    # -------------------------------------------------------------

    def vectorize_data(self, data, output_path):
        data_vectors = self.__get_data_vectors(data)
        self.save_as_w2v(data, data_vectors, output_path)

    def __get_data_vectors(self, data):
        vectors = np.zeros((len(data), self.vector_size))
        for i, word in enumerate(data):
            vectors[i, :] = self.model[word]
        return vectors

    # -------------------------------------------------------------
    # save
    # -------------------------------------------------------------

    @staticmethod
    def save_as_w2v(words: list, vectors: np.array, output_path: str):
        assert len(words) == len(vectors)
        with open(output_path, 'w', encoding='utf-8') as w:
            w.write(f"{vectors.shape[0]} {vectors.shape[1]}\n")
            for word, vector in zip(words, vectors):
                vector_line = " ".join(map(str, vector))
                w.write(f"{word.upper()} {vector_line}\n")

def process_data(input_file, output_file):
    with open(input_file, 'r', encoding='utf-8') as f:
        dataset = f.read().lower().split("\n")[:-1]
    ft_vec.vectorize_data(dataset, output_file)

if __name__ == '__main__':
    from helpers.utils import load_config

    config = load_config()
    ft_vec = FasttextVectorizer(config["vectorizer_path"])
    ruwordnet = RuWordnet(db_path=config["db_path"], ruwordnet_path=config["ruwordnet_path"], with_lemmas=False)
    noun_synsets = defaultdict(list)
    verb_synsets = defaultdict(list)
    for sense_id, synset_id, text in ruwordnet.get_all_senses():
        if synset_id.endswith("N"):
            noun_synsets[synset_id].append(text.lower())
        elif synset_id.endswith("V"):
            verb_synsets[synset_id].append(text.lower())
    exp_suffix = config["experiment_suffix"]
    os.makedirs(f"predictions/{exp_suffix}", exist_ok=True)
    os.makedirs(f"models/vectors/{exp_suffix}", exist_ok=True)
    synsets = {
        'nouns_synsets':noun_synsets,
        'verbs_synsets':verb_synsets,
              }
    for part in ["nouns", "verbs"]:
        vectors_dict = {
            f"ruwordnet_vectors_{part}_path" : f"../baselines/models/vectors/{exp_suffix}/ruwordnet_{part}_{exp_suffix}.txt",

            f"public_data_vectors_{part}_path" : f"../baselines/models/vectors/{exp_suffix}/{part}_public_{exp_suffix}.txt",
            f"private_data_vectors_{part}_path" : f"../baselines/models/vectors/{exp_suffix}/{part}_private_{exp_suffix}.txt",

            f"public_test_{part}_path": f"../data/public_test/{part}_public.tsv",
            f"private_test_{part}_path": f"../data/private_test/{part}_private.tsv",

            f"public_output_{part}_path": f"predictions/{exp_suffix}/predicted_public_{part}_{exp_suffix}.tsv",
            f"private_output_{part}_path": f"predictions/{exp_suffix}/predicted_private_{part}_{exp_suffix}.tsv",
        }

        ft_vec.vectorize_ruwordnet(synsets[f"{part}_synsets"], vectors_dict[f"ruwordnet_vectors_{part}_path"])
        process_data(vectors_dict[f"public_test_{part}_path"], vectors_dict[f"public_data_vectors_{part}_path"])
        process_data(vectors_dict[f"private_test_{part}_path"], vectors_dict[f"private_data_vectors_{part}_path"])

        config.update(vectors_dict) # Update: config with paths
        with open(sys.argv[1], "w", encoding="utf-8") as j:
            json.dump(obj=config, fp=j)

ext: py

max_stars_repo_name: alexey-pronkin/taxonomy-enrichment

max_issues_repo_name: alexey-pronkin/taxonomy-enrichment

max_forks_repo_name: alexey-pronkin/taxonomy-enrichment

reason: import numpy


text:
from dnachisel import (
    DnaOptimizationProblem,
    EnforceTranslation,
    EnforceGCContent,
    random_protein_sequence,
    reverse_translate,
    reverse_complement,
    translate,
    AvoidPattern,
    EnforceChanges,
    Location,
)
import pytest
import numpy

# Note: we are not providing a location for AvoidChanges: it applies globally

def test_EnforceTranslation():
    numpy.random.seed(1234)
    sequence = reverse_translate(random_protein_sequence(50, seed=123))
    problem = DnaOptimizationProblem(
        sequence=sequence,
        constraints=[AvoidPattern("AAA"), EnforceTranslation()],
    )
    assert not problem.all_constraints_pass()
    problem.resolve_constraints()
    assert problem.all_constraints_pass()

def test_EnforceTranslationReversed():
    numpy.random.seed(1234)
    sequence = reverse_translate(random_protein_sequence(50, seed=123))
    rev_sequence = reverse_complement(sequence)
    problem = DnaOptimizationProblem(
        sequence=rev_sequence,
        constraints=[
            AvoidPattern("AGC"),
            EnforceTranslation(location=(0, len(sequence), -1)),
        ],
        logger=None,
    )
    assert not problem.all_constraints_pass()
    problem.resolve_constraints()
    assert problem.all_constraints_pass()

def test_EnforceTranslation_error_location_not_3x():
    """Providing a location that is not multiple of 3 raises an error"""
    numpy.random.seed(1234)
    sequence = reverse_translate(random_protein_sequence(50, seed=123))
    with pytest.raises(ValueError) as err:
        _ = DnaOptimizationProblem(
            sequence=sequence,
            constraints=[EnforceTranslation(location=(0, 16))],
            logger=None,
        )
    assert "Location 0-16(+) has length 16" in str(err.value)

def test_EnforceTranslation_error_location_smaller_than_translation():
    """Providing a location that is not multiple of 3 raises an error"""
    numpy.random.seed(1234)
    sequence = reverse_translate(random_protein_sequence(15, seed=123))
    with pytest.raises(ValueError) as err:
        _ = DnaOptimizationProblem(
            sequence=sequence,
            constraints=[
                EnforceTranslation(
                    translation=random_protein_sequence(30, seed=111)
                )
            ],
            logger=None,
        )
    assert str(err.value).startswith("Window size")

def test_EnforceTranslation_bacterial_valine():
    table_name = "Bacterial"
    protein = "LLTMMVTTTTVMVL"
    protein_sequence = reverse_translate(protein, table=table_name)

    for first_codon_before, first_codon_after in [
        ("ATG", "ATG"),  # methionine stays the only methionine codon
        ("GTG", "GTG"),  # valine-start-codon stays the only valine-start-codon
    ]:
        sequence = first_codon_before + protein_sequence
        cds_constraint = EnforceTranslation(
            genetic_table="Bacterial", start_codon="keep"
        )
        problem = DnaOptimizationProblem(
            sequence=sequence,
            constraints=[cds_constraint],
            objectives=[EnforceChanges()],
            logger=None,
        )
        assert problem.constraints[0].translation == "MLLTMMVTTTTVMVL"
        problem.optimize()
        protein_after = translate(
            problem.sequence, table_name, assume_start_codon=True
        )
        assert protein_after == "M" + protein
        assert problem.sequence[:3] == first_codon_after

def test_EnforceTranslation_bacterial_valine_antisense():
    table_name = "Bacterial"
    protein = "LLTMMVTTTTVMVL"
    protein_sequence = reverse_translate(protein, table=table_name)
    for first_codon_before, first_codon_after in [
        ("ATG", "ATG"),  # methionine stays the only methionine codon
        ("GTG", "GTG"),  # valine-start-codon stays the only valine-start-codon
    ]:
        sequence = first_codon_before + protein_sequence
        cds_constraint = EnforceTranslation(
            genetic_table="Bacterial",
            start_codon="keep",
            location=Location(0, len(sequence), -1),
        )
        problem = DnaOptimizationProblem(
            sequence=reverse_complement(sequence),
            constraints=[cds_constraint],
            objectives=[EnforceChanges()],
            logger=None,
        )
        assert problem.constraints[0].translation == "MLLTMMVTTTTVMVL"
        problem.optimize()
        problem_sequence_rv = reverse_complement(problem.sequence)
        protein_after = translate(
            problem_sequence_rv, table_name, assume_start_codon=True
        )
        assert protein_after == "M" + protein
        assert problem_sequence_rv[:3] == first_codon_after

max_stars_repo_name: Edinburgh-Genome-Foundry/DnaChisel

max_issues_repo_name: Edinburgh-Genome-Foundry/DnaChisel

max_forks_repo_name: Edinburgh-Genome-Foundry/DnaChisel


text:
# -*- coding: utf-8 -*-
'''
This module defines :class:`EventArray`, an array of events. Introduced for
performance reasons.

:class:`EventArray` derives from :class:`BaseNeo`, from
:module:`neo.core.baseneo`.
'''

# needed for python 3 compatibility
from __future__ import absolute_import, division, print_function

import sys

import numpy as np
import quantities as pq

from neo.core.baseneo import BaseNeo, merge_annotations

PY_VER = sys.version_info[0]

class EventArray(BaseNeo):
    '''
    Array of events. Introduced for performance reasons.

    An :class:`EventArray` is prefered to a list of :class:`Event` objects.

    *Usage*::

        >>> from neo.core import EventArray
        >>> from quantities import s
        >>> import numpy as np
        >>>
        >>> evtarr = EventArray(np.arange(0, 30, 10)*s,
        ...                     labels=np.array(['trig0', 'trig1', 'trig2'],
        ...                                     dtype='S'))
        >>>
        >>> evtarr.times
        array([  0.,  10.,  20.]) * s
        >>> evtarr.labels
        array(['trig0', 'trig1', 'trig2'],
              dtype='|S5')

    *Required attributes/properties*:
        :times: (quantity array 1D) The time of the events.
        :labels: (numpy.array 1D dtype='S') Names or labels for the events.

    *Recommended attributes/properties*:
        :name: (str) A label for the dataset.
        :description: (str) Text description.
        :file_origin: (str) Filesystem path or URL of the original data file.

    Note: Any other additional arguments are assumed to be user-specific
            metadata and stored in :attr:`annotations`.

    '''

    _single_parent_objects = ('Segment',)
    _necessary_attrs = (('times', pq.Quantity, 1),
                       ('labels', np.ndarray, 1, np.dtype('S')))

    def __init__(self, times=None, labels=None, name=None, description=None,
                 file_origin=None, **annotations):
        '''
        Initialize a new :class:`EventArray` instance.
        '''
        BaseNeo.__init__(self, name=name, file_origin=file_origin,
                         description=description, **annotations)
        if times is None:
            times = np.array([]) * pq.s
        if labels is None:
            labels = np.array([], dtype='S')

        self.times = times
        self.labels = labels

        self.segment = None

    def __repr__(self):
        '''
        Returns a string representing the :class:`EventArray`.
        '''
        # need to convert labels to unicode for python 3 or repr is messed up
        if PY_VER == 3:
            labels = self.labels.astype('U')
        else:
            labels = self.labels
        objs = ['%s@%s' % (label, time) for label, time in zip(labels,
                                                               self.times)]
        return '<EventArray: %s>' % ', '.join(objs)

    def merge(self, other):
        '''
        Merge the another :class:`EventArray` into this one.

        The :class:`EventArray` objects are concatenated horizontally
        (column-wise), :func:`np.hstack`).

        If the attributes of the two :class:`EventArray` are not
        compatible, and Exception is raised.
        '''
        othertimes = other.times.rescale(self.times.units)
        times = np.hstack([self.times, othertimes]) * self.times.units
        labels = np.hstack([self.labels, other.labels])
        kwargs = {}
        for name in ("name", "description", "file_origin"):
            attr_self = getattr(self, name)
            attr_other = getattr(other, name)
            if attr_self == attr_other:
                kwargs[name] = attr_self
            else:
                kwargs[name] = "merge(%s, %s)" % (attr_self, attr_other)

        merged_annotations = merge_annotations(self.annotations,
                                               other.annotations)
        kwargs.update(merged_annotations)
        return EventArray(times=times, labels=labels, **kwargs)

max_stars_repo_name: tclose/python-neo

max_issues_repo_name: tclose/python-neo

max_forks_repo_name: tclose/python-neo


text:
import numpy as np
import datetime as dt
from datetime import datetime

from pyspark.sql.functions import from_json, col
from pyspark.sql.types import *
from pyspark.ml.feature import VectorAssembler
from pyspark.ml.feature import OneHotEncoderEstimator, StringIndexer
from pyspark.ml import Pipeline
from pyspark.ml.regression import LinearRegression
from pyspark.ml.regression import GBTRegressor

def extract(row, colname):
		return tuple(map(lambda x: row[x], row.__fields__)) + tuple(row[colname + "_ohe"].toArray().tolist())

def run(spark, message_dataframe, timetable, stations, ris, train_dict):

	jsonSchema = spark.read.json(message_dataframe.rdd.map(lambda r: r.payload)).schema
	trainEvent = message_dataframe.withColumn('json', from_json(col('payload'), jsonSchema)).select("station_id", "ris_status_id", "timetable_id", "json")

	trainEvent = trainEvent.select(trainEvent.station_id,
			trainEvent.ris_status_id,
			trainEvent.timetable_id,
			trainEvent.json.info.event_datetime.alias("datetime"),
			trainEvent.json.info.train_id.alias("train_id"))

	ris_status_id = trainEvent.collect()[0]["ris_status_id"]
	station_id = trainEvent.collect()[0]["station_id"]
	timetable_id = trainEvent.collect()[0]["timetable_id"]
	datetime_str =  trainEvent.collect()[0]["datetime"]

	ris_status_name = spark.sql("select name from ris where id == " + str(ris_status_id)).collect()[0]["name"]

	trainEvent.registerTempTable("trainEvent")
	trainEvent = spark.sql("select station_id, hour(datetime) as hour, dayofmonth(datetime) AS day, month(datetime) as month from trainEvent")

	feature_columns = ["station_id", "hour", "day", "month"]
	keep_columns = [c for c in train_dict["keepStationColumns"] if c != "dur_stop"]

	medians = dict(list(map(lambda x: (float(x.station_id), x.median), train_dict["medianTimePerStation"].collect())))
	stops_time = 0.0
	if ris_status_name == "Arrival":
		transTrainEvent = train_dict["pipelineIndexerModel"].transform(trainEvent)
		encodedTrainEvent = train_dict["oheStationModel"].transform(transTrainEvent)
		for colname in feature_columns:
			colIdx =  train_dict["stationNameDict"][colname]
			newCols = list(map(lambda x: str(x[1]).split(".")[0], colIdx))
			actualCol = encodedTrainEvent.columns
			allColNames = actualCol + newCols
			encodedTrainEvent = encodedTrainEvent.rdd.map(lambda x : extract(x, colname)).toDF(allColNames)
		encodedTrainEvent = encodedTrainEvent[keep_columns]
		encodedTrainEvent = train_dict["stationTrainAssembler"].transform(encodedTrainEvent)
		encodedTrainEvent = encodedTrainEvent.select("features")
		prediction = train_dict["stationTrainModel"].transform(encodedTrainEvent)
		stops_time += prediction.collect()[0].prediction

	station_rows = spark.sql("select station_id, station_order from timetable where timetable_id == " + timetable_id + " and station_order > (select station_order from timetable where station_id == '" + str(station_id) + "') and station_order < (select max(station_order) from timetable where timetable_id == " + timetable_id + ") order by station_order").collect()
	station_rows = list(map(lambda x : x.station_id, station_rows))

	for station in station_rows:
		if station in medians.keys():
			stops_time += medians[station]
		else:
			stops_time += spark.sql("select departure_second_of_the_day - arrival_second_of_the_day as stop from timetable where timetable_id == " + str(timetable_id) + " and station_id ==  " + str(int(station))).collect()[0]["stop"]

	next_route = spark.sql("select start_lon, start_lat, end_lon, end_lat, hour, day, month from (select longitude as start_lon, latitude as start_lat, hour(e.datetime) as hour, dayofmonth(e.datetime) AS day, month(e.datetime) as month, t.station_order as ord, e.timetable_id as tid from trainEvent as e, stations as s, timetable as t where e.station_id == s.id and t.timetable_id == e.timetable_id and t.station_id == e.station_id) as s, (select longitude as end_lon, latitude as end_lat, station_order as ord, timetable_id as tid from timetable as t, stations as s where t.station_id == s.id) as e where s.tid == e.tid and s.ord + 1 == e.ord")

	feature_columns = ["hour", "day", "month"]
	keep_columns = [c for c in train_dict["keepRouteColumns"] if c != "dur_route"]

	routes_time = 0.0
	next_route = train_dict["pipelineIndexerModel"].transform(next_route)
	encoded_train_route_data = train_dict["oheRouteModel"].transform(next_route)
	for colname in feature_columns:
        	colIdx = train_dict["routeNameDict"][colname]
	        newCols = list(map(lambda x: str(x[1]).split(".")[0], colIdx))
	        actualCol = encoded_train_route_data.columns
	        allColNames = actualCol + newCols
	        encoded_train_route_data = encoded_train_route_data.rdd.map(lambda x : extract(x, colname)).toDF(allColNames)

	encoded_train_route_data = encoded_train_route_data[keep_columns]
	encoded_train_route_data = train_dict["routeTrainAssembler"].transform(encoded_train_route_data)
	encoded_train_route_data = encoded_train_route_data.select("features")
	prediction = train_dict["routeTrainModel"].transform(encoded_train_route_data)
	routes_time += prediction.collect()[0].prediction

	rest_routes = spark.sql("select start_lon, start_lat, end_lon, end_lat from (select s.longitude as start_lon, s.latitude as start_lat, t.station_order as ord from timetable as t, stations as s, trainEvent as e where e.timetable_id == t.timetable_id and t.station_id == s.id and t.station_order > (select station_order from trainEvent as e, timetable as t where e.station_id == t.station_id and e.timetable_id == t.timetable_id)) as s, (select s.longitude as end_lon, s.latitude as end_lat, t.station_order as ord from timetable as t, stations as s, trainEvent as e where e.timetable_id == t.timetable_id and t.station_id == s.id and t.station_order > (select station_order from trainEvent as e, timetable as t where e.station_id == t.station_id and e.timetable_id == t.timetable_id)) as e where s.ord + 1 = e.ord")

	rest_routes = train_dict["routePartialTrainAssembler"].transform(rest_routes)
	rest_routes = rest_routes.select('features')
	predictions = train_dict["routePartialTrainModel"].transform(rest_routes)
	routes_time += predictions.groupBy().sum().collect()[0][0]

	total_time = stops_time + routes_time
	eta = str((datetime.strptime(datetime_str, '%Y-%m-%d %H:%M:%S') + dt.timedelta(0, seconds = total_time)))
	print(eta)
	return(eta)

max_stars_repo_name: selisproject/bda

max_issues_repo_name: selisproject/bda

max_forks_repo_name: selisproject/bda


text:
import cv2
import random
import os

import numpy as np
import matplotlib.pyplot as plt

import torch
from torch.utils.data import Dataset
from torchvision import transforms

class ShanghaiTechPartA(Dataset):
    def __init__(self, root, shuffle=False, transform=None, downsample=1):
        self.root = root
        self.shuffle = shuffle
        self.transform = transform
        self.downsample = downsample

        self.image_names = [filename for filename in os.listdir(self.root)]
        self.n_samples = len(self.image_names)

        if self.shuffle:
            random.shuffle(self.image_names)

    def __len__(self):
        return self.n_samples

    def __getitem__(self, index):
        assert index <= len(self), 'index range error'
        img_name = self.image_names[index]
        # Read image and normalize its pixels to [0,1]
        img = plt.imread(os.path.join(self.root,img_name)) / 255
        # Expand grayscale image to three channel.
        if len(img.shape) == 2:
            img = img[:,:,np.newaxis]
            img = np.concatenate((img,img,img),2)

        # Read ground truth density-map
        density_map = np.load(os.path.join(self.root.replace('images','density_maps'),img_name.replace('.jpg','.npy')))

        # Downsample image and density-map to match model's input
        if self.downsample >1:
            rows = int(img.shape[0] // self.downsample)
            cols = int(img.shape[1] // self.downsample)
            img = cv2.resize(img,(cols*self.downsample, rows*self.downsample))
            img = img.transpose((2,0,1)) # convert to order (channel,rows,cols)
            density_map = cv2.resize(density_map, (cols,rows))
            density_map = density_map[np.newaxis,:,:] * self.downsample * self.downsample
            # transform image and density_map to tensors
            img_tensor = torch.tensor(img, dtype=torch.float)
            density_map_tensor = torch.tensor(density_map, dtype=torch.float)
            # Apply any other transformation
            if self.transform is not None:
                img_tensor = self.transform(img_tensor)

        return img_tensor, density_map_tensor

# Test code
if __name__== "__main__":
    root = '/Users/pantelis/Downloads/archive/ShanghaiTech/part_A/train_data/images'
    dataset = ShanghaiTechPartA(root,
                                transform=transforms.Normalize(
                                            mean=[0.485, 0.456, 0.406],
                                            std=[0.229, 0.224, 0.225]),
                                downsample=8)
    index = random.randint(0, len(dataset))
    img, dmap = dataset[index]
    print(index, img.shape, dmap.shape)

max_stars_repo_name: PanPapag/Context-Aware-Crowd-Counting

max_issues_repo_name: PanPapag/Context-Aware-Crowd-Counting

max_forks_repo_name: PanPapag/Context-Aware-Crowd-Counting


text:
#### Created on 2016/11/08 by Taufiq ###

import matplotlib.pyplot as plt
import matplotlib.animation as animation
import matplotlib.patches as mpatches
import pandas as pd
import numpy as np
from datetime import datetime

col = ['red', 'green', 'blue', 'yellow', 'black', 'green', 'red', 'red', 'green', 'blue', 'yellow', 'green',
       'blue', 'yellow', 'green', 'blue', 'blue', 'red', 'blue', 'yellow', 'blue', 'blue', 'yellow', 'red', 'yellow',
       'blue', 'blue', 'blue', 'yellow', 'blue', 'green', 'yellow', 'green', 'green', 'blue', 'yellow', 'yellow',
       'blue',
       'yellow', 'blue', 'blue', 'blue', 'green', 'green', 'blue', 'blue', 'green', 'blue', 'green', 'yellow', 'blue',
       'blue', 'yellow', 'yellow', 'red', 'green', 'green', 'red', 'red', 'red', 'red', 'green', 'red', 'green',
       'yellow',
       'red', 'red', 'blue', 'red', 'red', 'red', 'red', 'blue', 'blue', 'blue', 'blue', 'blue', 'red', 'blue', 'blue',
       'blue', 'yellow', 'red', 'green', 'blue', 'blue', 'red', 'blue', 'red', 'green', 'black', 'yellow', 'blue',
       'blue',
       'green', 'red', 'red', 'yellow', 'yellow', 'yellow', 'red', 'green', 'green', 'yellow', 'blue', 'green', 'blue',
       'blue', 'red', 'blue', 'green', 'blue', 'red', 'green', 'green', 'blue', 'blue', 'green', 'red', 'blue', 'blue',
       'green', 'green', 'red', 'red', 'blue', 'red', 'blue', 'yellow', 'blue', 'green', 'blue', 'green', 'yellow',
       'yellow', 'yellow', 'red', 'red', 'red', 'blue', 'blue']

class GraphData:
    def __init__(self, file, title, xlabel, ylabel, grid, x, y, threshold):
        print("Graph initializing....")
        self.file = file
        self.title = title
        self.xlabel = xlabel
        self.ylabel = ylabel
        self.grid = grid
        self.x = x
        self.y = y
        self.threshold = threshold

        print("GD constructor initialized with:", file, title, xlabel, ylabel, grid, x, y, threshold)

    def CSV_Handler(self):
        print("Handling CSV file....")
        # delimiter = ',' and '\t' and ':'
        self.df = pd.read_csv(self.file, sep=',', header=0, engine='python')
        print(self.df)
        self.graph()

    def file_type(self):
        print("file_type processing....")
        if self.file.lower().endswith(('.txt', '.csv')):
            self.CSV_Handler()
        elif self.file.lower().endswith(('.xlrs', '.xlr', '.xlsx')):
            # Excel_Handler.__init__(self.file)
            print("Support for Excel coming Soon!")
        else:
            print("Unsupported Filetype!")
        print("[Process Finished]")

    def details(self):
        plt.xlabel(self.xlabel)
        plt.ylabel(self.ylabel)
        plt.title(self.title)
        #plt.title("Histogram of IQ: " r'$\mu = 100, \sigma$ =15')
        plt.minorticks_on()
        plt.grid(b=self.grid, which='minor', color='w', linestyle='-')
        plt.grid(b=self.grid, which='major', color='#737e8e', linestyle='-')

        plt.rc('axes', edgecolor='r')

    def graph(self):
        plt.style.use(['fivethirtyeight', 'ggplot'])
        x = self.df[self.x]
        y = self.df[self.y]

        np_pop = np.array(y)
        np_pop = np_pop * 8

        self.details()
        Color = 'c'
        plt.plot(x, y, marker='o', color=Color)
        plt.grid(self.grid)
        plt.legend(handles=[mpatches.Patch(color=Color, label=self.title)])
        if self.threshold < 1:
            print("threshold is less than 1 ")
        else:
            plt.axhline(y=self.threshold, xmin=0, xmax=1, hold=None, color='m')
        plt.show()
        plt.savefig(self.title + " " + datetime.now().strftime('%Y-%m-%d %H:%M') + '.png')
        self.details()
        plt.scatter(x, y, s=np_pop, c=col, alpha=0.8)
        plt.grid(self.grid)
        if self.threshold < 1:
            print("threshold is less than 1 ")
        else:
            plt.axhline(y=self.threshold, xmin=0, xmax=1, hold=None, color='m')
        plt.show()
        plt.savefig(self.title + " " + datetime.now().strftime('%Y-%m-%d %H:%M') + " " + '.png')
        self.details()
        plt.hist(y, bins=12)
        plt.grid(self.grid)
        plt.show()
        plt.savefig(self.title + " " + datetime.now().strftime('%Y-%m-%d %H:%M') + '.png')

class LiveGraph:
    def __init__(self, file):
        self.file = file
        self.fig = plt.figure()
        self.ax1 = self.fig.add_subplot(1, 1, 1)

    def animate(self, i):
        plt.style.use('fivethirtyeight')
        graph_data = open(self.file, 'r').read()
        lines = graph_data.split('\n')
        xs = []
        ys = []
        for line in lines:
            if len(line) > 1:
                x, y = line.split(',')
                xs.append(x)
                ys.append(y)
        self.ax1.clear()
        self.ax1.plot(xs, ys)

    def showupdate(self):
        ani = animation.FuncAnimation(self.fig, self.animate, interval=1000)
        plt.show()


text:
import networkx as nx
import numpy as np
import pandas as pd
import pytest

@pytest.fixture
def edgelist1():
    return [[1, 2, 3], [4], [5, 6], [6, 7, 8]]

@pytest.fixture
def edgelist2():
    return [[1, 2], [3, 4], [4, 5, 6]]

@pytest.fixture
def edgelist3():
    return [[1, 2, 3], [3, 4], [4, 5, 6]]

@pytest.fixture
def edgelist4():
    return [[1, 2, 3], [2, 3, 4, 5], [3, 4, 5]]

@pytest.fixture
def edgelist5():
    return [[0, 1, 2, 3], [4], [5, 6], [6, 7, 8]]

@pytest.fixture
def edgelist6():
    return [[0, 1, 2], [1, 2, 3], [2, 3, 4]]

@pytest.fixture
def edgelist7():
    return [[0, 1, 2], [1, 2, 3], [2, 3, 4], [4]]

@pytest.fixture
def dict5():
    return {0: [0, 1, 2, 3], 1: [4], 2: [5, 6], 3: [6, 7, 8]}

@pytest.fixture
def incidence5():
    I = np.zeros((9, 4))
    I[0, 0] = I[1, 0] = I[2, 0] = I[3, 0] = I[4, 1] = 1
    I[5, 2] = I[6, 2] = I[6, 3] = I[7, 3] = I[8, 3] = 1
    return I

@pytest.fixture
def dataframe5():
    data = [
        [0, 0],
        [1, 0],
        [2, 0],
        [3, 0],
        [4, 1],
        [5, 2],
        [6, 2],
        [6, 3],
        [7, 3],
        [8, 3],
    ]
    return pd.DataFrame(data)

@pytest.fixture
def bipartite_graph1():
    G = nx.Graph()
    G.add_nodes_from([1, 2, 3, 4], bipartite=0)
    G.add_nodes_from(["a", "b", "c"], bipartite=1)
    G.add_edges_from([(1, "a"), (1, "b"), (2, "b"), (2, "c"), (3, "c"), (4, "a")])
    return G

@pytest.fixture
def bipartite_graph2():
    G = nx.Graph()
    G.add_nodes_from([1, 2, 3, 4], bipartite=0)
    G.add_nodes_from(["a", "b", "c"], bipartite=2)
    G.add_edges_from([(1, "a"), (1, "b"), (2, "b"), (2, "c"), (3, "c"), (4, "a")])
    return G

@pytest.fixture
def bipartite_graph3():
    G = nx.Graph()
    G.add_nodes_from([1, 2, 3, 4])
    G.add_nodes_from(["a", "b", "c"])
    G.add_edges_from([(1, "a"), (1, "b"), (2, "b"), (2, "c"), (3, "c"), (4, "a")])
    return G

@pytest.fixture
def bipartite_graph4():
    # this is to test when the bipartite condition is not true
    G = nx.Graph()
    G.add_nodes_from([1, 2, 3, 4], bipartite=0)
    G.add_nodes_from(["a", "b", "c"], bipartite=1)
    G.add_edges_from([(1, "a"), ("a", "b"), (2, "b"), (2, "c"), (3, "c"), (4, "a")])
    return G

@pytest.fixture
def attr1():
    return {"color": "red", "name": "horse"}

@pytest.fixture
def attr2():
    return {"color": "blue", "name": "pony"}

@pytest.fixture
def attr3():
    return {"color": "yellow", "name": "zebra"}

max_stars_repo_name: ComplexGroupInteractions/xgi

max_issues_repo_name: ComplexGroupInteractions/xgi

max_forks_repo_name: ComplexGroupInteractions/xgi

reason: import numpy,import networkx


text:
import numpy as np
import matplotlib.pyplot as plt
import librosa
import librosa.display
import os
from collections import deque
import pandas as pd

WINDOW_SIZE = 50
WINDOW_DIST = 25

fail = []

flag = False

def filter(song, sr):
  D = librosa.stft(song)
  D_harmonic, D_percussive = librosa.decompose.hpss(D)
  S = D_harmonic
  s, phase = librosa.magphase(S)
  S_filter = librosa.decompose.nn_filter(s,
                                       aggregate=np.median,
                                       metric='cosine',
                                       width=int(librosa.time_to_frames(2, sr=sr)))
  S_filter = np.minimum(s, S_filter)
  i_mask = librosa.util.softmask(S_filter,
                                2* (s - S_filter),
                                power=2)
  pure_instrumental = s * i_mask
  return pure_instrumental

def window(seq, n=WINDOW_SIZE):
    it = iter(seq)
    win = deque((next(it, None) for _ in range(n)), maxlen=n)
    yield win
    append = win.append
    for e in it:
        append(e)
        yield win

num_musicas = len(os.listdir('DataBase/Audio'))

for i in range(0, num_musicas):

    music_name = os.listdir('DataBase/Audio')[i]
    if music_name[num_musicas-5:] == ".wav":
        music_name = music_name[:num_musicas - 5]
    elif music_name[num_musicas-6:] == ".wav":
        music_name = music_name[:num_musicas - 6]

    if music_name == 'I Want A New Drug.wav' or flag == True:
        flag = True
        data, sr = librosa.load(f'DataBase/Audio/{music_name}', sr = None)

        D = filter(data, sr)

        n = D.shape[0]
        yout = librosa.fft_frequencies(sr=sr ,n_fft=1+(2 * (n - 1)) )

        m = D.shape[1]
        hop_length=512

        xout = librosa.frames_to_time(np.arange(m+1), sr=sr, hop_length=hop_length)

        print('\n-----------------------------------------------------------------------------------')
        print(f'Musica: {music_name}')
        print(f'Duracao: {xout[-1]}')
        print(f'Dimensões(f/t): {D.shape}')
        print('-----------------------------------------------------------------------------------\n')

        window_decibel = np.zeros((n,WINDOW_SIZE))

        x1 = WINDOW_DIST
        a = 0

        for w in window(xout, WINDOW_SIZE):

            if x1 >= WINDOW_DIST:
                try:

                    for k in range(0, n):
                        for l in range(0, len(w)):
                            window_decibel[k][l] = D[k][l + a]

                    plt.figure(figsize=(20,10))
                    librosa.display.specshow(librosa.amplitude_to_db(window_decibel), x_axis='time', y_axis='log', sr= sr, cmap = 'jet')
                    plt.xlabel("Time [s]")
                    plt.ylabel("Frequency [Hz]")
                    plt.colorbar(format = "%+2.0f db")
                    plt.savefig(f'DataBase/new_windows/{music_name}{a/WINDOW_DIST}.png', dpi = 70)
                    plt.close()

                    a += WINDOW_DIST
                    x1 = 0

                except:
                    fail.append(f'{music_name}{a/WINDOW_DIST}.png')
                    print(f'Failed to save: {music_name}{a/WINDOW_DIST}.png')

            x1 += 1

print(fail)


text:
#!/usr/bin/env python

import numpy as np
from pycrazyswarm import *
import rospy
from crazyflie_driver.msg import GenericLogData
import uav_trajectory
import argparse

file = None
HEIGHT = 1.0 #0.75

# "ctrltarget.x", "ctrltarget.y", "ctrltarget.z", "stateEstimate.x", "stateEstimate.y", "stateEstimate.z"
def logData(data):
    global file

    print(data.values)
    file.write(",".join([str(i) for i in data.values]) + "\n")

if __name__ == "__main__":
    parser = argparse.ArgumentParser()
    parser.add_argument("speed", type=str, help="speed of trajectory")
    parser.add_argument("output", type=str, help="output file name")
    args = parser.parse_args()

    swarm = Crazyswarm(parse_args = False)
    timeHelper = swarm.timeHelper
    allcfs = swarm.allcfs

    file = open(args.output, "w")

    rospy.Subscriber("/cf{}/log1/".format(allcfs.crazyflies[0].id), GenericLogData, logData)

    # swarm = Crazyswarm()
    # timeHelper = swarm.timeHelper
    # allcfs = swarm.allcfs

    traj1 = uav_trajectory.Trajectory()
    traj1.loadcsv("figure8.csv", swapxy=True)

    TIMESCALE = float(args.speed) #0.75
    for cf in allcfs.crazyflies:
        cf.uploadTrajectory(0, 0, traj1)

    allcfs.takeoff(targetHeight=HEIGHT, duration=3.0)
    timeHelper.sleep(2.5)
    for cf in allcfs.crazyflies:
        # pos = np.array(cf.initialPosition) + np.array([0, 0, 1.0])
        pos = np.array([0, 0, HEIGHT])
        cf.goTo(pos, 0, 2.0)
    timeHelper.sleep(2.5)

    allcfs.startTrajectory(0, timescale=TIMESCALE)
    timeHelper.sleep(traj1.duration * TIMESCALE + 2.0)
    # allcfs.startTrajectory(0, timescale=TIMESCALE, reverse=True)
    # timeHelper.sleep(traj1.duration * TIMESCALE + 2.0)

    allcfs.land(targetHeight=0.06, duration=2.0)
    timeHelper.sleep(2.0)

max_stars_repo_name: SumeetBatra/crazyswarm-nn

max_issues_repo_name: TaoChenOSU/crazyswarm-nn

max_forks_repo_name: TaoChenOSU/crazyswarm-nn


text:
import petsc4py
import sys

petsc4py.init(sys.argv)
from petsc4py import PETSc
import numpy as np
from dolfin import tic, toc
import HiptmairSetup
import PETScIO as IO
import scipy.sparse as sp
# import matplotlib.pylab as plt
import MatrixOperations as MO
import HiptmairSetup

class BaseMyPC(object):
    def setup(self, pc):
        pass
    def reset(self, pc):
        pass
    def apply(self, pc, x, y):
        raise NotImplementedError
    def applyT(self, pc, x, y):
        self.apply(pc, x, y)
    def applyS(self, pc, x, y):
        self.apply(pc, x, y)
    def applySL(self, pc, x, y):
        self.applyS(pc, x, y)
    def applySR(self, pc, x, y):
        self.applyS(pc, x, y)
    def applyRich(self, pc, x, y, w, tols):
        self.apply(pc, x, y)

class Matrix(object):
    def __init__(self):
        pass
    def create(self, mat):
        pass
    def destroy(self, mat):
        pass

class ShiftedCurlCurl(BaseMyPC):

    def __init__(self, Z, M, T):
        self.T = T
        self.Z = Z
        self.M = M

    def create(self, pc):
        print "Create"

    def setUp(self, pc):
        kspT = PETSc.KSP()
        kspT.create(comm=PETSc.COMM_WORLD)
        pcT = kspT.getPC()
        kspT.setType('preonly')
        pcT.setType('hypre')
        kspT.setOperators(self.T, self.T)
        options = PETSc.Options()
        # options['pc_hypre_boomeramg_cycle_type']  = "W"
        #     #
        options["pc_hypre_type"] = "boomeramg"
        options["pc_hypre_boomeramg_nodal_coarsen"] = 6
        options["pc_hypre_boomeramg_vec_interp_variant"] = 2
        kspT.setFromOptions()
        self.kspT = kspT

        kspM = PETSc.KSP()
        kspM.create(comm=PETSc.COMM_WORLD)
        pcM = kspM.getPC()
        kspM.setType('preonly')
        pcM.setType('hypre')
        kspM.setOperators(self.M, self.M)
        self.kspM = kspM

        print "setup"
    def apply(self, pc, x, y):
        x1 = x.duplicate()
        x2 = x.duplicate()
        x3 = x.duplicate()
        self.kspT.solve(x, x1)
        self.Z.mult(x1, x2)
        self.kspM.solve(x2, x3)
        y.array = x3.array

max_stars_repo_name: wathen/PhD

max_issues_repo_name: wathen/PhD

max_forks_repo_name: wathen/PhD

reason: import numpy,import scipy


text:
import numpy as np

from poliastro.core.util import circular_velocity

from .._jit import jit

@jit
def extra_quantities(k, a_0, a_f, inc_0, inc_f, f):
    """Extra quantities given by the Edelbaum (a, i) model.
    """
    V_0, beta_0_, _ = compute_parameters(k, a_0, a_f, inc_0, inc_f)
    delta_V_ = delta_V(V_0, beta_0_, inc_0, inc_f)
    t_f_ = delta_V_ / f

    return delta_V_, t_f_

@jit
def beta(t, *, V_0, f, beta_0):
    """Compute yaw angle (β) as a function of time and the problem parameters.
    """
    return np.arctan2(V_0 * np.sin(beta_0), V_0 * np.cos(beta_0) - f * t)

@jit
def beta_0(V_0, V_f, inc_0, inc_f):
    """Compute initial yaw angle (β) as a function of the problem parameters.
    """
    delta_i_f = abs(inc_f - inc_0)
    return np.arctan2(
        np.sin(np.pi / 2 * delta_i_f), V_0 / V_f - np.cos(np.pi / 2 * delta_i_f)
    )

@jit
def compute_parameters(k, a_0, a_f, inc_0, inc_f):
    """Compute parameters of the model.
    """
    delta_inc = abs(inc_f - inc_0)
    V_0 = circular_velocity(k, a_0)
    V_f = circular_velocity(k, a_f)
    beta_0_ = beta_0(V_0, V_f, inc_0, inc_f)

    return V_0, beta_0_, delta_inc

@jit
def delta_V(V_0, beta_0, inc_0, inc_f):
    """Compute required increment of velocity.
    """
    delta_i_f = abs(inc_f - inc_0)
    return V_0 * np.cos(beta_0) - V_0 * np.sin(beta_0) / np.tan(
        np.pi / 2 * delta_i_f + beta_0
    )

max_stars_repo_name: sundeshgupta/poliastro

max_issues_repo_name: sundeshgupta/poliastro

max_forks_repo_name: sundeshgupta/poliastro


text:
import scipy.sparse as sp
import numpy as np
import mkl

def estimate_memory_matrix_product(matrix, k):
    row = matrix.shape[0]
    d = row * k * 32 / 8 / 1024 / 1024
    maxk = 2147483647/(row)
    dmax = row * maxk * 32 / 8 / 1024 / 1024
    print('array data (float32) = %d MB, array row (int32) = %d MB, array col (int32)= %d MB'%(d,d,d))
    print('total = %d MB \t (%.2f GB)'%(3*d,3*d/1024))

def max_k_matrix_product(matrix):
    row = matrix.shape[0]
    maxk = 2147483647/(row)
    dmax = row * maxk * 32 / 8 / 1024 / 1024
    print('max value K due to integer range limit = %d (size in memory %.2f GB)'%(maxk,3*dmax/1024))

def get_mkl_threads():
    print('mkl max threads: %d'%mkl.get_max_threads())

def set_mkl_threads(num_threads=1):
    try:
        mkl.set_num_threads(num_threads)
        assert(mkl.get_max_threads()==num_threads)
    except:
        print('!!! ERROR: setting mkl num threads !!!')

def show_config_numpy(num_threads):
    np.__config__.show()

max_stars_repo_name: edervishaj/spotify-recsys-challenge

max_issues_repo_name: kiminh/spotify-recsys-challenge

max_forks_repo_name: kiminh/spotify-recsys-challenge


text:
"""
File contains functions used to train Linear Support Vector Machine.

Author: Dominik Macko
"""

from typing import Dict, Any, Union, Callable, Tuple, Optional

import numpy as np
import pandas as pd
from sklearn.svm import LinearSVC
import joblib

from androidcrypto.analysis.training.generic import train_model, train_gridsearchcv_model
from androidcrypto.analysis.training.config import TrainingTaskConfig
from androidcrypto.analysis.training.utils import is_multiclass

def linear_svm_best_params_surroundings(best_params: Dict[str, Any]) -> Dict[str, Any]:
    """Get best parameters surroundings for random forest."""

    C = best_params["C"]
    return {
        "C": [1/2 * C, C, 3/2 * C]
    }

def train_linear_svm(train_X: np.array,
                     train_y: np.array,
                     scoring: Union[str, Callable[[Any, np.array, np.array], int]]="f1_macro",
                     n_jobs: int=8,
                     verbose: int=3,
                     seed: int=42,
                     cv_splits: int=5,
                     ) -> Tuple[LinearSVC, pd.DataFrame]:
    """Trains linear support vector machine by searching for optimal l1 ratio.

    train_X - training set features
    train_y - training set targets
    scoring - scikit scoring function to use
    n_jobs - threads to use
    verbose - scikit verbose level
    seed - random seed to use

    returns (model, history dataframe)
    """

    grid = {
        "penalty": ["l1", "l2"],
        "loss": ["hinge", "squared_hinge"],
        "C": [0.001, 0.01, 0.1, 1.0, 10.0, 100.0]
    }
    return train_model(
        LinearSVC(
            dual=False,
            class_weight="balanced",
            random_state=seed,
            max_iter=1000
        ),
        train_gridsearchcv_model,
        train_gridsearchcv_model,
        grid,
        linear_svm_best_params_surroundings,
        train_X,
        train_y,
        scoring=scoring,
        n_jobs=n_jobs,
        verbose=verbose,
        cv_splits=cv_splits
    )

def train_svm_based_on_config(train_X: pd.DataFrame,
                              train_y: pd.Series,
                              test_X: Optional[pd.DataFrame],
                              task_config: TrainingTaskConfig,
                              cv_splits: int,
                              threads: int
                              ) -> Tuple[LinearSVC, Optional[pd.DataFrame]]:

    print("Training: SVM: Starting training.")
    if is_multiclass(train_y):
        metric = "f1_macro"
    else:
        metric = "f1"
    model, _ = train_linear_svm(train_X.values, train_y.values, scoring=metric, cv_splits=cv_splits, n_jobs=threads)
    if task_config.output_model_path:
        print("Training: SVM: Saving model.")
        joblib.dump(model, task_config.output_model_path)

    test_pred = None
    if test_X is not None:
        print("Training: SVM: Predicting test labels.")
        test_pred = model.predict(test_X)
        if task_config.output_prediction_path:
            print("Training: SVM: Saving predicted test labels.")
            pd.Series(test_pred).to_csv(task_config.output_prediction_path)

    return model, test_pred

max_stars_repo_name: adamjanovsky/AndroidMalwareCrypto

max_issues_repo_name: adamjanovsky/AndroidMalwareCrypto

max_forks_repo_name: adamjanovsky/AndroidMalwareCrypto