text:
"""

Originally adapted from https://github.com/cajohare/AxionLimits/
"""

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
import matplotlib.ticker as mticker
from mpl_toolkits.axes_grid1 import make_axes_locatable
import matplotlib.cm as cm
from scipy.stats import norm

from Dirs import axionlimits_dir

import os

rootdir = os.path.dirname(os.path.abspath(__file__)) + "/"
pltdir = rootdir + '../plots/'
pltdir_png = rootdir + '../plots/'

def MySaveFig(fig,pltname,pngsave=False):
    fig.savefig(pltdir+pltname+'.pdf',bbox_inches='tight')
    if pngsave:
        fig.savefig(pltdir_png+pltname+'.png',bbox_inches='tight')

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)]

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

def FigSetup(xlab=r'Dark photon mass $m_{\gamma^\prime}$ [eV]',ylab='Kinetic mixing $\epsilon$',\
             chi_min = 1.0e-18,chi_max = 1.0e0,\
             m_min = 1.0e-15,m_max = 1e5,\
             lw=2.5,lfs=40,tfs=25,tickdir='out',\
             Grid=False,Shape='Rectangular',mathpazo=True,\
             TopAndRightTicks=False,FrequencyAxis=True,FrequencyLabels=True,UnitAxis=True,f_rescale=1,\
            tick_rotation = 20,width=20,height=10,upper_tickdir='out', upper_xlabel=None):

    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=(14.2,14))
    elif Shape=='Custom':
        fig = plt.figure(figsize=(width,height))

    ax = fig.add_subplot(111)

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

    #ax.plot(m_dp,chi100,color="black",zorder=100)
    #ax.plot(m_dp,chi10,color="blue",zorder=100)
    #ax.plot(m_dp,chi1,color="orange",zorder=100)

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

    ax.set_yscale('log')
    ax.set_xscale('log')
    ax.set_xlim([m_min,m_max])
    ax.set_ylim([chi_min,chi_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=tick_rotation)

    if Grid:
        ax.grid(zorder=0)

    if FrequencyAxis:
        ax2 = ax.twiny()

        ax2.set_xscale('log')
        ax2.tick_params(which='major',direction=upper_tickdir,width=2.5,length=13,pad=7, labelsize=tfs)
        ax2.tick_params(which='minor',direction=upper_tickdir,width=1,length=10, labelsize=tfs)
        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)
        ax2.xaxis.set_major_locator(locmaj)
        ax2.xaxis.set_minor_locator(locmin)
        ax2.xaxis.set_minor_formatter(mpl.ticker.NullFormatter())

        if FrequencyLabels:
            ax2.set_xticks([1e0,1e3,1e6,1e9,1e12,1*241.8*1e12,1000*241.8*1e12])
            ax2.set_xticklabels(['Hz','kHz','MHz','GHz','THz','eV','keV'])
        ax2.set_xlim([m_min*241.8*1e12/f_rescale,m_max*241.8*1e12/f_rescale])

        if (upper_xlabel is not None):
            ax2.set_xlabel(upper_xlabel, color='white', fontsize=lfs)

        plt.sca(ax)
    return fig,ax

import matplotlib.patheffects as pe
pek=[pe.Stroke(linewidth=7, foreground='k'), pe.Normal()]

def Haloscopes(ax,col=[0.75, 0.2, 0.2],fs=17,projection=True,text_on=True):
    y2 = ax.get_ylim()[1]
    zo = 0.3

    HAYSTAC_col = 'indianred'
    CAPP_col = 'crimson'
    QUAX_col = 'r'
    ADMX_col = 'firebrick'

    # ADMX
    costh = sqrt(0.0025)
    B = 7.6
    dat = loadtxt(axionlimits_dir + "limit_data/AxionPhoton/ADMX.txt")
    dat[:,1] = 1e-9*dat[:,1]*(B/(1.444e-3*costh*dat[:,0]))
    plt.fill_between(dat[:,0],dat[:,1],y2=y2,edgecolor='none',facecolor=ADMX_col,zorder=0.1,lw=3)

    B = 6.8
    dat = loadtxt(axionlimits_dir + "limit_data/AxionPhoton/ADMX2018.txt")
    dat[:,1] = 1e-9*dat[:,1]*(B/(1.444e-3*costh*dat[:,0]))
    plt.fill_between(dat[:,0],dat[:,1],y2=y2,edgecolor='none',facecolor=ADMX_col,zorder=0.1)

    B = 7.6
    dat = loadtxt(axionlimits_dir + "limit_data/AxionPhoton/ADMX2019_1.txt")
    dat[:,1] = 1e-9*dat[:,1]*(B/(1.444e-3*costh*dat[:,0]))
    plt.fill_between(dat[:,0],dat[:,1],y2=y2,edgecolor='none',facecolor=ADMX_col,zorder=0.1)

    B = 7.6
    dat = loadtxt(axionlimits_dir + "limit_data/AxionPhoton/ADMX2019_2.txt")
    dat[:,1] = 1e-9*dat[:,1]*(B/(1.444e-3*costh*dat[:,0]))
    plt.fill_between(dat[:,0],dat[:,1],y2=y2,edgecolor='none',facecolor=ADMX_col,zorder=0.1)

#     B = 3.11
#     dat = loadtxt(axionlimits_dir + "limit_data/AxionPhoton/ADMX_Sidecar_AC.txt")
#     dat[:,1] = 1e-9*dat[:,1]*(B/(1.444e-3*costh*dat[:,0]))
#     plt.fill_between(dat[:,0],dat[:,1],y2=y2,edgecolor=ADMX_col,facecolor=ADMX_col,zorder=0.1)

#     B = 5.0
#     dat = loadtxt(axionlimits_dir + "limit_data/AxionPhoton/ADMX_SLIC.txt")
#     dat[:,1] = 1e-9*dat[:,1]*(B/(1.444e-3*costh*dat[:,0]))
#     plt.fill_between(dat[:,0],dat[:,1],y2=y2,edgecolor=ADMX_col,facecolor=ADMX_col,zorder=100)

    B = 9
    dat = loadtxt(axionlimits_dir + "limit_data/AxionPhoton/HAYSTAC_highres.txt")
    dat[:,1] = 1e-9*dat[:,1]*(B/(1.444e-3*costh*dat[:,0]))
    plt.fill_between(dat[:,0],dat[:,1],y2=y2,edgecolor='none',facecolor=HAYSTAC_col,zorder=0.1)
    dat = loadtxt(axionlimits_dir + "limit_data/AxionPhoton/HAYSTAC_2020_highres.txt")
    dat[:,1] = 1e-9*dat[:,1]*(B/(1.444e-3*costh*dat[:,0]))
    plt.fill_between(dat[:,0],dat[:,1],y2=y2,edgecolor='none',facecolor=HAYSTAC_col,zorder=0.1)

    # CAPP
    B = 7.3
    dat = loadtxt(axionlimits_dir + "limit_data/AxionPhoton/CAPP-1.txt")
    dat[:,1] = 1e-9*dat[:,1]*(B/(1.444e-3*costh*dat[:,0]))
    plt.fill_between(dat[:,0],dat[:,1],y2=y2,edgecolor='none',facecolor=CAPP_col,zorder=0.1)

    B = 7.8
    dat = loadtxt(axionlimits_dir + "limit_data/AxionPhoton/CAPP-2.txt")
    dat[:,1] = 1e-9*dat[:,1]*(B/(1.444e-3*costh*dat[:,0]))
    plt.fill_between(dat[:,0],dat[:,1],y2=y2,edgecolor='none',facecolor=CAPP_col,zorder=0.1)

    B = 7.9
    dat = loadtxt(axionlimits_dir + "limit_data/AxionPhoton/CAPP-3.txt")
    dat[:,1] = 1e-9*dat[:,1]*(B/(1.444e-3*costh*dat[:,0]))
    plt.fill_between(dat[:,0],dat[:,1],y2=y2,edgecolor='none',facecolor=CAPP_col,zorder=0.1)

    # CAPP-3 [KSVZ]
    dat_min = dat[argmin(dat[:,1]),:]
    dat_min[1] = dat_min[1]*costh/0.11
    plt.plot([dat_min[0],dat_min[0]],[1e-10,dat_min[1]],'-',color=CAPP_col,lw=1.5,zorder=0.1)

    B = 8.1
    dat = loadtxt(axionlimits_dir + "limit_data/AxionPhoton/QUAX.txt")
    dat[:,1] = 1e-9*dat[:,1]*(B/(1.444e-3*0.0046*dat[:,0]))
    plt.fill_between([dat[0,0],dat[0,0]],[y2,dat[0,1]],y2=y2,color=QUAX_col,zorder=0.1)

    if text_on:
        plt.text(1.9e-6,0.5e-14,r'{\bf ADMX}',fontsize=fs,color=ADMX_col,rotation=90,rotation_mode='anchor',ha='center',va='center')
        plt.text(0.8e-5,0.4e-13,r'{\bf CAPP}',fontsize=fs-2,color=CAPP_col,rotation=90,rotation_mode='anchor',ha='center',va='center')
        plt.text(0.19e-4,8e-15,r'{\bf HAYSTAC}',fontsize=fs-5,color=HAYSTAC_col,rotation=90,rotation_mode='anchor',ha='center',va='center')
        plt.text(0.55e-4,4e-11,r'{\bf QUAX}',fontsize=fs-5,color=QUAX_col,rotation=-90,rotation_mode='anchor',ha='center',va='center')

    return

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

    # Globular clusters
    HB_col = [0.01, 0.75, 0.24]
    HB = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/RG.txt")

    if (Higgsed):
        HB = np.concatenate((np.atleast_2d([1e-3, 1.1e-13*(0.1/e_D)]), np.atleast_2d([HB[0,0], 1.1e-13*(0.1/e_D)]), HB))

        #mask = (HB[:,0] < 100) & (HB[:,1] >  1.1e-13)
        #HB[mask,1] = HB[mask,1]*0 + 1.1e-13
    plt.plot(HB[:,0],HB[:,1],color='k',alpha=0.5,zorder=0.9,lw=2)
    plt.fill_between(HB[:,0],HB[:,1],y2=y2,edgecolor=None,facecolor=HB_col,zorder=0.9)

    # Globular clusters
    HB_col = 'DarkGreen'
    HB = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/HB.txt")
    if (Higgsed):
        HB = np.concatenate((np.atleast_2d([1e-3, 2e-13*(0.1/e_D)]), np.atleast_2d([HB[0,0], 2e-13*(0.1/e_D)]), HB))
        #mask = (HB[:,0] < 100) & (HB[:,1] >  2e-13)
        #HB[mask,1] = HB[mask,1]*0 + 2e-13
    plt.plot(HB[:,0],HB[:,1],color='k',alpha=0.5,zorder=0.95,lw=2)
    plt.fill_between(HB[:,0],HB[:,1],y2=y2,edgecolor=None,facecolor=HB_col,zorder=0.95)

    # Solar bound
    Solar_col = 'ForestGreen'
    Solar = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/Solar.txt")
    if (Higgsed):
        mask = (Solar[:,0] < 100) & (Solar[:,1] > 3.4e-13*(0.1/e_D))
        Solar[mask,1] = Solar[mask,1]*0 + 3.4e-13*(0.1/e_D)
    plt.plot(Solar[:,0],Solar[:,1],color='k',alpha=0.5,zorder=1,lw=2)
    plt.fill_between(Solar[:,0],Solar[:,1],y2=y2,edgecolor=None,facecolor=Solar_col,zorder=1)

    if not Higgsed:
        Solar = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/Solar-Global.txt")
        plt.plot(Solar[:,0],Solar[:,1]/Solar[:,0],color='k',alpha=0.5,zorder=1.01,lw=2)
        plt.fill_between(Solar[:,0],Solar[:,1]/Solar[:,0],y2=y2,edgecolor=None,facecolor=Solar_col,zorder=1.01)

    if text_on:
        plt.text(1e2*(1-0.01),1.5e-14*(1+0.05),r'{\bf Solar}',fontsize=fs,color='k',rotation=-39,rotation_mode='anchor',ha='center',va='center')
        plt.text(1e3*(1-0.01),0.7e-14*(1+0.05),r'{\bf HB}',fontsize=fs,color='k',rotation=-36,rotation_mode='anchor',ha='center',va='center')
        plt.text(0.8e4*(1-0.01),0.7e-14*(1+0.05),r'{\bf RG}',fontsize=fs,color='k',rotation=-35,rotation_mode='anchor',ha='center',va='center')
        plt.text(1e2,1.5e-14,r'{\bf Solar}',fontsize=fs,color='w',rotation=-39,rotation_mode='anchor',ha='center',va='center')
        plt.text(1e3,0.7e-14,r'{\bf HB}',fontsize=fs,color='w',rotation=-36,rotation_mode='anchor',ha='center',va='center')
        plt.text(0.8e4,0.7e-14,r'{\bf RG}',fontsize=fs,color='w',rotation=-35,rotation_mode='anchor',ha='center',va='center')
    return

def Xenon(ax,col='crimson',fs=23,text_on=True, f_DP = 1.0):
    y2 = ax.get_ylim()[1]
    dat = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/Xenon1T.txt")
    dat[:,1] = dat[:,1]*sqrt(0.3/(0.45*f_DP))

    plt.plot(1e3*dat[:,0],dat[:,1],color='k',alpha=0.5,zorder=0.5,lw=2)
    plt.fill_between(1e3*dat[:,0],dat[:,1],y2=y2,edgecolor=None,facecolor=col,zorder=0.5)

    dat = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/Xenon1T_S1S2.txt")
    dat[:,1] = dat[:,1]*sqrt(0.3/(0.45*f_DP))

    plt.plot(dat[:,0],dat[:,1],color='k',alpha=0.5,zorder=0.5,lw=2)
    plt.fill_between(dat[:,0],dat[:,1],y2=y2,edgecolor=None,facecolor=col,zorder=0.5)

    dat = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/Xenon100.txt")
    dat[:,1] = dat[:,1]*sqrt(0.3/(0.45*f_DP))

    plt.plot(dat[:,0],dat[:,1],color='k',alpha=0.5,zorder=0.5,lw=2)
    plt.fill_between(dat[:,0],dat[:,1],y2=y2,edgecolor=None,facecolor=col,zorder=0.5)
    if text_on:
        plt.text(8e2,3e-17,r'{\bf XENON}',fontsize=fs,color=col,rotation=0,rotation_mode='anchor',ha='center',va='center')

    return

def DAMIC(ax,col='salmon',fs=21,text_on=True, f_DP = 1.0):
    m1,y1 = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/DM_combined.txt",unpack=True)
    dat = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/DAMIC.txt")
    dat[:,1] = dat[:,1]*sqrt(0.3/(0.45*f_DP))

    y2 = interp(dat[:,0],m1,y1)
    dat[0,1] = y2[0]
    dat[-1,1] = y2[-1]
    plt.plot(dat[:,0],dat[:,1],color='k',alpha=1,zorder=1.05,lw=2)
    plt.fill_between(dat[:,0],dat[:,1],y2=y2,edgecolor=None,facecolor=col,zorder=1.05)
    if text_on:
        plt.text(6e-1,1.5e-14,r'{\bf DAMIC}',fontsize=fs,color=col,rotation=0,rotation_mode='anchor',ha='center',va='center')
        plt.plot([5e0,1e1],[3e-14,6e-14],'k-',lw=2.5,color=col)
    return

def FUNK(ax,col='red',fs=21,text_on=True, f_DP = 1.0):
    m1,y1 = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/DM_combined.txt",unpack=True)
    dat = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/FUNK.txt")
    dat[:,1] = dat[:,1]*sqrt(0.3/(0.45*f_DP))*sqrt(2/3/0.27)

    y2 = interp(dat[:,0],m1,y1)
    dat[0,1] = y2[0]/1.1
    dat[-1,1] = y2[-1]/1.1
    plt.plot(dat[:,0],dat[:,1],color='k',alpha=1,zorder=1.075,lw=2)
    plt.fill_between(dat[:,0],dat[:,1],y2=y2,edgecolor=None,facecolor=col,zorder=1.075)
    if text_on:
        plt.text(5e-1,1e-13,r'{\bf FUNK}',fontsize=fs,color=col,rotation=0,rotation_mode='anchor',ha='center',va='center')
        plt.plot([9e-1,3e0],[3e-13,1e-12],'k-',lw=2.5,color=col)
    return

def SENSEI(ax,col='firebrick',fs=21,text_on=True, f_DP = 1.0):
    #y2 = ax.get_ylim()[1]
    y2 = 1e-10
    dat = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/SENSEI.txt")
    dat[:,1] = dat[:,1]*sqrt(0.3/(0.45*f_DP))

    plt.plot(dat[:,0],dat[:,1],color='k',alpha=1,zorder=1.025,lw=2)
    plt.fill_between(dat[:,0],dat[:,1],y2=y2,edgecolor=None,facecolor=col,zorder=1.025)
    if text_on:
        plt.text(1e0,3e-14,r'{\bf SENSEI}',fontsize=fs,color=col,rotation=0,rotation_mode='anchor',ha='center',va='center')
        plt.plot([2.5e0,4e0],[3e-14,5.3e-14],'k-',lw=2.5,color=col)
    return

def Nanowire(ax,col='pink',fs=22,text_on=True):
    m1,y1 = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/DM_combined.txt",unpack=True)
    dat = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/WSi_Nanowire.txt")
    dat[:,1] = dat[:,1]*sqrt(0.3/0.45)
    y2 = interp(dat[:,0],m1,y1)
    dat[0,1] = y2[0]/1.1
    dat[-1,1] = y2[-1]/1.1
    plt.plot(dat[:,0],dat[:,1],color='k',alpha=1,zorder=1.3,lw=2)
    plt.fill_between(dat[:,0],dat[:,1],y2=y2,edgecolor=None,facecolor=col,zorder=1.3)
    if text_on:
        plt.text(5e-4,1e-10,r'{\bf WSi Nanowire}',fontsize=fs,color=col,rotation=0,rotation_mode='anchor',ha='center',va='center')
        plt.plot([9e-3,3e-3],[3e-10,9e-10],'k-',lw=2.5,color=col)
    return

def Tokyo(ax,col='darkred',fs=15,text_on=True):
    m1,y1 = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/DM_combined.txt",unpack=True)
    dat = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/Tokyo-Dish.txt")
    dat[:,1] = dat[:,1]*sqrt(2/3/0.5)
    y2 = interp(dat[:,0],m1,y1)
    dat[0,1] = y2[0]/1.1
    dat[-1,1] = y2[-1]/1.1
    plt.plot(dat[:,0],dat[:,1],color='k',alpha=1,zorder=1.4,lw=2)
    plt.fill_between(dat[:,0],dat[:,1],y2=y2,edgecolor=None,facecolor=col,zorder=1.4)

    dat = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/Tokyo-Knirck.txt")
    dat[:,1] = dat[:,1]*sqrt(1/3/0.17)
    plt.fill_between(dat[:,0],dat[:,1],y2=1e0,edgecolor='k',facecolor=col,zorder=1.09)

    dat = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/Tokyo-Tomita.txt")
    plt.plot([dat[1,0],dat[1,0]],[dat[1,1],1e0],'-',color=col,lw=3,zorder=0.2)
    if text_on:
        plt.text(2e-4,2.5e-10,r'{\bf Tokyo-3}',fontsize=fs,color=col,rotation=-90,rotation_mode='anchor',ha='center',va='center')
        plt.text(0.5e-3,3e-8,r'{\bf Tokyo-2}',fontsize=fs-2,color='k',rotation=90,rotation_mode='anchor',ha='center',va='center')
        plt.text(0.45e-1,4e-12,r'{\bf Tokyo-1}',fontsize=fs+4,color=col,rotation=0,rotation_mode='anchor',ha='center',va='center')
        plt.plot([3e-1,4e0],[5e-12,8e-12],'-',lw=2.5,color=col)
    return

def Jupiter(ax,col='Green',fs=17,text_on=True):
    y2 = ax.get_ylim()[1]
    dat = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/Jupiter.txt")
    plt.plot(dat[:,0],dat[:,1],color='k',alpha=1,zorder=2,lw=2)
    plt.fill_between(dat[:,0],dat[:,1],y2=y2,edgecolor=None,facecolor=col,zorder=2)
    if text_on:
        plt.text(1e-14*(1-0.02),3e-1*(1+0.07),r'{\bf Jupiter}',fontsize=fs,color='k',rotation=0,rotation_mode='anchor',ha='center',va='center')
        plt.text(1e-14,3e-1,r'{\bf Jupiter}',fontsize=fs,color='w',rotation=0,rotation_mode='anchor',ha='center',va='center')
    return

def Earth(ax,col='DarkGreen',fs=17,text_on=True):
    y2 = ax.get_ylim()[1]
    dat = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/Earth.txt")
    plt.plot(dat[:,0],dat[:,1],color='k',alpha=1,zorder=1.9,lw=2)
    plt.fill_between(dat[:,0],dat[:,1],y2=y2,edgecolor=None,facecolor=col,zorder=1.9)
    if text_on:
        plt.text(4e-13*(1-0.01),2e-1*(1+0.05),r'{\bf Earth}',fontsize=fs,color='k',rotation=0,rotation_mode='anchor',ha='center',va='center')
        plt.text(4e-13,2e-1,r'{\bf Earth}',fontsize=fs,color='w',rotation=0,rotation_mode='anchor',ha='center',va='center')
    return

def Crab(ax,col=[0.1,0.4,0.1],fs=17,text_on=True):
    y2 = ax.get_ylim()[1]
    dat = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/Crab.txt")
    plt.plot(dat[:,0],dat[:,1],color='k',alpha=1,zorder=1.09999,lw=2)
    plt.fill_between(dat[:,0],dat[:,1],y2=y2,edgecolor=None,facecolor=col,zorder=1.09999)

#     dat = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/Crab_2.txt")
#     plt.plot(dat[:,0],dat[:,1],color='k',alpha=1,zorder=1.9,lw=2)
#     plt.fill_between(dat[:,0],dat[:,1],y2=y2,edgecolor=None,facecolor=col,zorder=1.9)
    if text_on:
        plt.text(0.5e-6*(1-0.02),3e-1*(1+0.07),r'{\bf Crab}',fontsize=fs,color='k',rotation=0,rotation_mode='anchor',ha='center',va='center')
        plt.text(0.5e-6,3e-1,r'{\bf Crab}',fontsize=fs,color='w',rotation=0,rotation_mode='anchor',ha='center',va='center')

        plt.text(0.8e-6*(1-0.02),0.9e-1*(1+0.07),r'{\bf nebula}',fontsize=fs,color='k',rotation=0,rotation_mode='anchor',ha='center',va='center')
        plt.text(0.8e-6,0.9e-1,r'{\bf nebula}',fontsize=fs,color='w',rotation=0,rotation_mode='anchor',ha='center',va='center')

    return

def SHUKET(ax,col='maroon',fs=13,text_on=True,edge_on=False,lw=0.8):
    y2 = ax.get_ylim()[1]
    dat = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/SHUKET.txt")
    dat[:,1] = dat[:,1]*sqrt(0.3/0.45)*sqrt(1/3/0.0086)
    plt.fill_between(dat[:,0],dat[:,1],y2=y2,edgecolor=None,facecolor=col,zorder=0.2)
    if edge_on:
        plt.plot(dat[:,0],dat[:,1],'k-',lw=lw,zorder=0.2)
    if text_on:
        plt.text(3.5e-5,1.5e-12,r'{\bf SHUKET}',fontsize=fs,color=col,rotation=-90,rotation_mode='anchor',ha='center',va='center')
    return

def DarkEfield(ax,col='darkred',fs=17,text_on=True,edge_on=False,lw=0.8):
    y2 = ax.get_ylim()[1]
    dat = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/DarkEfield.txt")
    dat[:,1] = dat[:,1]*sqrt(0.3/0.45)*sqrt(1/3/0.079)
    plt.fill_between(dat[:,0],dat[:,1],y2=y2,edgecolor=None,facecolor=col,zorder=0.2)
    if edge_on:
        plt.plot(dat[:,0],dat[:,1],'k-',lw=lw,zorder=0.2)
    if text_on:
        plt.text(0.88e-7,0.2e-12,r'{\bf Dark}',fontsize=fs,color=col,rotation=90,rotation_mode='anchor',ha='center',va='center')
        plt.text(2e-7,0.2e-12,r'{\bf E-field}',fontsize=fs,color=col,rotation=90,rotation_mode='anchor',ha='center',va='center')
    return

def WISPDMX(ax,col='crimson',fs=12,text_on=True,edge_on=False,lw=0.8):
    y2 = ax.get_ylim()[1]
    dat = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/WISPDMX.txt")
    dat[:,1] = dat[:,1]*sqrt(0.3/0.45)*sqrt(1/3/0.079)
    plt.fill_between(dat[:,0],dat[:,1],y2=y2,edgecolor=None,facecolor=col,zorder=0.201)
    if edge_on:
        plt.plot(dat[:,0],dat[:,1],color='k',alpha=1,zorder=0.201,lw=lw)

    if text_on:
        plt.text(9e-7,4.1e-12,r'{\bf WISP}',fontsize=fs,color='w',rotation=0,rotation_mode='anchor',ha='center',va='center')
        plt.text(9e-7,1.8e-12,r'{\bf DMX}',fontsize=fs,color='w',rotation=0,rotation_mode='anchor',ha='center',va='center')

    return

def SQuAD(ax,col=[0.7,0,0],fs=12,text_on=True,lw=0.5,point_on=False,ms=10):
    y2 = ax.get_ylim()[1]
    dat = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/SQuAD.txt")
    dat[:,1] = dat[:,1]*sqrt(0.4/0.45)*sqrt(1/3/0.0025)
    plt.plot([dat[0,0],dat[0,0]],[y2,dat[0,1]],lw=lw,color=col,alpha=1,zorder=0.2)
    if point_on:
        plt.plot(dat[0,0],dat[0,1],'o',mfc=col,mec='k',mew=lw+1,zorder=0.2,markersize=ms)
    if text_on:
        plt.text(36e-6,0.25e-13,r'{\bf SQuAD}',fontsize=fs,color=col,rotation=-90,rotation_mode='anchor',ha='center',va='center')
    return

def DMPathfinder(ax,col='pink',fs=13,text_on=True):
    y2 = ax.get_ylim()[1]
    dat = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/DM-Pathfinder.txt")
    dat[:,1] = dat[:,1]*sqrt(1/0.028)
    plt.plot([dat[0,0],dat[0,0]],[y2,dat[0,1]],lw=2,color=col,alpha=1,zorder=0.6)
    if text_on:
        plt.text(2.1e-9,0.5e-8,r'{\bf DM}',fontsize=fs,color=col,rotation=0,rotation_mode='anchor',ha='center',va='center')
        plt.text(2.1e-9,0.2e-8,r'{\bf Pathfinder}',fontsize=fs,color=col,rotation=0,rotation_mode='anchor',ha='center',va='center')

    return

def DarkMatter(ax,Witte_col='royalblue',Caputo_col='dodgerblue',Arias_col='navy',fs=20,projection=True,text_on=True):
    y2 = ax.get_ylim()[1]
    zo = 0.3

    # Combined limits
    dat = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/DM_combined.txt")
    plt.plot(dat[:,0],dat[:,1],'-',color='w',alpha=1,zorder=zo+0.1,lw=2.5,path_effects=pek)
    plt.fill_between(dat[:,0],dat[:,1],y2=y2,edgecolor=None,facecolor='lightgray',zorder=zo,alpha=1.0)
    plt.plot([1e-16,dat[0,0]],[dat[0,1],dat[0,1]],'--',color='w',alpha=1,zorder=zo+0.1,lw=2.5,path_effects=pek)
    plt.fill_between([1e-16,dat[0,0]],[dat[0,1],dat[0,1]],y2=y2,edgecolor=None,facecolor='lightgray',zorder=zo+0.1,alpha=1.0)
    plt.plot(dat[40:,0],dat[40:,1],'--',color='w',alpha=1,lw=2.5,zorder=1000,solid_capstyle='round')

    # Individual limits
    #dat2 = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/Cosmology_Witte_inhomogeneous.txt")
    #dat4 = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/Cosmology_Caputo_HeII.txt",delimiter=',')
    #dat5 = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/Cosmology_Arias.txt")

    #plt.fill_between(dat2[:,0],dat2[:,1],y2=y2,edgecolor='k',facecolor=Witte_col,zorder=0.305,alpha=0.8)
    #plt.fill_between(dat4[:,0],dat4[:,1],y2=y2,edgecolor='k',facecolor=Caputo_col,zorder=0.305,alpha=0.8)
    #plt.fill_between(dat5[:,0],dat5[:,1],y2=y2,edgecolor='k',facecolor=Arias_col,zorder=0.306,alpha=1)

    """
    if text_on:
        plt.gcf().text(0.21,0.42,r'{\bf DPDM} HeII',fontsize=15,color='w',ha='center')
        plt.gcf().text(0.21,0.4,r'Reionisation',fontsize=15,color='w',ha='center')
        plt.gcf().text(0.21,0.38,r'(Caputo et al.)',fontsize=13,color='w',ha='center')

        plt.gcf().text(0.297,0.37,r'{\bf DPDM}',fontsize=17,color='w',ha='center')
        plt.gcf().text(0.297,0.35,r'(Witte et al.)',fontsize=13,color='w',ha='center')

        plt.gcf().text(0.44,0.48,r'{\bf DPDM}',fontsize=18,color='w',ha='center')
        plt.gcf().text(0.44,0.46,r'(Arias et al.)',fontsize=16,color='w',ha='center')

        plt.arrow(0.24, 0.26, 0, -0.045, transform=fig.transFigure,figure=fig,
          length_includes_head=True,lw=2.5,
          head_width=0.012, head_length=0.028, overhang=0.13,
          edgecolor='k',facecolor='w',clip_on=False)
        plt.text(1e-12,0.4e-16,r'{\bf Dark}',fontsize=27,ha='center')
        plt.text(1e-12,0.7e-17,r'{\bf photon DM}',fontsize=27,ha='center')
    """
    return

def COBEFIRAS(ax,col=[0.1,0.2,0.5],text_on=True):
    y2 = ax.get_ylim()[1]
    dat3 = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/COBEFIRAS.txt",delimiter=',')
    plt.fill_between(dat3[:,0],dat3[:,1],y2=y2,edgecolor='k',facecolor=col,zorder=0.5,alpha=1)
    if text_on:
        plt.gcf().text(0.21,0.72,r'{\bf COBE/FIRAS}',fontsize=22,color='w',ha='center')
        plt.gcf().text(0.21,0.69,r'$\gamma \rightarrow X$',fontsize=22,color='w',ha='center')
    return

def IGM(ax,col='seagreen',fs=17,text_on=True):
    y2 = ax.get_ylim()[1]
    dat = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/IGM.txt")
    plt.plot(dat[:,0],dat[:,1],color='k',alpha=0.5,zorder=0.49,lw=2)

    plt.fill_between(dat[:,0],dat[:,1],y2=y2,edgecolor=None,facecolor=col,zorder=0.49)

    if text_on:
        plt.text(1.2e-12*(1-0.05),0.125e-7*(1+0.07),r'{\bf DPDM} (IGM)',fontsize=fs,color='k',rotation=-38,rotation_mode='anchor',ha='center',va='center')
        plt.text(1.2e-12,0.125e-7,r'{\bf DPDM} (IGM)',fontsize=fs,color='w',rotation=-38,rotation_mode='anchor',ha='center',va='center')
    return

def LeoT(ax,col='mediumseagreen',fs=19,text_on=True):
    y2 = ax.get_ylim()[1]
    dat = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/LeoT.txt")
    plt.plot(dat[:,0],dat[:,1],color='k',alpha=0.48,zorder=0.306,lw=2)
    plt.fill_between(dat[:,0],dat[:,1],y2=y2,edgecolor=None,facecolor=col,zorder=0.306)

    if text_on:
        plt.text(1.1e-13*(1-0.07),1.29e-9*(1+0.1),r'{\bf DPDM} (Leo T)',fontsize=fs,color='k',rotation=-37,rotation_mode='anchor',ha='center',va='center')
        plt.text(1.1e-13,1.29e-9,r'{\bf DPDM} (Leo T)',fontsize=fs,color='w',rotation=-37,rotation_mode='anchor',ha='center',va='center')
    return

def LSW(ax,text_on=True):
    y2 = ax.get_ylim()[1]
    dat = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/SPring-8.txt")
    plt.plot(dat[:,0],dat[:,1],color='k',alpha=1,zorder=1.1001,lw=2)
    plt.fill_between(dat[:,0],dat[:,1],y2=y2,edgecolor=None,facecolor=[0.45, 0.05, 0.1],zorder=1.1001)

    dat = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/ALPS.txt")
    plt.plot(dat[:,0],dat[:,1],color='k',alpha=1,zorder=1.091,lw=2)
    plt.fill_between(dat[:,0],dat[:,1],y2=y2,edgecolor=None,facecolor=[0.55, 0.0, 0.16],zorder=1.091)

    dat = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/LSW_UWA.txt")
    plt.plot(dat[:,0],dat[:,1],color='k',alpha=1,zorder=1.09,lw=2)
    plt.fill_between(dat[:,0],dat[:,1],y2=y2,edgecolor=None,facecolor=[0.6, 0.0, 0.2],zorder=1.09)

    dat = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/LSW_ADMX.txt")
    plt.plot(dat[:,0],dat[:,1],color='k',alpha=1,zorder=1.089,lw=2)
    plt.fill_between(dat[:,0],dat[:,1],y2=y2,edgecolor=None,facecolor=[0.65, 0.1, 0.24],zorder=1.089)

#     dat = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/LSW_CERN.txt")
#     plt.plot(dat[:,0],dat[:,1],color='k',alpha=1,zorder=1.089,lw=2)
#     plt.fill_between(dat[:,0],dat[:,1],y2=y2,edgecolor=None,facecolor=[0.65, 0.15, 0.2],zorder=1.089)

    dat = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/CROWS.txt")
    plt.plot(dat[:,0],dat[:,1],color='k',alpha=1,zorder=1.08,lw=2)
    plt.fill_between(dat[:,0],dat[:,1],y2=y2,edgecolor=None,facecolor=[0.7, 0.2, 0.2],zorder=1.08)

    if text_on:
        plt.text(0.4e-6,0.15e-3,r'{\bf LSW-ADMX}',fontsize=17,color='w',rotation=-55,rotation_mode='anchor',ha='center',va='center')
        plt.text(1e-5,5e-5,r'{\bf LSW-UWA}',fontsize=14,color='w',rotation=-55,rotation_mode='anchor',ha='center',va='center')

        plt.text(0.6e0*(1-0.02),0.9e-4*(1+0.08),r'{\bf LSW-SPring-8}',fontsize=15,color='k',rotation=0,rotation_mode='anchor',ha='center',va='center')
        plt.text(0.6e0,0.9e-4,r'{\bf LSW-SPring-8}',fontsize=15,color='w',rotation=0,rotation_mode='anchor',ha='center',va='center')

        plt.text(1.2e-4*(1-0.02),0.9e-5*(1+0.08),r'{\bf ALPS}',fontsize=25,color='k',rotation=-54,rotation_mode='anchor',ha='center',va='center')
        plt.text(1.2e-4,0.9e-5,r'{\bf ALPS}',fontsize=25,color='w',rotation=-54,rotation_mode='anchor',ha='center',va='center')

        plt.text(0.75e-7*(1-0.01),9.9e-5*(1+0.05),r'{\bf CROWS}',fontsize=24,color='k',rotation=-54,rotation_mode='anchor',ha='center',va='center')
        plt.text(0.75e-7,9.9e-5,r'{\bf CROWS}',fontsize=24,color='w',rotation=-54,rotation_mode='anchor',ha='center',va='center')
    return

def Coulomb(ax,text_on=True):
    y2 = ax.get_ylim()[1]
    dat = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/Cavendish.txt")
    plt.plot(dat[:,0],dat[:,1],color='k',alpha=1,zorder=1.07,lw=2)
    plt.fill_between(dat[:,0],dat[:,1],y2=y2,edgecolor=None,facecolor=[0.7,0,0],zorder=1.07)

    dat = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/PlimptonLawton.txt")
    plt.plot(dat[:,0],dat[:,1],color='k',alpha=1,zorder=1.071,lw=2)
    plt.fill_between(dat[:,0],dat[:,1],y2=y2,edgecolor=None,facecolor='crimson',zorder=1.071)

    dat = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/Spectroscopy.txt")
    plt.plot(dat[:,0],dat[:,1],color='k',alpha=1,zorder=1.11,lw=2)
    plt.fill_between(dat[:,0],dat[:,1],y2=y2,edgecolor=None,facecolor=[0.4, 0.0, 0.13],zorder=1.11)

    dat = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/AFM.txt")
    plt.plot(dat[:,0],dat[:,1],color='k',alpha=1,zorder=1.5,lw=2)
    plt.fill_between(dat[:,0],dat[:,1],y2=y2,edgecolor=None,facecolor=[0.4, 0.2, 0.2],zorder=1.5)
    if text_on:
        plt.text(2e-10*(1-0.02),0.3e-1*(1+0.08),r'{\bf Plimpton-Lawton}',fontsize=18,color='k',rotation=-35,rotation_mode='anchor',ha='center',va='center')
        plt.text(2e-10,0.3e-1,r'{\bf Plimpton-Lawton}',fontsize=18,color='w',rotation=-35,rotation_mode='anchor',ha='center',va='center')

        plt.text(3e1*(1-0.02),3e-1*(1+0.08),r'{\bf AFM}',fontsize=20,color='k',rotation=0,rotation_mode='anchor',ha='center',va='center')
        plt.text(3e1,3e-1,r'{\bf AFM}',fontsize=20,color='w',rotation=0,rotation_mode='anchor',ha='center',va='center')

        plt.text(0.7e-8*(1-0.02),3e-6*(1+0.08),r'{\bf Cavendish-Coulomb}',fontsize=23,color='k',rotation=-35,rotation_mode='anchor',ha='center',va='center')
        plt.text(0.7e-8,3e-6,r'{\bf Cavendish-Coulomb}',fontsize=23,color='w',rotation=-35,rotation_mode='anchor',ha='center',va='center')

        plt.text(0.2e2*(1-0.01),1e-3*(1+0.08),r'{\bf Spectroscopy}',fontsize=23,color='k',rotation=-30,rotation_mode='anchor',ha='center',va='center')
        plt.text(0.2e2,1e-3,r'{\bf Spectroscopy}',fontsize=23,color='w',rotation=-30,rotation_mode='anchor',ha='center',va='center')

    return

def CAST(ax,col='maroon',fs=27,text_on=True):
    y2 = ax.get_ylim()[1]
    dat = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/CAST.txt")
    plt.plot(dat[:,0],dat[:,1],color='k',alpha=1,zorder=1.1,lw=2)
    plt.fill_between(dat[:,0],dat[:,1],y2=y2,edgecolor=None,facecolor=col,zorder=1.1)
    if text_on:
        plt.text(1e3*(1-0.01),0.8e-6*(1+0.08),r'{\bf CAST}',fontsize=fs,color='k',rotation=0,rotation_mode='anchor',ha='center',va='center')
        plt.text(1e3,0.8e-6,r'{\bf CAST}',fontsize=fs,color='w',rotation=0,rotation_mode='anchor',ha='center',va='center')
    return

def SHIPS(ax,col='indianred',fs=20,text_on=True):
    y2 = ax.get_ylim()[1]
    dat = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/SHIPS.txt")
    dat[:,1] = dat[:,1]/dat[:,0]
    plt.plot(dat[:,0],dat[:,1],color='k',alpha=1,zorder=1.09,lw=2)
    plt.fill_between(dat[:,0],dat[:,1],y2=y2,edgecolor=None,facecolor=col,zorder=1.09)
    if text_on:
        plt.text(0.6e-1*(1-0.05),0.08e-8*(1+0.1),r'{\bf SHIPS}',fontsize=fs,color='k',rotation=-30,rotation_mode='anchor',ha='center',va='center')
        plt.text(0.6e-1,0.08e-8,r'{\bf SHIPS}',fontsize=fs,color='w',rotation=-30,rotation_mode='anchor',ha='center',va='center')
    return

def TEXONO(ax,col=[0.5, 0.0, 0.13],fs=15,text_on=True):
    y2 = ax.get_ylim()[1]
    dat = loadtxt(axionlimits_dir + "limit_data/DarkPhoton/TEXONO.txt")
    plt.plot(dat[:,0],dat[:,1],color='k',alpha=1,zorder=1.101,lw=2)
    plt.fill_between(dat[:,0],dat[:,1],y2=y2,edgecolor=None,facecolor=col,zorder=1.101)
    if text_on:
        plt.text(1.2e2*(1-0.01),0.1e-4*(1+0.08),r'{\bf TEXONO}',fontsize=fs,color='k',rotation=0,rotation_mode='anchor',ha='center',va='center')
        plt.text(1.2e2,0.1e-4,r'{\bf TEXONO}',fontsize=fs,color='w',rotation=0,rotation_mode='anchor',ha='center',va='center')
    return

max_stars_repo_name: bradkav/DarkAxionPortal

max_issues_repo_name: bradkav/DarkAxionPortal

max_forks_repo_name: bradkav/DarkAxionPortal

reason: import numpy,from numpy,from scipy


text:
#
# Copyright (C) 2018-2020 Pico Technology Ltd. See LICENSE file for terms.
#
# PS3000 Series (A API) STREAMING MODE EXAMPLE
# This example demonstrates how to call the ps3000a driver API functions in order to open a device, setup 2 channels and collects streamed data (1 buffer).
# This data is then plotted as mV against time in ns.

import ctypes
import numpy as np
from picosdk.ps3000a import ps3000a as ps
import matplotlib.pyplot as plt
from picosdk.functions import adc2mV, assert_pico_ok
import time

# Create chandle and status ready for use
chandle = ctypes.c_int16()
status = {}

# Open PicoScope 5000 Series device
status["openunit"] = ps.ps3000aOpenUnit(ctypes.byref(chandle), None)

try:
    assert_pico_ok(status["openunit"])
except: # PicoNotOkError:

    powerStatus = status["openunit"]

    if powerStatus == 286:
        status["changePowerSource"] = ps.ps3000aChangePowerSource(chandle, powerStatus)
    elif powerStatus == 282:
        status["changePowerSource"] = ps.ps3000aChangePowerSource(chandle, powerStatus)
    else:
        raise

    assert_pico_ok(status["changePowerSource"])

enabled = 1
disabled = 0
analogue_offset = 0.0

# Set up channel A
# handle = chandle
# channel = PS3000A_CHANNEL_A = 0
# enabled = 1
# coupling type = PS3000A_DC = 1
# range = PS3000A_2V = 7
# analogue offset = 0 V
channel_range = ps.PS3000A_RANGE['PS3000A_2V']
status["setChA"] = ps.ps3000aSetChannel(chandle,
                                        ps.PS3000A_CHANNEL['PS3000A_CHANNEL_A'],
                                        enabled,
                                        ps.PS3000A_COUPLING['PS3000A_DC'],
                                        channel_range,
                                        analogue_offset)
assert_pico_ok(status["setChA"])

# Set up channel B
# handle = chandle
# channel = PS3000A_CHANNEL_B = 1
# enabled = 1
# coupling type = PS3000A_DC = 1
# range = PS3000A_2V = 7
# analogue offset = 0 V
status["setChB"] = ps.ps3000aSetChannel(chandle,
                                        ps.PS3000A_CHANNEL['PS3000A_CHANNEL_B'],
                                        enabled,
                                        ps.PS3000A_COUPLING['PS3000A_DC'],
                                        channel_range,
                                        analogue_offset)
assert_pico_ok(status["setChB"])

# Size of capture
sizeOfOneBuffer = 500
numBuffersToCapture = 10

totalSamples = sizeOfOneBuffer * numBuffersToCapture

# Create buffers ready for assigning pointers for data collection
bufferAMax = np.zeros(shape=sizeOfOneBuffer, dtype=np.int16)
bufferBMax = np.zeros(shape=sizeOfOneBuffer, dtype=np.int16)

memory_segment = 0

# Set data buffer location for data collection from channel A
# handle = chandle
# source = PS3000A_CHANNEL_A = 0
# pointer to buffer max = ctypes.byref(bufferAMax)
# pointer to buffer min = ctypes.byref(bufferAMin)
# buffer length = maxSamples
# segment index = 0
# ratio mode = PS3000A_RATIO_MODE_NONE = 0
status["setDataBuffersA"] = ps.ps3000aSetDataBuffers(chandle,
                                                     ps.PS3000A_CHANNEL['PS3000A_CHANNEL_A'],
                                                     bufferAMax.ctypes.data_as(ctypes.POINTER(ctypes.c_int16)),
                                                     None,
                                                     sizeOfOneBuffer,
                                                     memory_segment,
                                                     ps.PS3000A_RATIO_MODE['PS3000A_RATIO_MODE_NONE'])
assert_pico_ok(status["setDataBuffersA"])

# Set data buffer location for data collection from channel B
# handle = chandle
# source = PS3000A_CHANNEL_B = 1
# pointer to buffer max = ctypes.byref(bufferBMax)
# pointer to buffer min = ctypes.byref(bufferBMin)
# buffer length = maxSamples
# segment index = 0
# ratio mode = PS3000A_RATIO_MODE_NONE = 0
status["setDataBuffersB"] = ps.ps3000aSetDataBuffers(chandle,
                                                     ps.PS3000A_CHANNEL['PS3000A_CHANNEL_B'],
                                                     bufferBMax.ctypes.data_as(ctypes.POINTER(ctypes.c_int16)),
                                                     None,
                                                     sizeOfOneBuffer,
                                                     memory_segment,
                                                     ps.PS3000A_RATIO_MODE['PS3000A_RATIO_MODE_NONE'])
assert_pico_ok(status["setDataBuffersB"])

# Begin streaming mode:
sampleInterval = ctypes.c_int32(250)
sampleUnits = ps.PS3000A_TIME_UNITS['PS3000A_US']
# We are not triggering:
maxPreTriggerSamples = 0
autoStopOn = 1
# No downsampling:
downsampleRatio = 1
status["runStreaming"] = ps.ps3000aRunStreaming(chandle,
                                                ctypes.byref(sampleInterval),
                                                sampleUnits,
                                                maxPreTriggerSamples,
                                                totalSamples,
                                                autoStopOn,
                                                downsampleRatio,
                                                ps.PS3000A_RATIO_MODE['PS3000A_RATIO_MODE_NONE'],
                                                sizeOfOneBuffer)
assert_pico_ok(status["runStreaming"])

actualSampleInterval = sampleInterval.value
actualSampleIntervalNs = actualSampleInterval * 1000

print("Capturing at sample interval %s ns" % actualSampleIntervalNs)

# We need a big buffer, not registered with the driver, to keep our complete capture in.
bufferCompleteA = np.zeros(shape=totalSamples, dtype=np.int16)
bufferCompleteB = np.zeros(shape=totalSamples, dtype=np.int16)
nextSample = 0
autoStopOuter = False
wasCalledBack = False

def streaming_callback(handle, noOfSamples, startIndex, overflow, triggerAt, triggered, autoStop, param):
    global nextSample, autoStopOuter, wasCalledBack
    wasCalledBack = True
    destEnd = nextSample + noOfSamples
    sourceEnd = startIndex + noOfSamples
    bufferCompleteA[nextSample:destEnd] = bufferAMax[startIndex:sourceEnd]
    bufferCompleteB[nextSample:destEnd] = bufferBMax[startIndex:sourceEnd]
    nextSample += noOfSamples
    if autoStop:
        autoStopOuter = True

# Convert the python function into a C function pointer.
cFuncPtr = ps.StreamingReadyType(streaming_callback)

# Fetch data from the driver in a loop, copying it out of the registered buffers and into our complete one.
while nextSample < totalSamples and not autoStopOuter:
    wasCalledBack = False
    status["getStreamingLastestValues"] = ps.ps3000aGetStreamingLatestValues(chandle, cFuncPtr, None)
    if not wasCalledBack:
        # If we weren't called back by the driver, this means no data is ready. Sleep for a short while before trying
        # again.
        time.sleep(0.01)

print("Done grabbing values.")

# Find maximum ADC count value
# handle = chandle
# pointer to value = ctypes.byref(maxADC)
maxADC = ctypes.c_int16()
status["maximumValue"] = ps.ps3000aMaximumValue(chandle, ctypes.byref(maxADC))
assert_pico_ok(status["maximumValue"])

# Convert ADC counts data to mV
adc2mVChAMax = adc2mV(bufferCompleteA, channel_range, maxADC)
adc2mVChBMax = adc2mV(bufferCompleteB, channel_range, maxADC)

# Create time data
time = np.linspace(0, (totalSamples) * actualSampleIntervalNs, totalSamples)

# Plot data from channel A and B
plt.plot(time, adc2mVChAMax[:])
plt.plot(time, adc2mVChBMax[:])
plt.xlabel('Time (ns)')
plt.ylabel('Voltage (mV)')
plt.show()

# Stop the scope
# handle = chandle
status["stop"] = ps.ps3000aStop(chandle)
assert_pico_ok(status["stop"])

# Disconnect the scope
# handle = chandle
status["close"] = ps.ps3000aCloseUnit(chandle)
assert_pico_ok(status["close"])

# Display status returns
print(status)

ext: py

max_stars_repo_name: Liambeguin/picosdk-python-wrappers

max_issues_repo_name: Liambeguin/picosdk-python-wrappers

max_forks_repo_name: Liambeguin/picosdk-python-wrappers

reason: import numpy


text:
import random

import numpy as np
import tensorflow as tf

# Height of our input data.
HEIGHT = 10

# The size of each mini-batch.
BATCH_SIZE = 2

# Create a 10x3 matrix in numpy; this lives in main memory (*not* the GPU).
np_data = np.array(range(30), dtype=np.float32).reshape(10, 3)

print('numpy data (in main memory)')
print('---------------------------')
print(np_data)
print()

# Number of epochs to train for.
EPOCHS = 100

# Shuffle the indexes of mini-batches, so that the mini-batches are generated
# in a random order. This helps break locality in the structure of the training
# dataset, which can help with overfitting.
INDEXES = list(range(HEIGHT // BATCH_SIZE))
random.shuffle(INDEXES)

# Copy the numpy data into TF memory as a constant var; this will be copied
# exactly one time into the GPU (if one is available).
tf_data = tf.constant(np_data, dtype=tf.float32)

# The index to use when generating our mini-batch.
ix = tf.placeholder(shape=(), dtype=tf.int32)

# The mini-batch of data we'll work on.
batch = tf.slice(tf_data, [BATCH_SIZE * ix, 0], [BATCH_SIZE, -1])

# The output of the Tensorflow graph.
outp = tf.reduce_sum(tf.square(batch))

with tf.Session() as sess:
    sess.run(tf.global_variables_initializer())
    for epoch in range(EPOCHS):
        for i in INDEXES:
            # Run the computation. The only data in the feed_dict is a single
            # 32-bit integer we supply here. All of the data needed for the
            # mini-batch already lives in GPU memory, and doesn't need to be
            # copied from main memory.
            b, o = sess.run([batch, outp], feed_dict={ix: i})
            print('epoch = {}, ix = {}'.format(epoch, i))
            print('batch: {}'.format(b))
            print('output: {}'.format(o))

max_stars_repo_name: eklitzke/tf-slice

max_issues_repo_name: eklitzke/tf-slice

max_forks_repo_name: eklitzke/tf-slice


text:
# -*- coding:utf-8 -*-
# author: bgl
# weixin dingyuehao: Learn X
# date: 2019/7/1

import argparse
import glob
import os
import uuid

import cv2
import numpy as np
from PIL import Image, ImageDraw, ImageFont

def binarize(im):
    ret, dst = cv2.threshold(im, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
    return dst

def draw_polylines(im, pts):
    pts = np.array(pts, dtype=np.int32).reshape((-1, 1, 2))
    cv2.polylines(im, [pts], True, (0, 255, 255))
    return im

def blur(img, sigma=1, kernel=3):
    return cv2.GaussianBlur(img, (kernel, kernel), sigma)

def warp_affine_point(points, perspective_mat):
    points = np.float32(points).reshape([-1, 2])
    out_point = np.matmul(points, perspective_mat[:, :2].T) + perspective_mat[:, 2].T
    out_point = np.array(out_point, dtype=np.float32).reshape([-1, 2])
    return out_point

def rotate_no_loss(img, angle):
    h, w = img.shape[:2]
    center = (w // 2, h // 2)
    M = cv2.getRotationMatrix2D(center, angle, 1.0)
    in_points = [0, 0, w, 0, w, h, 0, h]
    out_points = warp_affine_point(in_points, M)
    out_points = np.reshape(np.array(out_points, dtype=np.float32), [4, 2])
    minx, miny = np.min(out_points, axis=0)
    w, h = np.max(out_points, axis=0) - np.min(out_points, axis=0)
    out_points -= np.array([minx, miny])

    src = np.reshape(np.float32(in_points), [4, 2])
    dst = np.reshape(np.float32(out_points), [4, 2])
    M = cv2.getPerspectiveTransform(src, dst)
    rotated = cv2.warpPerspective(img, M, (w, h), flags=cv2.INTER_CUBIC, borderMode=cv2.BORDER_CONSTANT)
    return rotated, angle, out_points

def get_text_bmp(text, font, font_size, color):
    font = ImageFont.truetype(font, size=font_size, encoding="utf-8")
    size = font.getsize(text)

    im = np.zeros([size[1], size[0], 4], dtype=np.uint8)
    im[..., 3] = 0
    pilim = Image.fromarray(im)
    draw = ImageDraw.Draw(pilim)

    pos = [0, 0]
    draw.text(pos, text, fill=color, font=font)
    im = np.array(pilim)
    return im

def put_bmp(im, bmp, pos, size):
    bmp = cv2.resize(bmp, tuple(size), cv2.INTER_AREA)
    bmp[..., 3] = binarize(bmp[..., 3])
    target = im[pos[1]:pos[1] + size[1], pos[0]:pos[0] + size[0]]
    for i in range(3):
        target[..., i] = np.where(bmp[..., 3] < 50, target[..., i], bmp[..., i])
    im[pos[1]:pos[1] + size[1], pos[0]:pos[0] + size[0]] = target
    return im

def generate(background_dir, out_dir, count, show=False):
    path_list = glob.glob(os.path.join(background_dir, '*.jpg'))
    path_list += glob.glob(os.path.join(background_dir, '*.png'))
    image_dir = os.path.join(out_dir, 'image')
    label_dir = os.path.join(out_dir, 'label')
    if not os.path.exists(image_dir):
        os.makedirs(image_dir)
    if not os.path.exists(label_dir):
        os.makedirs(label_dir)

    for i in range(count):
        path = path_list[np.random.randint(len(path_list))]

        image = cv2.imread(path)
        h, w = image.shape[:2]
        x = np.random.randint(w // 2)
        y = np.random.randint(h // 2)
        x1 = np.random.randint(w // 2, w)
        y1 = np.random.randint(h // 2, h)
        image = image[y:y1, x:x1]

        image, box_list = put_texts(image)

        name = str(uuid.uuid1())
        label_file = os.path.join(out_dir, 'label', name + '.txt')
        label_file = open(label_file, mode='w', encoding='utf8')
        for box in box_list:
            if show:
                image = draw_polylines(image, box)
            box = np.array(box, dtype=np.float32).reshape([-1, 2]) / np.array(image.shape[:2][::-1])
            box = box.reshape([-1])
            label_file.write(' '.join([str(b) for b in box]) + '\n')

        image_file = os.path.join(out_dir, 'image', name + '.jpg')
        cv2.imwrite(image_file, image)
        if i % 100 == 0:
            print('processed: ', i)

def put_texts(im):
    atoz = ''.join([chr(c) for c in range(ord('a'), ord('z'))])
    AtoZ = ''.join([chr(c) for c in range(ord('A'), ord('Z'))])
    atozAtoZ = atoz + AtoZ
    line_num = 11
    w = 520
    h = 400
    font_list = ['LiberationMono-Bold.ttf', 'NotoMono-Regular.ttf', 'LiberationSans-Bold.ttf', 'Sawasdee-Bold.ttf',
                 'FreeSansBold.ttf']
    im = cv2.resize(im, (w, h))
    line_height = im.shape[0] // line_num
    min_word_per_line = 0
    max_word_per_line = 9
    min_space_width = 20
    max_space_width = 30
    min_font_size = 18
    max_font_size = 25
    min_right_padding = 10
    min_char_per_word = 2
    max_char_per_word = 9

    box_list = []
    for i in range(line_num):
        y = i * line_height
        x = np.random.randint(h // 2)
        for j in range(np.random.randint(min_word_per_line, max_word_per_line)):
            xi = x + np.random.randint(min_space_width, max_space_width)
            word = ''.join([atozAtoZ[np.random.randint(len(atozAtoZ))] for _ in
                            range(np.random.randint(min_char_per_word, max_char_per_word))])
            font_size = np.random.randint(min_font_size, max_font_size)
            font = font_list[np.random.randint(len(font_list))]

            color = (np.random.randint(0, 255), np.random.randint(0, 255), np.random.randint(0, 255),
                     np.random.randint(128, 255))
            yi = y + np.random.randint(line_height // 10)

            bmp = get_text_bmp(word, font, font_size, color)
            bmp, angle, out_points = rotate_no_loss(bmp, np.random.randint(-10, 10))
            s = bmp.shape
            if s[1] + xi > w - min_right_padding:
                continue
            if s[0] + yi > h - min_right_padding:
                continue
            size = (s[1], s[0])
            im = put_bmp(im, bmp, [xi, yi], size)
            four_point = np.reshape(np.array(out_points), [4, 2]) + np.array([xi, yi])

            x = xi + size[0]
            box = four_point.reshape([-1]).tolist()
            box_list.append(box)
    im = blur(im)

    return im, box_list

if __name__ == '__main__':
    parser = argparse.ArgumentParser()
    parser.add_argument("--count", type=int, default=2000)
    parser.add_argument("--back_dir", required=True, help='background image directory')
    parser.add_argument("--out_dir", required=True, help='output result directory')
    args = parser.parse_args()
    count = args.count
    show = False

    generate(args.back_dir, args.out_dir, args.count, show=show)


text:
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
from datetime import datetime

class DataReader():
    def __init__(self, face_data, person_data):
        """Process raw object detection data and provides different plot
        outputs

        Arguments:
            face_data (df)
            person_data (df)

        Attributes:

        """
        self.raw_face_data = face_data
        self.raw_person_data = person_data
        self.clean_face_data = self.process_data(face_data)
        self.clean_person_data = self.process_data(person_data)
        self.start_date = None
        self.end_date = None
        self.x_axis_ts = self.get_timeframe()
        self.x_axis_timeofday = self.ts_to_timeofday(self.x_axis_ts)

        self.htr_data = self.get_htr_data()

    def process_data(self, raw_df):
        raw_gb = raw_df.groupby(by=['ts', 'id'])
        raw_gb_mean = raw_gb.mean().reset_index()

        # add object centroid to the df
        raw_gb_mean['X'] = raw_gb_mean.apply(
            lambda x: self.get_object_centroid(x, 'X'), axis=1)
        raw_gb_mean['Y'] = raw_gb_mean.apply(
            lambda x: self.get_object_centroid(x, 'Y'), axis=1)

        # add time of day column so that we don't need to read ts
        raw_gb_mean = self.add_timeofday_col(raw_gb_mean)
        return raw_gb_mean

    def get_object_centroid(self, row, axis='X'):
        """Iterates through each row to return the average of start and end
        Pandas apply function
        """
        if axis == 'X':
            return (row['startX'] + row['endX']) / 2
        elif axis == 'Y':
            return (row['startY'] + row['endY']) / 2
        else:
            print("'axis' must be 'X' or 'Y'")
            raise

    def add_timeofday_col(self, df):
        df['dt'] = pd.to_datetime(df['ts'], unit='s')
        df['time_of_day'] = df['dt'].apply(lambda x: x.strftime('%H:%M:%S'))
        df = df.drop('dt', axis=1)
        return df

    def ts_to_timeofday(self, l_ts):
        timeofday = []
        for ts in l_ts:
            timeofday.append(
                datetime.utcfromtimestamp(ts).strftime('%H:%M:%S'))
        assert len(timeofday) == len(l_ts)
        return timeofday

    def get_timeframe(self):
        """Returns min and max ts from processed data"""
        min_ts_face = self.clean_face_data['ts'].min()
        min_ts_person = self.clean_person_data['ts'].min()
        max_ts_face = self.clean_face_data['ts'].max()
        max_ts_person = self.clean_person_data['ts'].min()
        return np.arange(min(min_ts_face, min_ts_person),
                         max(max_ts_face, max_ts_person), 1)

    def get_htr_data(self):
        count_faces = self.clean_face_data.groupby(by='ts').nunique()
        count_persons = self.clean_person_data.groupby(by='ts').nunique()

        infill_count_faces = [count_faces.loc[ts, 'id']
                              if ts in count_faces.index
                              else 0 for ts in self.x_axis_ts]
        infill_count_persons = [count_persons.loc[ts, 'id']
                                if ts in count_persons.index
                                else 0 for ts in self.x_axis_ts]

        htr_df = pd.DataFrame(data={
            'faces': infill_count_faces,
            'persons': infill_count_persons}, index=self.x_axis_timeofday)

        return htr_df

    def save_htr_data(self, path=None, method='csv'):
        """Write processed data to csv or sql."""
        pass

    def stackarea(self, figsize=(6, 6), timestep=5):
        """Returns a stacked area plot with persons on the bottom"""
        fig, ax = plt.subplots(figsize=figsize)
        X = self.x_axis_timeofday
        Y1 = self.htr_data['persons']
        Y2 = self.htr_data['faces']
        ax.stackplot(X, np.array([Y1, Y2]), labels=['persons', 'faces'],
                     edgecolor='white')
        xticks = ax.get_xticks()
        ax.set_xticks(xticks[::timestep])
        plt.legend()
        plt.show()

    def stackbar(self, figsize=(6, 6), timestep=5):
        """Returns a stacked bar plot with persons on the bottom"""
        fig, ax = plt.subplots(figsize=figsize)
        X = self.x_axis_timeofday
        Y1 = self.htr_data['persons']
        Y2 = self.htr_data['faces']
        ax.bar(X, Y1, color='#5392ff')
        ax.bar(X, Y2, color='#ff5c49', bottom=Y1)
        xticks = ax.get_xticks()
        ax.set_xticks(xticks[::timestep])
        plt.show()

    def plotbar(self, metric='face', figsize=(6, 6), timestep=5):
        """

        Arguments:
            metric (str): 'face' or 'persons'
        """
        fig, ax = plt.subplots(figsize=figsize)
        X = self.x_axis_timeofday

        if metric == 'face':
            Y = self.htr_data['faces']
        elif metric == 'person':
            Y = self.htr_data['persons']
        else:
            print("'metric' must be 'face' or 'person'")
            raise
        ax.bar(X, Y, color='#5392ff')
        xticks = ax.get_xticks()
        ax.set_xticks(xticks[::timestep])
        plt.show()

    def plotline(self, metric='face', figsize=(6, 6), timestep=5):
        """

        Arguments:
            metric (str): 'face' or 'persons'
        """
        fig, ax = plt.subplots(figsize=figsize)
        X = self.x_axis_timeofday
        if metric == 'face':
            Y = self.htr_data['faces']
        elif metric == 'person':
            Y = self.htr_data['persons']
        else:
            print("'metric' must be 'face' or 'person'")
            raise
        ax.plot(X, Y, color='#5392ff')
        xticks = ax.get_xticks()
        ax.set_xticks(xticks[::timestep])
        plt.show()


text:
# MIT License - Copyright Petri Laarne and contributors
# See the LICENSE.md file included in this source code package

"""Tests for ennemi._estimate_single_mi() and friends."""

import math
from math import log
import numpy as np
from scipy.special import gamma, psi
from scipy.stats import gamma as gamma_dist
import unittest
from ennemi._entropy_estimators import _estimate_single_entropy,\
    _estimate_single_mi, _estimate_conditional_mi,\
    _estimate_semidiscrete_mi, _estimate_conditional_semidiscrete_mi

class TestEstimateSingleEntropy(unittest.TestCase):

    def test_univariate_gaussian(self) -> None:
        cases = [ (1, 100, 2, 0.2),
                  (1, 200, 3, 0.05),
                  (1, 2000, 3, 0.02),
                  (0.5, 2000, 3, 0.02),
                  (2.0, 2000, 1, 0.002),
                  (2.0, 2000, 3, 0.02),
                  (2.0, 2000, 30, 0.02), ]
        for (sd, n, k, delta) in cases:
            with self.subTest(sd=sd, n=n, k=k):
                rng = np.random.default_rng(0)
                x = rng.normal(0, sd, size=n)

                actual = _estimate_single_entropy(x, k=k)
                expected = 0.5 * log(2 * math.pi * math.e * sd**2)
                self.assertAlmostEqual(actual, expected, delta=delta)

    def test_uniform(self) -> None:
        cases = [ (0, 1, 1000, 3, 0.05),
                  (1, 2, 1000, 3, 0.05),
                  (-1, 1, 1000, 3, 0.05),
                  (-0.1, 0.1, 1000, 3, 0.05), ]
        for (a, b, n, k, delta) in cases:
            with self.subTest(a=a, b=b, n=n, k=k):
                rng = np.random.default_rng(1)
                x = rng.uniform(a, b, size=n)

                actual = _estimate_single_entropy(x, k=k)
                expected = log(b - a)
                self.assertAlmostEqual(actual, expected, delta=delta)

    def test_bivariate_gaussian(self) -> None:
        cases = [ (0, 1, 200, 3, 0.09),
                  (0, 2, 2000, 3, 0.03),
                  (0.2, 1, 2000, 3, 0.03),
                  (0.2, 2, 2000, 5, 0.03),
                  (0.6, 1, 2000, 1, 0.02),
                  (0.6, 0.5, 2000, 3, 0.04),
                  (0.9, 1, 2000, 3, 0.04),
                  (-0.5, 1, 2000, 5, 0.03), ]
        for (rho, var1, n, k, delta) in cases:
            with self.subTest(rho=rho, var1=var1, n=n, k=k):
                rng = np.random.default_rng(2)
                cov = np.array([[var1, rho], [rho, 1]])
                data = rng.multivariate_normal([0, 0], cov, size=n)

                actual = _estimate_single_entropy(data, k=k)
                expected = 0.5 * log(np.linalg.det(2 * math.pi * math.e * cov))
                self.assertAlmostEqual(actual, expected, delta=delta)

    def test_4d_gaussian(self) -> None:
        rng = np.random.default_rng(3)
        cov = np.array([
            [ 1.0,  0.5,  0.6, -0.2],
            [ 0.5,  1.0,  0.7, -0.5],
            [ 0.6,  0.7,  2.0, -0.1],
            [-0.2, -0.5, -0.1,  0.5]])
        data = rng.multivariate_normal([0, 0, 0, 0], cov, size=2000)

        actual = _estimate_single_entropy(data, k=3)
        expected = 0.5 * log(np.linalg.det(2 * math.pi * math.e * cov))
        self.assertAlmostEqual(actual, expected, delta=0.05)

    def test_gamma_exponential(self) -> None:
        # As in the MI test, the analytical result is due to doi:10.1109/18.825848.
        #
        # x1      ~ Gamma(rate, shape)
        # x2 | x1 ~ Exp(t * x1)
        rng = np.random.default_rng(4)
        r = 1.2
        s = 3.4
        t = 0.56

        x1 = rng.gamma(shape=s, scale=1/r, size=1000)
        x2 = rng.exponential(x1 * t)
        data = np.asarray([x1, x2]).T

        raw = _estimate_single_entropy(data)
        trans = _estimate_single_entropy(np.log(data))

        # The estimate with unlogarithmed data is very bad
        expected = 1 + s - s*psi(s) + log(gamma(s)) - log(t)
        self.assertAlmostEqual(raw, expected, delta=0.65)
        self.assertAlmostEqual(trans, expected, delta=0.01)

class TestEstimateSingleMi(unittest.TestCase):

    def test_bivariate_gaussian(self) -> None:
        cases = [ (0, 40, 3, 0.1),
                  (0, 200, 3, 0.06),
                  (0, 2000, 3, 0.005),
                  (0, 2000, 5, 0.006),
                  (0, 2000, 20, 0.003),
                  (0.5, 200, 3, 0.05),
                  (0.5, 200, 5, 0.02),
                  (0.5, 2000, 3, 0.02),
                  (-0.9, 200, 3, 0.05),
                  (-0.9, 2000, 3, 0.05),
                  (-0.9, 2000, 5, 0.02), ]
        for (rho, n, k, delta) in cases:
            with self.subTest(rho=rho, n=n, k=k):
                rng = np.random.default_rng(0)
                cov = np.array([[1, rho], [rho, 1]])

                data = rng.multivariate_normal([0, 0], cov, size=n)
                x = data[:,0]
                y = data[:,1]

                actual = _estimate_single_mi(x, y, k=k)
                expected = -0.5 * log(1 - rho**2)
                self.assertAlmostEqual(actual, expected, delta=delta)

    def test_sum_of_exponentials(self) -> None:
        # We define X ~ Exp(1), W ~ Exp(2) and Y = X + W.
        # Now by arXiv:1609.02911, Y has known, closed-form entropy.
        cases = [ (1, 2), (0.2, 0.3), (3, 3.1) ]
        for (a, b) in cases:
            with self.subTest(a=a, b=b):
                rng = np.random.default_rng(20200302)
                x = rng.exponential(1/a, 1000)
                w = rng.exponential(1/b, 1000)
                y = x + w

                actual = _estimate_single_mi(x, y, k=5)
                expected = np.euler_gamma + log((b-a)/a) + psi(b/(b-a))

                self.assertAlmostEqual(actual, expected, delta=0.025)

    def test_independent_uniform(self) -> None:
        # We have to use independent random numbers instead of linspace,
        # because the algorithm has trouble with evenly spaced values
        rng = np.random.default_rng(1)
        x = rng.uniform(0.0, 1.0, 1024)
        y = rng.uniform(0.0, 1.0, 1024)

        actual = _estimate_single_mi(x, y, k=8)
        actual2 = _estimate_single_mi(y, x, k=8)
        self.assertAlmostEqual(actual, 0, delta=0.04)
        self.assertAlmostEqual(actual, actual2, delta=0.00001)

    def test_independent_transformed_uniform(self) -> None:
        # Very non-uniform density, but MI should still be zero
        rng = np.random.default_rng(1)
        x = rng.uniform(0.0, 10.0, 1024)
        y = np.exp(rng.uniform(0.0, 1.0, 1024))

        actual = _estimate_single_mi(x, y, k=8)
        self.assertAlmostEqual(actual, 0, delta=0.02)

    def test_gamma_exponential(self) -> None:
        # Kraskov et al. mention that this distribution is hard to estimate
        # without logarithming the values.
        # The analytical result is due to doi:10.1109/18.825848.
        #
        # x1      ~ Gamma(rate, shape)
        # x2 | x1 ~ Exp(t * x1)
        rng = np.random.default_rng(2)
        r = 1.2
        s = 3.4
        t = 0.56

        x1 = rng.gamma(shape=s, scale=1/r, size=1000)
        x2 = rng.exponential(x1 * t)

        raw = _estimate_single_mi(x1, x2)
        trans = _estimate_single_mi(np.log(x1), np.log(x2))

        expected = psi(s) - np.log(s) + 1/s
        self.assertAlmostEqual(raw, expected, delta=0.04)
        self.assertAlmostEqual(trans, expected, delta=0.005)

class TestEstimateConditionalMi(unittest.TestCase):

    def test_gaussian_with_independent_condition(self) -> None:
        # In this case, the results should be same as in ordinary MI
        cases = [ (0.5, 200, 3, 0.03),
                  (0.75, 400, 3, 0.01),
                  (-0.9, 4000, 5, 0.03), ]
        for (rho, n, k, delta) in cases:
            with self.subTest(rho=rho, n=n, k=k):
                rng = np.random.default_rng(0)
                cov = np.array([[1, rho], [rho, 1]])

                data = rng.multivariate_normal([0, 0], cov, size=n)
                x = data[:,0]
                y = data[:,1]
                cond = rng.uniform(0, 1, size=n)

                actual = _estimate_conditional_mi(x, y, cond, k=k)
                expected = -0.5 * log(1 - rho**2)
                self.assertAlmostEqual(actual, expected, delta=delta)

    def test_gaussian_with_condition_equal_to_y(self) -> None:
        # MI(X;Y | Y) should be equal to 0
        rng = np.random.default_rng(4)
        cov = np.array([[1, 0.6], [0.6, 1]])

        data = rng.multivariate_normal([0, 0], cov, size=314)
        x = data[:,0]
        y = data[:,1]

        actual = _estimate_conditional_mi(x, y, y, k=4)
        self.assertAlmostEqual(actual, 0.0, delta=0.001)

    def test_three_gaussians(self) -> None:
        # First example in doi:10.1103/PhysRevLett.99.204101,
        # we know the analytic expression for conditional MI of a three-
        # dimensional Gaussian random variable. Here, the covariance matrix
        # is not particularly interesting. The expected CMI expression
        # contains determinants of submatrices.
        rng = np.random.default_rng(5)
        cov = np.array([[1, 1, 1], [1, 4, 1], [1, 1, 9]])

        data = rng.multivariate_normal([0, 0, 0], cov, size=1000)

        actual = _estimate_conditional_mi(data[:,0], data[:,1], data[:,2])
        expected = 0.5 * (log(8) + log(35) - log(9) - log(24))
        self.assertAlmostEqual(actual, expected, delta=0.015)

    def test_four_gaussians(self) -> None:
        # As above, but now the condition is two-dimensional.
        # The covariance matrix is defined by transforming a standard normal
        # distribution (u1, u2, u3, u4) as follows:
        #   x  = u1,
        #   y  = u2 + u3 + 2*u4,
        #   z1 = 2*u1 + u3,
        #   z2 = u1 + u4.
        # Unconditionally, x and y are independent, but conditionally they aren't.
        rng = np.random.default_rng(25)
        cov = np.array([[1, 0, 2, 1],
                        [0, 6, 1, 2],
                        [2, 1, 5, 2],
                        [1, 2, 2, 2]])

        # The data needs to be normalized for estimation accuracy,
        # and the sample size must be quite large
        data = rng.multivariate_normal([0, 0, 0, 0], cov, size=8000)
        data = data / np.sqrt(np.var(data, axis=0))

        actual = _estimate_conditional_mi(data[:,0], data[:,1], data[:,2:])
        expected = 0.64964
        self.assertAlmostEqual(actual, expected, delta=0.04)

class TestEstimateSemiDiscreteMi(unittest.TestCase):

    def test_independent_variables(self) -> None:
        cases = [ (2, 200, 3, 0.04),
                  (2, 400, 1, 0.02),
                  (2, 400, 3, 0.02),
                  (2, 800, 8, 0.02),
                  (4, 2000, 2, 0.01) ]
        for (discrete_count, n, k, delta) in cases:
            with self.subTest(count=discrete_count, n=n, k=k):
                rng = np.random.default_rng(50)
                x = rng.normal(0.0, 1.0, size=n)
                y = rng.choice(np.arange(discrete_count), size=n)

                mi = _estimate_semidiscrete_mi(x, y, k)
                self.assertAlmostEqual(max(mi, 0.0), 0.0, delta=delta)

    def test_two_disjoint_uniforms(self) -> None:
        # Y takes two equally probable values, and then X is sampled
        # from two disjoint distributions depending on Y.
        # Therefore I(X;Y) = H(Y) = log(2).
        rng = np.random.default_rng(51)
        y = rng.choice([0, 2], size=800)
        x = rng.uniform(y, y+1)

        mi = _estimate_semidiscrete_mi(x, y)
        self.assertAlmostEqual(mi, log(2), delta=0.02)

    def test_three_disjoint_uniforms(self) -> None:
        # As above, but with three equally probable values for Y.
        rng = np.random.default_rng(51)
        y = rng.choice([0, 2, 5], size=800)
        x = rng.uniform(y, y+1)

        mi = _estimate_semidiscrete_mi(x, y)
        self.assertAlmostEqual(mi, log(3), delta=0.02)

    def test_two_overlapping_uniforms(self) -> None:
        # Here there are two values for Y, but the associated X intervals overlap.
        # Additionally, one of the values is more likely than the other.
        rng = np.random.default_rng(52)
        y = rng.choice([0, 0.7, 0.7], size=2000)
        x = rng.uniform(y, y+1)

        mi = _estimate_semidiscrete_mi(x, y)
        expected = log(3)*7/30 + log(1)*9/30 + log(3/2)*14/30
        self.assertAlmostEqual(mi, expected, delta=0.05)

class TestEstimateConditionalSemiDiscreteMi(unittest.TestCase):

    def test_independent_vars_and_condition(self) -> None:
        cases = [ (2, 200, 3, 0.02),
                  (2, 400, 1, 0.05),
                  (2, 400, 3, 0.02),
                  (2, 800, 8, 0.02),
                  (4, 2000, 2, 0.02) ]
        for (discrete_count, n, k, delta) in cases:
            with self.subTest(count=discrete_count, n=n, k=k):
                rng = np.random.default_rng(50)
                x = rng.normal(0.0, 1.0, size=n)
                y = rng.choice(np.arange(discrete_count), size=n)
                z = rng.normal(0.0, 1.0, size=n)

                mi = _estimate_conditional_semidiscrete_mi(x, y, z, k)
                self.assertAlmostEqual(max(mi, 0.0), 0.0, delta=delta)

    def test_irrelevant_condition(self) -> None:
        # Two disjoint uniforms, but with independent condition
        rng = np.random.default_rng(51)
        y = rng.choice([0, 2], size=800)
        x = rng.uniform(y, y+1)
        z = rng.beta(2, 3, size=800)

        mi = _estimate_conditional_semidiscrete_mi(x, y, z)
        self.assertAlmostEqual(mi, log(2), delta=0.02)

    def test_condition_equal_to_x(self) -> None:
        # Y = sign(X), and I(Y;X | X) = 0
        rng = np.random.default_rng(52)
        x = rng.normal(0.0, 1.0, size=800)
        y = np.sign(x)

        # Consistency check: the unconditional MI should be equal to y entropy
        uncond = _estimate_semidiscrete_mi(x, y)
        self.assertAlmostEqual(uncond, log(2), delta=0.01)

        mi = _estimate_conditional_semidiscrete_mi(x, y, x)
        self.assertAlmostEqual(mi, 0.0, delta=0.02)

    def test_condition_increases_mi(self) -> None:
        # X, Y are normal, Z = X + Y
        # W = sign(X)
        # I(Z; W) < I(X; W), but I(Z; W | Y) = I(X; W)
        rng = np.random.default_rng(53)
        x = rng.normal(0.0, 1.0, size=1500)
        y = rng.normal(0.0, 0.3, size=1500)
        z = x + y
        w = np.sign(x)

        # Consistency check
        uncond_zw = _estimate_semidiscrete_mi(z, w)
        uncond_xw = _estimate_semidiscrete_mi(x, w)
        self.assertLess(uncond_zw, uncond_xw - 0.2)
        self.assertAlmostEqual(uncond_xw, log(2), delta=0.01)

        # Conditioning should increase MI
        cond_zw = _estimate_conditional_semidiscrete_mi(z, w, y)
        self.assertAlmostEqual(cond_zw, log(2), delta=0.02)

    def test_multiple_levels(self) -> None:
        # X is a Gamma random variable, Y takes levels based on X,
        # and W is X with added Gaussian noise. Conditioning on noise increases MI.
        rng = np.random.default_rng(54)
        x = rng.gamma(shape=1.5, scale=1.0, size=2000)
        z = rng.normal(size=x.shape)
        w = x + z

        # The 1e-4 level would cause issues with the continuous-continuous algorithm
        # as it would be picked up in neighbor searches on the y=0 plane
        y = np.zeros(x.shape)
        y[x < 0.5] = 1e-4
        y[x > 2.0] = 4

        uncond = _estimate_semidiscrete_mi(w, y, k=1)
        cond = _estimate_conditional_semidiscrete_mi(w, y, z, k=1)

        # The expected MI is the discrete entropy of Y
        p_low = gamma_dist.cdf(0.5, a=1.5)
        p_high = 1 - gamma_dist.cdf(2.0, a=1.5)
        p_mid = 1 - p_low - p_high
        expected = -log(p_low)*p_low - log(p_mid)*p_mid - log(p_high)*p_high

        self.assertLess(uncond, cond - 0.6)
        self.assertAlmostEqual(cond, expected, delta=0.06)

# Test our custom implementation of the digamma function
from ennemi._entropy_estimators import _psi

class TestPsi(unittest.TestCase):

    def test_comparison_with_scipy(self) -> None:
        # Zero
        self.assertEqual(_psi(0), np.inf)

        # Small values
        for x in range(1, 20):
            with self.subTest(x=x):
                self.assertAlmostEqual(_psi(x), psi(x), delta=0.0001)

        # Large values
        for x in range(2, 1000):
            with self.subTest(x=10*x):
                self.assertAlmostEqual(_psi(x*10), psi(x*10), delta=0.0001)

max_stars_repo_name: atentas/ennemi

max_issues_repo_name: atentas/ennemi

max_forks_repo_name: atentas/ennemi

reason: import numpy,from scipy


text:
import numpy as np
from teilab.utils import subplots_create
from teilab.statistics import mann_whitney_u_test
fig, axes = subplots_create(ncols=3, figsize=(18,4), style="matplotlib")
A = np.array([1.83, 1.50, 1.62, 2.48, 1.68, 1.88, 1.55, 3.06, 1.30, 2.01, 3.11])
B = np.array([0.88, 0.65, 0.60, 1.05, 1.06, 1.29, 1.06, 2.14, 1.29])
for ax,alternative in zip(axes,["less","two-sided","greater"]):
    mann_whitney_u_test(A, B, alternative=alternative, plot=True, alpha=.1, ax=ax)
fig.show()


text:
import cv2
import numpy as np
import time

now = time.time()
img2 = cv2.imread('output.png')

str=''
arr=[]
arr1 =[]
c=0
def two_dim_conv(x):
    c=0
    for i in x:
        for j in i:
            for k in j:
                if c>img2[0][0][0]:
                    break
                arr.append(k)
                c+=1
    return arr

two_dim_conv(img2)
l = arr[0]
string_arr = arr[1:l+1]

for i in string_arr:
    str += chr(i)

print(str)
print("Intercercepted message length >>> ",end=' ')
print(img2[0][0][0])
print("Total time taken to Decode >>> ",end=' ')
print(time.time()-now)


text:
import argparse
import copy
import datetime
import os
import random
import socket
import string
import sys
import cv2
import matplotlib.pyplot as plt
import numpy as np
import torch
import GPUtil

def get_gpus(maxLoad=0.5, maxMemory=0.5):
    """Returns the available GPUs."""
    deviceIDs = GPUtil.getAvailable(
        order='first',
        limit=8,
        maxLoad=maxLoad,
        maxMemory=maxMemory,
        includeNan=False,
        excludeID=[],
        excludeUUID=[])
    return deviceIDs

def get_chunks(lst, n):
    """Returns list of n elements, constaining a sublist."""
    size = len(lst) // n
    chunks = []
    for i in range(0, len(lst), size):
        chunks.append(lst[i:i + size])
    return chunks

def to_numpy(x):
    if isinstance(x, torch.Tensor):
        return x.detach().cpu().numpy()
    else:
        return x

def to_torch(x):
    return torch.FloatTensor(x)

def dict_to_torch_device(stuff, torch_device):
    """Set everything in the dict (2 layers deep) to a specified torch device.
    """
    if type(stuff) is dict:
        for k, v in stuff.items():
            stuff[k] = dict_to_torch_device(v, torch_device)
    if type(stuff) is torch.Tensor:
        return stuff.to(torch_device)
    else:
        return stuff

def dict_to_cpu(stuff):
    if type(stuff) is dict:
        for k, v in stuff.items():
            stuff[k] = dict_to_cpu(v)
    if type(stuff) is torch.Tensor:
        return stuff.detach().cpu()
    else:
        return stuff

def dict_to_torch(d):
    for k, v in d.items():
        if isinstance(v, dict):
            d[k] = dict_to_torch(v)
        elif isinstance(v, (np.ndarray, np.generic)):
            d[k] = torch.from_numpy(v)
    return d

def dict_to_numpy(d):
    for k, v in d.items():
        if isinstance(v, dict):
            d[k] = dict_to_numpy(v)
        elif isinstance(v, torch.Tensor):
            d[k] = to_numpy(v)
    return d

class dotdict(dict):
    """dot.notation access to dictionary attributes"""
    __getattr__ = dict.get
    __setattr__ = dict.__setitem__
    __delattr__ = dict.__delitem__

def get_random_string(length=8):
    letters = string.ascii_lowercase
    result_str = ''.join(random.choice(letters) for i in range(length))
    return result_str

def gettimedatestring():
    return datetime.datetime.now().strftime("%m-%d-%H-%M-%S")

def select_gpus(gpus_arg):
    # so that default gpu is one of the selected, instead of 0
    if len(gpus_arg) > 0:
        os.environ['CUDA_VISIBLE_DEVICES'] = gpus_arg
        gpus = list(range(len(gpus_arg.split(','))))
    else:
        os.environ['CUDA_VISIBLE_DEVICES'] = ''
        gpus = []
    print('CUDA_VISIBLE_DEVICES={}'.format(os.environ['CUDA_VISIBLE_DEVICES']))
    return gpus

def str2bool(v):
    assert type(v) is str
    if v.lower() in ('yes', 'true', 't', 'y', '1'):
        return True
    elif v.lower() in ('no', 'false', 'f', 'n', '0'):
        return False
    else:
        raise argparse.ArgumentTypeError('Boolean (yes, true, t, y or 1, lower or upper case) string expected.')

def listdir(folder, prepend_folder=False, extension=None, type=None):
    assert type in [None, 'file', 'folder'], "Type must be None, 'file' or 'folder'"
    files = [k for k in os.listdir(folder) if (True if extension is None else k.endswith(extension))]
    if type == 'folder':
        files = [k for k in files if os.path.isdir(folder + '/' + k)]
    elif type == 'file':
        files = [k for k in files if not os.path.isdir(folder + '/' + k)]
    if prepend_folder:
        files = [folder + '/' + f for f in files]
    return files

def get_hostname():
    return socket.gethostname()

def stop():
    import sys
    print("called goat.stop()")
    sys.exit()

def pose_to_homo(pose):
    """TODO(ethan): write this for torch too
    """
    if pose.shape == (3, 4):
        pose = np.concatenate([pose, np.zeros_like(pose[:1])], -2)
        pose[3, 3] = 1
    assert pose.shape == (4,4)
    return pose

max_stars_repo_name: ethanweber/goat

max_issues_repo_name: ethanweber/goat

max_forks_repo_name: ethanweber/goat


text:
import argparse
import torch
import os
import platform
import sys
import math
import cv2
import numpy as np
from api.setup import SingleImageAlphaPose, get_args

def evaluate(save_img, save_keypts, save_json):
	outputpath = "examples/res/"
	if not os.path.exists(outputpath + '/vis'):
		os.makedirs(outputpath + '/vis')

	# Create Network.
	args, cfg = get_args()
	demo = SingleImageAlphaPose(args, cfg)

	# Give path of the image.
	im_name = 'examples/demo/80.png'
	image = cv2.cvtColor(cv2.imread(im_name), cv2.COLOR_BGR2RGB)

	# Estimate key-points and scores.
	pose = demo.process(im_name, image)

	# Access key-points and store them.
	key_points = pose['result'][0]['keypoints'].detach().cpu().numpy()
	score = pose['result'][0]['kp_score'].detach().cpu().numpy()
	score_binary = score > 0.4
	key_points = key_points * score_binary
	key_points = key_points.astype(np.int32)

	if save_keypts: np.savez(os.path.join(outputpath, os.path.basename(im_name).split('.')[0]+'.npz'), key_points = key_points)

	img = cv2.cvtColor(cv2.imread(im_name), cv2.COLOR_BGR2RGB)
	img = demo.vis(img, pose)   # visulize the pose result
	if save_img: cv2.imwrite(os.path.join(outputpath, 'vis', os.path.basename(im_name)), cv2.cvtColor(img, cv2.COLOR_BGR2RGB))

	# if you want to vis the img:
	# import matplotlib.pyplot as plt
	# plt.imshow(img)
	# plt.show()

	# write the result to json:
	# if save_json:
	# 	result = [pose]
	# 	demo.writeJson(result, outputpath, form=args.format, for_eval=args.eval)

if __name__ == "__main__":
	save_img = True
	save_keypts = True
	save_json = False
	evaluate(save_img, save_keypts, save_json)


text:
"""This module contains simple helper functions """

from __future__ import print_function
import torch
import numpy as np
from PIL import Image
import os

def compute_psnr(img1, img2, max_pixel_value=255.):
    """ Image Peak Signal to Noise Ratio
        Reference: https://dsp.stackexchange.com/questions/38065/peak-signal-to-noise-ratio-psnr-in-python-for-an-image
    """

    mse = torch.mean( (img1 - img2) ** 2)
    if mse:
        psnr = 10. * torch.log10((max_pixel_value**2)/mse)

        return psnr.item()
    else:
        return 1000.

def tensor2im(input_image, imtype=np.uint8):
    """"Converts a Tensor array into a numpy image array.

    Parameters:
        input_image (tensor) --  the input image tensor array of shape (1, 3/1, h, w)
        imtype (type)        --  the desired type of the converted numpy array
    """
    if not isinstance(input_image, np.ndarray):
        if isinstance(input_image, torch.Tensor):  # get the data from a variable
            image_tensor = input_image.data
        else:
            return input_image
        image_numpy = image_tensor[0].cpu().float().numpy()  # convert it into a numpy array

        if image_numpy.shape[0] == 1:  # grayscale to RGB
            image_numpy = np.tile(image_numpy, (3, 1, 1))
        image_numpy = (np.transpose(image_numpy, (1, 2, 0)) * 255.0)  # post-processing: transpose and scaling
    else:  # if it is a numpy array, do nothing
        image_numpy = input_image
    return image_numpy.astype(imtype)

def save_image(image_numpy, image_path, aspect_ratio=1.0):
    """Save a numpy image to the disk

    Parameters:
        image_numpy (numpy array) -- input numpy array
        image_path (str)          -- the path of the image
    """

    image_pil = Image.fromarray(image_numpy)
    h, w, _ = image_numpy.shape

    if aspect_ratio > 1.0:
        image_pil = image_pil.resize((h, int(w * aspect_ratio)), Image.BICUBIC)
    if aspect_ratio < 1.0:
        image_pil = image_pil.resize((int(h / aspect_ratio), w), Image.BICUBIC)
    image_pil.save(image_path)

def mkdirs(paths):
    """create empty directories if they don't exist

    Parameters:
        paths (str list) -- a list of directory paths
    """
    if isinstance(paths, list) and not isinstance(paths, str):
        for path in paths:
            mkdir(path)
    else:
        mkdir(paths)

def mkdir(path):
    """create a single empty directory if it didn't exist

    Parameters:
        path (str) -- a single directory path
    """
    if not os.path.exists(path):
        os.makedirs(path)

max_stars_repo_name: pkmandke/low-light-enhance

max_issues_repo_name: pkmandke/low-light-enhance

max_forks_repo_name: pkmandke/low-light-enhance


text:
from pathlib import Path
from typing import Union
from warnings import warn
from tempfile import NamedTemporaryFile
from subprocess import Popen, PIPE

from numpy import nan
from pandas import read_csv
from pandas import DataFrame, concat, Series
from pandas.errors import EmptyDataError
from rpy2.robjects import r
from rpy2.robjects.packages import importr

from helpers.temp import create_tmp_dir
from data_sources.molecular_signatures_db import MolecularSignaturesDatabase
from enhanced_multiprocessing.cache_manager import multiprocess_cache_manager

from ..models import Profile
from ..models.with_controls import ExpressionWithControls
from . import scoring_function
from .gsea import combine_gsea_results

db = MolecularSignaturesDatabase()
GSVA_CACHE = None
multiprocess_cache_manager.add_cache(globals(), 'GSVA_CACHE', 'dict')

gsva_tmp_dir = create_tmp_dir('gsva')
vanilla_R = ['R', '--vanilla', '--quiet']

def gsva(
    expression: Union[ExpressionWithControls, Profile], gene_sets_path: str, method: str = 'gsva',
    single_sample=False, permutations=1000, mx_diff=True, cores=1, _cache=True, limit_to_gene_sets=False,
    verbose=False
):
    """
    Excerpt from GSVA documentation:
        An important argument of the gsva() function is the flag mx.diff which is set to TRUE by default.

        Under this default setting, GSVA enrichment scores are calculated using Equation 5, and therefore, are
        more amenable by analysis techniques that assume the data to be normally distributed.  When setting
        mx.diff=FALSE , then Equation 4 is employed, calculating enrichment in an analogous way to classical
        GSEA which typically provides a bimodal distribution of GSVA enrichment scores for each gene.
    """

    if not single_sample and permutations:
        raise warn('permutations are not supported when not single_sample')

    key = (expression.hashable, method, gene_sets_path)

    if key in GSVA_CACHE:
        return GSVA_CACHE[key]

    if single_sample:
        assert isinstance(expression, Profile)
        joined = DataFrame(
            concat([expression.top.up, expression.top.down]),
        )
        joined.columns = ['condition']
        joined['control'] = 0
        joined.index = joined.index.astype(str)

        expression_classes = Series(['case', 'control'])
        expression = joined

    else:
        joined = expression.joined

        joined = DataFrame(joined)
        joined.index = joined.index.astype(str)

        nulls = joined.isnull().any(axis=0).reset_index(drop=True)
        if nulls.any():
            print(f'Following columns contain nulls and will be skipped: {list(joined.columns[nulls])}')
        joined = joined[joined.columns[~nulls]]

        expression_classes = expression.classes.loc[~nulls.reset_index(drop=True)]
        expression = joined

    mx_diff = 'T' if mx_diff else 'F'
    procedure = 'gene_permutation' if single_sample else 'bayes'
    cwd = Path(__file__).parent

    with NamedTemporaryFile(mode='w', prefix=gsva_tmp_dir) as f_expression, NamedTemporaryFile(prefix=gsva_tmp_dir) as f_result:
        expression.to_csv(f_expression)
        script = f"""
        source("{cwd}/gsva.R")
        expression = read.csv('{f_expression.name}', row.names=1)
        expression_classes = c{tuple(expression_classes)}
        gene_sets = readRDS('{gene_sets_path}')

        result = gsva.with_probabilities(
            expression, expression_classes, gene_sets, '{procedure}',
            method = '{method}', mx.diff={mx_diff}, include_control=F, cores={cores},
            limit_to_gene_sets={'c' + str(tuple(limit_to_gene_sets)) if limit_to_gene_sets is not False else 'F'}, progress=F
            {', permutations = ' + str(permutations) if procedure == 'permutations' else ''}
        )
        write.csv(result, '{f_result.name}')
        """

        process = Popen(vanilla_R, stdin=PIPE, stdout=PIPE, stderr=PIPE)
        r = process._stdin_write(script.encode())
        if verbose:
            from helpers.streams import handle_streams
            from signature_scoring.evaluation import display

            handlers = {'out': display, 'err': warn}
            handle_streams(process, handlers)
        else:
            process.wait()
        try:
            result = read_csv(f_result.name, index_col=0)
        except EmptyDataError:
            result = DataFrame()

    if _cache:
        GSVA_CACHE[key] = result
    return result

def create_gsva_scorer(
    gene_sets='c2.cp.kegg', id_type='entrez', grouping='by_substance',
    q_value_cutoff=0.1, na_action='fill_0', method='gsva', single_sample=False,
    permutations=None, mx_diff=True, custom_multiprocessing=False
):

    # as long as dummy controls are not included (include_control=F) there is no problem, otherwise:
    # if single_sample and method == 'plage':
    #    warn('PLAGE is not suitable for single sample testing')

    gmt_path = db.resolve(gene_sets, id_type)
    gsea_base = importr('GSEABase')
    gene_sets_r = gsea_base.getGmt(gmt_path)

    # transform to named list from GeneSetCollection class object
    gene_sets_r = gsea_base.geneIds(gene_sets_r)

    gene_sets_file = NamedTemporaryFile(delete=False, prefix=gsva_tmp_dir)
    r['saveRDS'](gene_sets_r, file=gene_sets_file.name)

    input = Profile if single_sample else ExpressionWithControls

    def gsva_score(disease: input, compound: input, cores=1):
        if not custom_multiprocessing:
            multiprocess_cache_manager.respawn_cache_if_needed()

        disease_gene_sets = gsva(
            disease, gene_sets_path=gene_sets_file.name, method=method, single_sample=single_sample,
            permutations=permutations, mx_diff=mx_diff, cores=cores
        )

        disease_gene_sets.drop(disease_gene_sets[disease_gene_sets['fdr_q-val'] > q_value_cutoff].index, inplace=True)

        assert len(disease_gene_sets.index)
        signature_gene_sets = gsva(
            compound, gene_sets_path=gene_sets_file.name, method=method, single_sample=single_sample,
            permutations=permutations, mx_diff=mx_diff, _cache=False, cores=cores,
            limit_to_gene_sets=list(disease_gene_sets.index)
        )
        if signature_gene_sets.empty:
            return nan

        joined = combine_gsea_results(disease_gene_sets, signature_gene_sets, na_action)
        assert not (set(signature_gene_sets.index) - set(disease_gene_sets.index))
        return joined.score.mean()

    gsva_score.__name__ = (
        method +
        f'_{permutations}' +
        f'_mx_diff:{mx_diff}' +
        ('_single_sample' if single_sample else '')
    )

    return scoring_function(gsva_score, input=input, grouping=grouping, custom_multiprocessing=custom_multiprocessing)

reason: from numpy


text:
from .data import load_cam_output, uniform_stratify_data
import unittest
import pandas as pd
import numpy as np

class TestData(unittest.TestCase):
    def setUp(self):
        self.bad_path = "/blah/de/blah"
        return

    def test_load_cam_output(self):
        with self.assertRaises(FileNotFoundError):
            load_cam_output(self.bad_path)

    def test_uniform_stratify_data(self):
        samples = 50000
        random_seed = 1232
        np.random.seed(random_seed)
        output_labels = pd.DataFrame({"two_labels": np.random.randint(0, 2, samples),
                                      "three_labels": np.random.randint(-1, 2, samples)})
        scaled_output_values = pd.DataFrame({"two_labels": np.zeros(samples),
                                             "three_labels": np.zeros(samples)})
        category_size = 500
        num_bins = 10
        bins = {}
        for col in output_labels.columns:
            nonzero_samples = output_labels[col] != 0
            scaled_output_values.loc[nonzero_samples, col] = np.random.normal(size=np.count_nonzero(nonzero_samples))
            bins[col] = {}
            labels = np.unique(output_labels[col])
            for label in labels:
                if label != 0:
                    label_vals = scaled_output_values.loc[output_labels[col] == label, col].values
                    print(label_vals.min(), label_vals.max())
                    bins[col][label] = np.linspace(label_vals.min(), label_vals.max() + 0.1, num_bins)
        sampling_indices = uniform_stratify_data(output_labels, scaled_output_values, category_size, bins)
        for out_label in output_labels.columns:
            for label in sampling_indices[out_label].keys():
                self.assertEqual(sampling_indices[out_label][label].size, category_size)

max_stars_repo_name: NCAR/mlmicrophysics

max_issues_repo_name: NCAR/mlmicrophysics

max_forks_repo_name: NCAR/mlmicrophysics