text:
# -*- coding: utf-8 -*-
"""
Created on Mon Aug 31 01:04:06 2015

@author: diego

"""

import numpy as np
import pylab as plt
import glob,subprocess,os
import del_dot_xi as ddxi
import pyNRO_1_mode as pyNRO
import scipy.integrate as scint
import legendre_interp as lint
import seaborn as sns
from timeit import default_timer as timer

from main_modules import find_vel,find_index,find_sigma,find_name,index_by_name,emode,phaser,phaser2,find_dr_ofmax

vis_path = os.getcwd()

global homedir
if (os.path.isfile(vis_path+"/lachesis"))==False:
    homedir = "/home/castaned/Documents/"
else:
    homedir = "/home/castaned/"

def run_visc_pert(model,vel,mode,par,sigma,reese,force_f,minL,phase=0.,ampl=1.,tlmax=0.,tcut=False):
    # Info for del dot xi calculation:------------------
    # NRO Mode filename:
    #modefname="MODE1"
    global kind

    norm_f = True
    ssc = 0.02
    if model[0]=="2p5": ssc = 0.01
    if kind=="g":
        scale = ssc
    else:
        scale = ssc/(sigma**2)
    scale =0.001*drtempmax
    depth = 10 #radial zones down from the surface
    forcenro = False
    #---------------------------------------------------

    #clic_run = True

    v = find_vel(model,vel)
    folder = homedir+"ROTORCmodels/"+par+"/M"+model[0]+"_V"+v+"/"

    rotorc_f = homedir+"From_Bob/Delta_Scuti_2010/"+model[0]+"Msun/"+model[0]+"Msun_V"+v+"/"

    bob_bin = homedir+"From_Bob/clotho_disc10_bin/"

    static_m = homedir+"ROTORCmodels/visibilities/"

    #check if visibility file from pulset exits for the selected model:

    #v_file = glob.glob(rotorc_f+"visibility_file")

    ###### FULL MODE FILE GENERATION:

    temp_freqs = np.genfromtxt(folder+"temp_freqs")
    tfreq = temp_freqs[mode-1]
    if force_f==True:
        tfreq = sigma
        print sigma
    #tfreq = 1.59692
    #print pyNRO.run_nro(tfreq,folder,model[0],v,par,mode)

    where =static_m+model[0]+'Msun/V'+v+"/MODE_"+par+"_"+str(mode)

    if not os.path.exists(static_m+model[0]+'Msun/'+'V'+v):
        os.makedirs(static_m+model[0]+'Msun/'+'V'+v)

    if not os.path.exists(static_m+model[0]+'Msun/'+'V'+v+"/MODE_"+par+"_"+str(mode)):
        os.makedirs(where)

    if os.path.isfile(where+'/MODE_'+par+'_'+str(mode))==False or forcenro==True:
        #check if visibility file from pulset exits for the selected model:

        #v_file = glob.glob(rotorc_f+"visibility_file")
        if True:
            os.chdir(rotorc_f)
            r_mod_name = glob.glob("*_ZAMS")
            subprocess.call(["cp",r_mod_name[0],"orcmod_pulset"])
            subprocess.call([bob_bin+"pulsetnonadb.exe"])
            print "Generated visibility_file in "+rotorc_f
            subprocess.call([bob_bin+"pulset_gammas.exe"]) #Generates Cmod with gammas
            subprocess.call(["cp","Cmod",folder+"Dmod_"+model[0]+"M_V"+v])
            print "Generated new Dmod"
            #print(glob.glob("*_ZAMS"))
            os.chdir(vis_path)

        print "Dmod with GAMMAS generated!"

        nro_stat = 1
        idx_iter = 0
        while nro_stat==1:
        #RUN NRO!
            nro_stat = pyNRO.run_nro(tfreq-idx_iter*1e-4,folder,model[0],v,par,mode,homedir)
            idx_iter += 1

        subprocess.call(['mv',folder+'MODE_'+par+'_'+str(mode),where+'/MODE_'+par+'_'+str(mode)])
        print "Mode file generation complete!"
    else:
        print "Mode file found! Not running NRO..."

    ###### del_DOT_xi>
    #modeloc = static_m+model[0]+'Msun/V'+vels[vel]+"/"
    modeloc = where+"/"
    modefname = 'MODE_'+par+'_'+str(mode)
    #modefname = 'MODE13'
    global s_model
    s_model = np.genfromtxt(glob.glob(static_m+model[0]+'Msun/'+model[0]+'Msun_V'+v+"*")[0])

    global xi_r_rot,xi_t_rot,dt_t_rot,zg_rot
    global xi_r,xi_t,dt_t,zg,r,zp,cs,xi_dot_g
    global xi_r_n,xi_t_n,dt_t_n,zg_n

    xi_r,xi_t,dt_t,zg,r,zp,sig,cs,xi_dot_g = ddxi.calcdeldotxi(par,model,vel,modeloc,modefname)

    xi_r_n,xi_t_n,dt_t_n,zg_n,dr_n = ddxi.norm_and_scale(xi_r,xi_t,dt_t,r,zg,norm_f,scale,depth,reese,sig,par)

    if phase!=0:
        #dt_t_n = phaser(dt_t_n,phase,1.,tlmax)
        #xi_r_n = phaser(xi_r_n,0.,1.,tlmax)
        dt_t_n,psi_T,psi_L,mx = phaser2(dt_t_n,phase,np.max(np.abs(dt_t[-depth])),np.max(np.abs(xi_r[-depth])))
        xi_r_n = phaser(xi_r_n,0.,1.,mx+np.pi)
        print "Theta L_max-----------> ",np.rad2deg(mx),"deg"
        np.savetxt(modeloc+"phases.txt",[psi_T,psi_L,np.rad2deg(mx+np.pi)],header="psi_T[deg],psi_L[deg],theta_max_L[deg]")

    a_r = scint.trapz(dt_t_n[-1,:])

    inv_f = 1.
    if a_r<0:
        print "-T area"
        #xi_r_n *= -1.
        #xi_t_n *= -1.
        #dt_t_n *= -1.
        #zg_n *= -1.
        inv_f = -1.

    if minL==True:
        inv_f *= -1.

    #xi_r_rot,xi_t_rot,dt_t_rot,zg_rot = ddxi.to_rotorc(xi_r_n,xi_t_n,dt_t_n,zg_n)

    if par == "EVEN":
        xi_r_fine = lint.leg_interp(xi_r_n[:,:],8,"EVEN")
        xi_t_fine = lint.leg_interp(xi_t_n[:,:],8,"EVEN")
        dt_t_fine = lint.leg_interp(dt_t_n[:,:],8,"EVEN")
    else:
        xi_r_fine = lint.leg_interp(xi_r_n[:,:],8,"OE")
        xi_t_fine = lint.leg_interp(xi_t_n[:,:],8,"OE")
        dt_t_fine = lint.leg_interp(dt_t_n[:,:],8,"OE")

    global rs_n,rs_rot
    rs_n = np.empty(xi_r_fine.shape)
    for d in range(depth):
        rs_n[d,:] = np.interp(np.linspace(0,90,100),np.linspace(10,80,8),r[-d-1,:])

#    dr_n = xi_r_fine*rs_n
#    for i in range(len(dr_n[:,0])):
#        imax = np.argmax(np.abs(dr_n[i,:]))
#        dr_n[i,:] = dr_n[i,:]/(dr_n[i,imax]/2.28534026)

    p_angles = np.empty(np.shape(xi_r_n))
    p_angles_fine = np.empty(np.shape(xi_r_fine))
    rot_ang = np.arange(4.5,90.,9)
    xi_r_rot = np.empty((depth,len(rot_ang)))
    xi_t_rot = np.empty((depth,len(rot_ang)))
    dt_t_rot = np.empty((depth,len(rot_ang)))
    rs_rot = np.empty((depth,len(rot_ang)))
    for d in range(depth):
        p_angles[d,:] = (np.linspace(10,80,8))*(1.+xi_t_n[d,:])
        p_angles_fine[d,:] = (np.linspace(0,90,100))*(1.+xi_t_fine[d,:])
        xi_r_rot[d,:] = np.interp(rot_ang,p_angles_fine[d,:],xi_r_fine[d,:])
        xi_t_rot[d,:] = np.interp(rot_ang,p_angles_fine[d,:],xi_t_fine[d,:])
        dt_t_rot[d,:] = np.interp(rot_ang,p_angles_fine[d,:],dt_t_fine[d,:])
        rs_rot[d,:] = np.interp(rot_ang,p_angles_fine[d,:],rs_n[d,:])

#
#    #plt.plot(p_angles_fine[2,:],xi_t_fine[2,:])
#    if dr_n[0,-1]>0:
#        ivrt = 1.
#    else:
#        ivrt = -1.
#    vr = ivrt*dr_n[0,:]
#    #vr = ivrt*dt_t_fine[0,:]
#    #lblb = find_vel(model,vel)+" km s$^{-1}$ "+find_name(model,vel,par,mode)
#    lblb = find_vel(model,vel)+" km s$^{-1}$ "
#    #lblb = r"$\ell$ = "+find_name(model,vel,par,mode).split()[0]
#    etsy = "-"
#    if find_name(model,vel,par,mode).split()[0] == "2": etsy = "-"
#    if par == "EVEN":
#        #plt.plot(p_angles_fine[2,:],ivrt*lint.leg_interp(vr[:,:],8,"EVEN")[0,:],label=lblb)
#        plt.plot(p_angles_fine[2,:],vr,ls=etsy,label=lblb)
#
#    else:
#        #plt.plot(p_angles_fine[2,:],ivrt*lint.leg_interp(vr[:,:],8,"OE")[0,:],label=lblb)
#        plt.plot(p_angles_fine[2,:],vr,ls=etsy,label=lblb)
#    #plt.plot(p_angles_fine[2,:],lint.leg_interp(xi_r[-depth::,:],8,"OE")[2,:],label=find_name(model,vel,par,mode))
#    #plt.plot(p_angles[2,:],dt_t[-depth+2,:],"o",mfc="white",mec="k",mew=1)
#    #plt.grid()
#    plt.xlim(0,90)
#    plt.yticks()
#    plt.xlabel("Colatitude [deg]")
#    plt.ylabel(r"$\delta$R [R$\odot$]")
#    #plt.ylabel(r"$\delta$T/T")
#    plt.legend(loc="best")
#    #plt.title("Mode:" + find_name(model,vel,par,mode))
#    #ax.yaxis.set_major_formatter(majorFormatter)
#    plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
#    #plt.title(r"M="+model[0]+"M$_{\odot}$, V="+v+"km s$^{-1}$")
#    plt.title(r"M="+model[0]+"M$_{\odot}$")

    print "Generating fine clic model..."
    pmodels_fine = gen_fine_clic(p_angles_fine,xi_r_fine,xi_t_fine,dt_t_fine,rs_n,dr_n,s_model,model,depth,inv_f,par,tcut)

    ############ find the perturbed models:
    global omega
    G = 6.67259e-8
    mass = (float((model[0]).replace("p",".")))*1.99e33
    theta = np.deg2rad(s_model[:,0])
    radius = s_model[:,1]*6.96e10
    vrot = s_model[:,4]*1e5
    omega = vrot/(radius*np.sin(theta))
    omega = np.average(omega[1:-1])

    #g_r = (G*mass/(radius**2))-((omega**2)*(radius*np.sin(theta)))*np.sin(theta)
    #g_theta = ((omega**2)*(radius*np.sin(theta)))*np.cos(theta)

    #geff = (g_r**2 + g_theta**2)**0.5

    pert_r = np.empty((len(s_model[:,0]),depth))
    pert_t = np.empty((len(s_model[:,0]),depth))
    pmodels = []
    for i in range(depth):
        pert_r[:,i] = s_model[:,1]*(1.+inv_f*xi_r_rot[-i,:])
        pert_t[:,i] = s_model[:,2]*(1.+inv_f*dt_t_rot[-i,:])
        g_pert = (G*mass/((pert_r[:,i]*6.96e10)**2))-((omega**2)*((pert_r[:,i]*6.96e10)*np.sin(theta)))*np.sin(theta)
        for j in range(len(g_pert)):
            if np.log10(g_pert[j])>4.33: g_pert[j]=10**(4.33)
        tmodel = 1.*s_model
        tmodel[:,1] = pert_r[:,i]
        tmodel[:,2] = pert_t[:,i]
        tmodel[:,3] = np.log10(g_pert)
        pmodels.append(tmodel)

    if par=="ODD":
        pert_r = np.empty((2*len(s_model[:,0]),depth))
        pert_t = np.empty((2*len(s_model[:,0]),depth))

        pmodels = []
        odd_angles = np.arange(4.5,180.,9)
        for i in range(depth):
            pert_r[0:10,i] = s_model[:,1]+inv_f*xi_r_rot[-i,:]*rs_rot[-i,:]
            pert_t[0:10,i] = s_model[:,2]*(1.+inv_f*dt_t_rot[-i,:])
            pert_r[10:20,i] = s_model[::-1,1]*(1.-inv_f*xi_r_rot[-i,::-1])
            pert_t[10:20,i] = s_model[::-1,2]*(1.-inv_f*dt_t_rot[-i,::-1])
            g_pert = (G*mass/((pert_r[:,i]*6.96e10)**2))-((omega**2)*((pert_r[:,i]*6.96e10)*np.sin(odd_angles)))*np.sin(odd_angles)
            for j in range(len(g_pert)):
                if np.log10(g_pert[j])>4.33: g_pert[j]=10**(4.33)
            tmodel = np.empty((20,5))
            tmodel[:,0] = odd_angles[:]
            tmodel[:,1] = pert_r[:,i]
            tmodel[:,2] = pert_t[:,i]
            #tmodel[0:10,3] = s_model[:,3]
            #tmodel[10:20,3] = s_model[::-1,3]
            tmodel[:,3] = np.log10(g_pert)
            tmodel[0:10,4] = s_model[:,4]
            tmodel[10:20,4] = s_model[::-1,4]

            pmodels.append(tmodel)

#    plt.plot(pmodels[-1][:,0],s_model[:,1],"--",color="0.5")
#
#    global old_modes
#    import seaborn as sns
#    sns.set(style="white",rc={"figure.figsize": (8, 8),'axes.labelsize': 16,
#                              'ytick.labelsize': 12,'xtick.labelsize': 12,
#                              'legend.fontsize': 16,'axes.titlesize':18,'font.size':14})
#    #plt.plot(pmodels[2][0:-1,0],np.diff(pmodels[2][:,1])/(pmodels[2][1,0]-pmodels[2][0,0]),label=r"$\ell=$"+find_name(model,vel,par,mode).strip().split()[0])
#    #plt.plot(pmodels[2][:,0],pmodels[2][:,1],label="old way")
#    #plt.plot(pmodels_fine[2][0:-1,0],np.diff(pmodels_fine[2][:,1])/(pmodels_fine[2][1,0]-pmodels_fine[2][0,0]))
#    plt.plot(pmodels_fine[2][:,0],pmodels_fine[2][:,2],"-",lw=1.5,label=r"$\ell=$"+find_name(model,vel,par,mode).strip().split()[0])
#    #plt.plot(s_model[:,0],s_model[:,1],label="Static")
#    #plt.vlines(90,9100,9300)
#    o_lims = plt.ylim()
#    #plt.vlines(90,min(pmodels_fine[2][:,2])-100,max(pmodels_fine[2][:,2])+100)
#    plt.grid()
#    plt.xlim(0,180)
#    #plt.ylim(o_lims[0],o_lims[1])
#    plt.xlabel("Colatitude [deg]")
#    plt.ylabel("Radius [M$_{\odot}$]")
#    #plt.ylabel("Temperature [K]")
#    plt.legend(loc="best")
#    #plt.title("Mode: " + find_name(model,vel,par,mode))
    #plt.title(r"Perturbed T$_{\mathrm{eff}}$ - M="+model[0]+"M$_{\odot}$, V="+v+"km s$^{-1}$")

    return modeloc,pmodels,pmodels_fine

def gen_fine_clic(p_angles_fine,xi_r,xi_t,dt_t,rs,dr_n,s_model,model,depth,inv_f,par,tcut):
    clic_theta_grid = 200
    clic_theta = np.linspace(0,180,clic_theta_grid)
    global clic_xi_r,clic_xi_t,clic_dt_t

    clic_xi_r = np.zeros((depth,clic_theta_grid))
    clic_xi_t = np.zeros((depth,clic_theta_grid))
    clic_dt_t = np.zeros((depth,clic_theta_grid))
    clic_rs = np.zeros((depth,clic_theta_grid))
    if tcut:
        cut=True
        print "--- CUT ---"
    else:
        cut=False
    for d in range(depth):
        if cut:
            xi_r[d,0:12] = 1.*xi_r[d,12]
            xi_t[d,0:12] = 1.*xi_t[d,12]
            dt_t[d,0:12] = 1.*dt_t[d,12]
            rs[d,0:12] = 1.*rs[d,12]
            if par=="EVEN":
                xi_r[d,-9::] = 1.*xi_r[d,-10]
                xi_t[d,-9::] = 1.*xi_t[d,-10]
                dt_t[d,-9::] = 1.*dt_t[d,-10]
                rs[d,-9::] = 1.*rs[d,-10]
        clic_xi_r[d,0:100] = np.interp(clic_theta[0:100],np.linspace(0,90,100),xi_r[d,:])
        clic_xi_t[d,0:100] = np.interp(clic_theta[0:100],np.linspace(0,90,100),xi_t[d,:])
        clic_dt_t[d,0:100] = np.interp(clic_theta[0:100],np.linspace(0,90,100),dt_t[d,:])
        clic_rs[d,0:100] = np.interp(clic_theta[0:100],np.linspace(0,90,100),rs[d,:])

    for i in range(depth):
        if par=="EVEN":
            clic_xi_r[i,100::] = clic_xi_r[i,::-1][100::]
            clic_xi_t[i,100::] = clic_xi_t[i,::-1][100::]
            clic_dt_t[i,100::] = clic_dt_t[i,::-1][100::]
        else:
            clic_xi_r[i,100::] = -1.*clic_xi_r[i,::-1][100::]
            clic_xi_t[i,100::] = -1.*clic_xi_t[i,::-1][100::]
            clic_dt_t[i,100::] = -1.*clic_dt_t[i,::-1][100::]

        clic_rs[i,100::] = clic_rs[i,::-1][100::]

    #plt.plot(clic_theta[0:200],clic_rs[-2,0:200])

    ############ find the perturbed models:
    # Take static ROTORC model to CLIC grid:
    global s_model_c
    s_model_c = np.zeros((clic_theta_grid,len(s_model[0,:])))
    s_model_c[:,0] = clic_theta
    for i in range(1,len(s_model[0,:])):
        s_model_c[0:100,i] = np.interp(clic_theta[0:100],s_model[:,0],s_model[:,i])
        s_model_c[100::,i] = s_model_c[::-1,i][100::]

    G = 6.67259e-8
    mass = (float((model[0]).replace("p",".")))*1.99e33
    theta = np.deg2rad(s_model_c[:,0])
    radius = s_model_c[:,1]*6.96e10
    vrot = s_model_c[:,4]*1e5
    omega = vrot/(radius*np.sin(theta))
    omega = np.average(omega[50:-50])

    a_r = scint.trapz(clic_dt_t[-1,:])

#    inv_f = 1.
#    if a_r<0:
#        print "fine -T area"
#        inv_f = -1.

    global pmodels_fine
    pert_r = np.empty((len(s_model_c[:,0]),depth))
    pert_t = np.empty((len(s_model_c[:,0]),depth))
    pmodels_fine = []
    for i in range(depth):
        #pert_r[:,i] = s_model_c[:,1]*(1.+inv_f*clic_xi_r[-i,:])
        pert_r[:,i] = s_model_c[:,1]+inv_f*clic_rs[-i,:]*clic_xi_r[-i,:]
        pert_t[:,i] = s_model_c[:,2]*(1.+inv_f*clic_dt_t[-i,:])
        for j in range(len(pert_t[:,i])):
            if pert_t[j,i]<7500.: pert_t[j,i]=7500.
        g_pert = (G*mass/((pert_r[:,i]*6.96e10)**2))-((omega**2)*((pert_r[:,i]*6.96e10)*np.sin(theta)))*np.sin(theta)
        for j in range(len(g_pert)):
            if np.log10(g_pert[j])>4.33: g_pert[j]=10**(4.33)
        tmodel = 1.*s_model_c
        """
        pert_r[0:11,:] = pert_r[11,:]
        pert_r[190:200,:] = pert_r[189,:]

        pert_t[0:11,:] = pert_t[11,:]
        pert_t[190:200,:] = pert_t[189,:]
        """
        tmodel[:,1] = pert_r[:,i]
        tmodel[:,2] = pert_t[:,i]
        tmodel[:,3] = np.log10(g_pert)
        pmodels_fine.append(tmodel)

    #plt.plot(clic_dt_t[0,:])
    return pmodels_fine

def make_contour(r,var,model,vel,par,mode,x):
    global new_r,new_zp
    import seaborn as sns
    sns.set(style="white",rc={"figure.figsize": (8, 8),'axes.labelsize': 16,
                              'ytick.labelsize': 12,'xtick.labelsize': 12,'axes.titlesize':18})

    t_var = np.empty(var[x::,:].shape)
    t_r = np.empty(t_var.shape)
    for i in range(len(r[0,:])):
        t_var[:,i] = var[x::,i]
        t_r[:,i] = np.linspace(r[0,i],r[-1,i],len(t_var[:,i]))

    var = 1.*t_var
    r = 1.*t_r

    new_r = lint.leg_interp(r[0:-2,:],8,"EVEN")
    new_zp = lint.leg_interp(var[0:-2,:],8,"EVEN")
    #levels = np.linspace(np.min(new_zp),np.max(new_zp), 40)

    #####REMOVE 0-10deg at the pole:
    new_r = new_r[:,11::]
    new_zp = new_zp[:,11::]
    theta = np.linspace(0,np.deg2rad(90),89)
    #########################################

    #theta = np.linspace(0,np.deg2rad(90),100)
    newt,n_r = np.meshgrid(theta,new_r[:,0])

#    np.savetxt("/home/castaned/x.dat",(new_r*np.sin(newt)))
#    np.savetxt("/home/castaned/y.dat",(new_r*np.cos(newt)))
#    np.savetxt("/home/castaned/z.dat",(new_zp))

    plt.contourf((new_r*np.sin(newt)),(new_r*np.cos(newt)),new_zp, 100, cmap=plt.cm.gnuplot,vmax=np.max(new_zp), vmin=np.min(new_zp))
    plt.contourf((new_r*np.sin(newt)),-(new_r*np.cos(newt)),new_zp, 100, cmap=plt.cm.gnuplot,vmax=np.max(new_zp), vmin=np.min(new_zp))
    plt.contourf(-(new_r*np.sin(newt)),-(new_r*np.cos(newt)),new_zp, 100, cmap=plt.cm.gnuplot,vmax=np.max(new_zp), vmin=np.min(new_zp))
    plt.contourf(-(new_r*np.sin(newt)),(new_r*np.cos(newt)),new_zp, 100, cmap=plt.cm.gnuplot,vmax=np.max(new_zp), vmin=np.min(new_zp))
    #CSl = plt.contour((new_r*np.sin(newt)),(new_r*np.cos(newt)),new_zp, 20, colors="k")
    #CS = plt.contourf((new_r*np.cos(newt)), (new_r*np.sin(newt)), new_zp, cmap=plt.cm.Spectral,levels=levels)
    #plt.axes().set_aspect('equal')
    #plt.xlim(np.ceil(r[-1,-1]))
    plt.ylim(-plt.xlim()[1],plt.xlim()[1])
    #plt.xlabel("Radius [R$_{\odot}$]")
    #plt.ylabel("Radius [R$_{\odot}$]")
    plt.axis('off')
    try:
        m = (find_name(model,vel,par,mode)).strip()
    except:
        m = "freq_%.5f" % find_sigma(model,vel,par,mode)
    v = find_vel(model,vel)

    plt.title("M"+model[0]+"_V"+v+" - mode: "+ m)
    return

def save_pert(modeloc,model,vel,par,mode,name):
    global xi_r_rot,xi_t_rot,dt_t_rot,zg_rot
    global xi_r,xi_t,dt_t,zg,r
    global xi_r_n,xi_t_n,dt_t_n,zg_n

    print os.getcwd()
    f = open(modeloc+"M"+model[0]+"_V"+str(vel)+"_"+name+"_perturbations","w")
    f.write("M"+model[0]+", V="+str(vel)+", Mode freq = %.5f\n" % find_sigma(model,vel,par,mode))
    f.write("n_angle, r, xi_r, xi_t, dT/T\n")
    for i in range(len(r[:,0])-10,len(r[:,0])):
        for j in range(len(xi_r[-1,:])):
            #print "%i %8.5f %8.5f %8.5f %8.5f\n"%(j+1,r[i,j],xi_r[i,j],xi_t[i,j],dt_t[i,j])
            f.write("%i %8.5f %8.5f %8.5f %8.5f\n"%(j+1,r[i,j],xi_r[i,j],xi_t[i,j],dt_t[i,j]))

    f.close()
    return

def list_gen(ells,nmin,nmax,tpe):
    ls = []
    for j in ells:
        for i in range(nmin,nmax+1):
            ls.append(str(j)+" "+tpe+str(i))

    return ls
"""

------------------------------------------------------
------------------------------------------------------

"""
def run_pipeline(modes,sub_dir="",minL=False,mm="1p875",vv=0,tph=0,drscale=None,tcut=False):
    global start,model,vel,old_modes,modeloc,pmodels,pmodels_fine,modes_not_found,mds
    #MODE info:

    model = [mm]
    vel = vv #index, not velocity!
    #modes = list_gen([0,1,2,3,4,5,6,7,8,9,10],5,5,"p")
    #list_gen([0,1,2,3],1,1,"p")
    #modes = list_gen([2],8,8,"p")
    #["1 p1","1 p2","1 p3","1 p4","1 p5", "1 p6", "1 p7", "1 p8", "1 p9"]
    #["0 1H","0 2H","0 3H","0 4H","0 5H","0 6H","0 7H","0 8H","0 9H","0 10H"]
    #par = emode(mm,vv,modes[0]).parity

    clic = False
    new_mags = False
    save_perturbation = False
    only_mags = False

    incl = list(np.arange(0,95,10))
    #incl = [90]

    reese = False
    force_f = False
    f_freq = 1.22

    phase = tph
    ampl = 1./2.2
    tlmax = 4.14 #radians (sigma*t of phased Lmax)

    #### contour plotting option:
    plot_contour = True
    excl = 20 #How many core zones to exclude
    ####

    freqs = [6.61043]

    dt_grid = 250

    old_modes = list(modes)
    mode_by_freq = False
    mode_by_name = False
    if type(modes[0]) == float:
        mode_by_freq = True

    if type(modes[0]) == str:
        mode_by_name = True

    if mode_by_freq==True:
            modes = find_index(freqs,model,vel,par)

    parity = []
    modesx = []
    if mode_by_name==True:
        for i in range(len(modes)):
            if int(modes[i].split()[0])%2==0:
                par = "EVEN"
            else:
                par = "ODD"

            parity.append(par)
            idx_tmp = index_by_name(model,vel,par,modes[i])
            print modes[i],"->",idx_tmp
            modesx.append(idx_tmp)

    start = timer()
    klm = 0
    rem = 0

    modes_not_found = []
    for mode in modesx:
        par = emode(mm,vv,mds[modesx.index(mode)]).parity
        if mode == 999:
            modes_not_found.append(old_modes[modes.index(mode)])
            continue

        rem += 1

        if mode_by_name == True:
            par = parity[klm]
            klm +=1

        if force_f==True:
            sigma=f_freq
        else:
            #sigma = find_sigma(model,vel,par,mode)
            sigma = emode(mm,vv,old_modes[modesx.index(mode)]).frequency

        #Find the perturbed models!

        modeloc,pmodels,pmodels_fine = run_visc_pert(model,vel,mode,par,sigma,reese,force_f,minL,phase,ampl,tlmax,tcut)
        modeloc += sub_dir
        if not os.path.exists(modeloc):
            os.makedirs(modeloc)

        print modeloc
        print mds[modesx.index(mode)]
        nrzone = len(dt_t[:,0])

        if plot_contour == True:
            make_contour(r,-zp*cs,model,vel,par,mode,excl)

        if save_perturbation==True:
            save_pert(modeloc,model,vel,par,mode,old_modes[modesx.index(mode)])

        for i in range(len(pmodels)):
            np.savetxt(modeloc+"model_MODE_"+str(mode)+"_r"+str(nrzone - i),pmodels[i],'%.16e')
            np.savetxt(modeloc+"fine_model_MODE_"+str(mode)+"_r"+str(nrzone - i),pmodels_fine[i],'%.16e')

        rzone = 2 # Radial zone from the surface
        rzone = 490
        print rzone

        if clic==True:
            #Ugly way of running pyCLIC:
            clic_folder = homedir+"pyCLIC/test/"

            if only_mags==False:
                subprocess.call(['mkdir',modeloc+'tmp'])
                os.chdir(modeloc+'tmp')
                tmpstr = "cp -r "+clic_folder+"* ."
                subprocess.call(tmpstr,shell=True)

                for i in incl:

                    if new_mags==True:
                        if (os.path.isfile(modeloc+"magnitudes"+"_i"+str(i)+"_MODE_"+str(mode)))==True:
                            subprocess.call(['rm',modeloc+"magnitudes"+"_i"+str(i)+"_MODE_"+str(mode)])
                        else:
                            print "No mag. file to delete!"
                        if (os.path.isfile(modeloc+"walraven_magnitudes"+"_i"+str(i)+"_MODE_"+str(mode)))==True:
                            subprocess.call(['rm',modeloc+"walraven_magnitudes"+"_i"+str(i)+"_MODE_"+str(mode)])
                        else:
                            print "No walraven_mag. file to delete!"

                    subprocess.call(['cp',modeloc+"model_MODE_"+str(mode)+"_r"+str(rzone),modeloc+'tmp'])
                    subprocess.call(['cp',modeloc+"fine_model_MODE_"+str(mode)+"_r"+str(rzone),modeloc+'tmp'])

                    if (os.path.isfile(modeloc+"magnitudes"+"_i"+str(i)+"_MODE_"+str(mode)))==False:
                        myfile = open(modeloc+"magnitudes"+"_i"+str(i)+"_MODE_"+str(mode), "w")
                        myfile.write('#%5s %5s %9s %11s %11s %11s %11s %11s\n' % ("incl","n_r","lum","u","b","v","r","i"))
                        myfile.close()
                        subprocess.call(['cp',modeloc+"magnitudes"+"_i"+str(i)+"_MODE_"+str(mode),modeloc+'tmp'])
                    else:
                        subprocess.call(['cp',modeloc+"magnitudes"+"_i"+str(i)+"_MODE_"+str(mode),modeloc+'tmp'])

                    if (os.path.isfile(modeloc+"walraven_magnitudes"+"_i"+str(i)+"_MODE_"+str(mode)))==False:
                        myfile = open(modeloc+"walraven_magnitudes"+"_i"+str(i)+"_MODE_"+str(mode), "w")
                        myfile.write('#%5s %5s %9s %11s %11s %11s %11s %11s\n' % ("incl","n_r","lum","w","u","l","b","v"))
                        myfile.close()
                        subprocess.call(['cp',modeloc+"walraven_magnitudes"+"_i"+str(i)+"_MODE_"+str(mode),modeloc+'tmp'])
                    else:
                        subprocess.call(['cp',modeloc+"walraven_magnitudes"+"_i"+str(i)+"_MODE_"+str(mode),modeloc+'tmp'])

                    os.chdir(modeloc+'tmp')
                    #sys.path.append(modeloc+'tmp')
                    #print sys.path
                    import pyCLIC_main as pyclic
                    #sys.path.append(vis_path)
                    os.chdir(vis_path)
                    import color_magnitudes as cmag
                    os.chdir(modeloc+'tmp')

                    pyclic.run_CLIC("model_MODE_"+str(mode)+"_r"+str(rzone),[i],False,3000.,7499.,10.0,par,dt_grid,fine_model=True,savexi=True)

                    cmag.calc_mags('outputflux_i'+str(i)+'.final',[i],mode,rzone)
                    cmag.calc_walraven('outputflux_i'+str(i)+'.final',[i],mode,rzone)

                    subprocess.call(['cp','outputflux_i'+str(i)+'.final',modeloc+"outputflux_i"+str(i)+"_MODE_"+str(mode)+"_r"+str(rzone)])
                    subprocess.call(['cp',"walraven_magnitudes"+"_i"+str(i)+"_MODE_"+str(mode),modeloc])
                    subprocess.call(['cp',"magnitudes"+"_i"+str(i)+"_MODE_"+str(mode),modeloc])

                    subprocess.call(['rm',"model_MODE_"+str(mode)+"_r"+str(rzone)])
                    subprocess.call(['rm','outputflux_i'+str(i)+'.final'])
                    subprocess.call(['rm',"magnitudes"+"_i"+str(i)+"_MODE_"+str(mode)])
                    subprocess.call(['rm',"walraven_magnitudes"+"_i"+str(i)+"_MODE_"+str(mode)])

                    os.chdir(vis_path)

                subprocess.call(['rm',"-r",modeloc+'tmp'])

            else:
                os.chdir(vis_path)
                import color_magnitudes as cmag
                os.chdir(modeloc)
                for i in incl:

                    if new_mags==True:
                        if (os.path.isfile(modeloc+"magnitudes"+"_i"+str(i)+"_MODE_"+str(mode)))==True:
                            subprocess.call(['rm',modeloc+"magnitudes"+"_i"+str(i)+"_MODE_"+str(mode)])
                        else:
                            print "No mag. file to delete!"
                        if (os.path.isfile(modeloc+"walraven_magnitudes"+"_i"+str(i)+"_MODE_"+str(mode)))==True:
                            subprocess.call(['rm',modeloc+"walraven_magnitudes"+"_i"+str(i)+"_MODE_"+str(mode)])
                        else:
                            print "No walraven_mag. file to delete!"

                    if (os.path.isfile("magnitudes"+"_i"+str(i)+"_MODE_"+str(mode)))==False:
                        myfile = open("magnitudes_MODE_"+str(mode), "w")
                        myfile.write('#%5s %5s %9s %11s %11s %11s %11s %11s\n' % ("incl","n_r","lum","u","b","v","r","i"))
                        myfile.close()

                    if (os.path.isfile(modeloc+"walraven_magnitudes"+"_i"+str(i)+"_MODE_"+str(mode)))==False:
                        myfile = open(modeloc+"walraven_magnitudes"+"_i"+str(i)+"_MODE_"+str(mode), "w")
                        myfile.write('#%5s %5s %9s %11s %11s %11s %11s %11s\n' % ("incl","n_r","lum","w","u","l","b","v"))
                        myfile.close()

                    #cmag.calc_mags("outputflux_i"+str(i)+"_MODE_"+str(mode)+"_r"+str(rzone),[i],mode,rzone)
                    cmag.calc_walraven("outputflux_i"+str(i)+"_MODE_"+str(mode)+"_r"+str(rzone),[i],mode,rzone)

                os.chdir(vis_path)
        print mode
#        plt.plot(0.02*xi_t[:,3]/np.max(m4p4_xi_r[-10,:]),"-",alpha=0.8,label=old_modes[modes.index(mode)])
#        plt.grid()
#        plt.xlabel("Radial Zone")
#        plt.ylabel(r"$\delta$T/T")
#        plt.ylabel(r"$\xi_{r}$")
#        plt.title(r"M=1.875M$\odot$ - V=0 - i=40deg")
#        plt.legend(loc="best")

sns.set(style="white",rc={"figure.figsize": (8, 8),'axes.labelsize': 16,
                              'ytick.labelsize': 12,'xtick.labelsize': 12,
                              'legend.fontsize': 12,'axes.titlesize':18,'font.size':14})
start = timer()
#varr = [0,4,6,7,8]
#mselect = "1p875"
#kind = "f"
#r_ord = 0
#mds = list_gen([0,2,3,4,5],r_ord,r_ord,kind)
#mmaxarr = []
#global tempmax
#tempmax = np.array([0,0,0])
#for m in mds:
#    mmax = find_dr_ofmax(emode(mselect,0,m),10,range(10))
#    t = mmax[np.argmax(mmax[:,-1])]
#    if abs(t[2])>tempmax[2]:
#        tempmax = t
#    mmaxarr.append(mmax)
global drtempmax
#drtempmax = 2.28534e+00 #n=5
#drtempmax = 3.3914441477217716 # n=2
drtempmax = 3.111 # n=0
#drtempmax = 4.4445 #8.889 # n=5 M2p5

for mnopq in [3]:
    mselect = "2"
    kind = "p"
    r_ord = mnopq
    varr = [0]
    mds = list_gen([6],r_ord,r_ord,kind)

    #with sns.color_palette("rainbow",n_colors=9):
    for i in varr:
         run_pipeline(mds,sub_dir="",mm=mselect,vv=i)
#        run_pipeline(mds,sub_dir="minL/",minL=True,mm=mselect,vv=i)
##################
#    for j in varr:
#        run_pipeline(mds,sub_dir="phi_max/",mm=mselect,vv=j,tph=123.)
#        run_pipeline(mds,sub_dir="phi_min/",minL=True,mm=mselect,vv=j,tph=123.)
#################
#    for j in varr:
#        run_pipeline(mds,sub_dir="cut_max/",mm=mselect,vv=j,tph=0.,tcut=True)
#        run_pipeline(mds,sub_dir="cut_min/",minL=True,mm=mselect,vv=j,tph=0.,tcut=True)

#if kind !="g":
#    plt.title(r"M="+mselect+"M$_{\odot}$ - $\ell$ = "+mds[0].split()[0] + ", n = "+mds[0].split()[1].strip("p").strip("h").strip("f").strip("F").strip("g"))
#else:
#    plt.title(r"M="+mselect+"M$_{\odot}$ - $\ell$ = "+mds[0].split()[0] + ", n = -"+mds[0].split()[1].strip("p").strip("h").strip("f").strip("F").strip("g"))

#ody = plt.ylim()
#odx = plt.xlim()
#plt.vlines(10,ody[0],ody[1],linestyles="--",alpha=0.5)
#plt.ylim(ody[0],ody[1])
#plt.hlines(0,odx[0],odx[1],linestyles="--",alpha=0.5)
#plt.xlim(odx[0],odx[1])
#plt.title(r"M="+mselect+"M$_{\odot}$ - V = 0 km/s, n = "+mds[0].split()[1].strip("p").strip("h").strip("f").strip("F").strip("g"))
plt.grid()
end = timer()
#print modes_not_found
print "time:",end-start

max_stars_repo_name: dfcastap/visuite

max_issues_repo_name: dfcastap/visuite

max_forks_repo_name: dfcastap/visuite


text:
# -*- coding: utf-8 -*-
import cv2
import numpy

src = cv2.imread('test.png')
mask = cv2.imread('manmade_mask.png')
src = cv2.resize(src,None,fx=0.6, fy=0.6, interpolation = cv2.INTER_CUBIC)
mask = cv2.resize(mask,None,fx=0.6, fy=0.6, interpolation = cv2.INTER_CUBIC)
save = numpy.zeros(src.shape, numpy.uint8) #创建一张空图像用于保存

for row in range(src.shape[0]):
    for col in range(src.shape[1]):
        for channel in range(src.shape[2]):
            if mask[row, col, channel] == 0:
                val = 0
            else:
                reverse_val = 255 - src[row, col, channel]
                val = 255 - reverse_val * 256 / mask[row, col, channel]
                if val < 0: val = 0
            save[row, col, channel] = val

cv2.imshow('src', src)
cv2.imshow('mask', mask)
cv2.imshow('save', save)
cv2.waitKey(0)
cv2.destroyAllWindows()

ext: py

max_stars_repo_name: QGrain/My-Awesome-Tools

max_issues_repo_name: QGrain/My-Awesome-Tools

max_forks_repo_name: QGrain/My-Awesome-Tools

reason: import numpy


text:
# !/usr/bin/hfo_env python3
# encoding utf-8
import argparse
import json
import os
import random

import numpy as np

import settings
from agents.utils import ServerDownError, get_vertices_around_ball
from agents.plastic_dqn_v1.hfo_env.game_interface import HFOAttackingPlayer
from agents.plastic_dqn_v1.agent.dqn import DQN
from agents.plastic_dqn_v1.hfo_env.actions.complex import Actions
from agents.plastic_dqn_v1.hfo_env.features.plastic_features import PlasticFeatures
from agents.plastic_dqn_v1 import config

"""
This module is used to test the agent previous trained
"""

class GameMetrics:
    def __init__(self):
        self.num_ep = 0
        self.num_wins = 0
        self.num_wrong_actions = 0
        self.num_correct_actions = 0
        self.num_games_touched_ball = 0
        self.num_games_passed_ball = 0

    def restart(self):
        self.num_ep = 0
        self.num_wins = 0
        self.num_wrong_actions = 0
        self.num_correct_actions = 0
        self.num_games_touched_ball = 0
        self.num_games_passed_ball = 0

    def inc_num_ep(self):
        self.num_ep += 1

    def inc_num_wins(self):
        self.num_wins += 1

    def inc_num_games_touched_ball(self):
        self.num_games_touched_ball += 1

    def inc_num_games_passed_ball(self):
        self.num_games_passed_ball += 1

    def inc_num_wrong_actions(self):
        self.num_wrong_actions += 1

    def inc_num_correct_actions(self):
        self.num_correct_actions += 1

    def export(self, num_episodes: int, verbose: bool = True):
        correct_actions_rate = \
            round(self.num_correct_actions /
                  (self.num_correct_actions + self.num_wrong_actions),
                  2)
        avr_touched_ball_rate = \
            round(self.num_games_touched_ball / num_episodes, 2)
        avr_passed_ball_rate = round(
            self.num_games_passed_ball / max(self.num_games_touched_ball, 1),
            2)
        win_rate = self.num_wins / num_episodes
        if verbose:
            print(f"[Game Metrics: Summary] WIN rate = {win_rate}; "
                  f"Correct Actions Rate = {correct_actions_rate}; "
                  f"Touched Ball rate = {avr_touched_ball_rate}; "
                  f"Pass Ball rate={avr_passed_ball_rate};")
        return win_rate, correct_actions_rate, avr_touched_ball_rate, \
               avr_passed_ball_rate

    def export_to_dict(self, num_episodes: int, verbose: bool = True):
        win_rt, corr_act_rt, avr_ball_rt, avr_pass_rt = \
            self.export(num_episodes=num_episodes, verbose=verbose)
        return dict(win_rate=win_rt, correct_actions_rate=corr_act_rt,
                    avr_touched_ball_rate=avr_ball_rt,
                    avr_passed_ball_rate=avr_pass_rt)

class Player:
    def __init__(self, num_opponents: int, num_teammates: int, model_file: str,
                 epsilon: int = 1,  port: int = 6000):
        # Game Interface:
        self.game_interface = HFOAttackingPlayer(num_opponents=num_opponents,
                                                 num_teammates=num_teammates,
                                                 port=port)
        self.game_interface.connect_to_server()
        # Features Interface:
        self.features = PlasticFeatures(num_op=num_opponents,
                                        num_team=num_teammates)
        # Actions Interface:
        self.actions = Actions(num_team=num_teammates, features=self.features,
                               game_interface=self.game_interface)
        # Agent instance:
        self.epsilon = epsilon
        self.dqn = DQN.load(load_file=model_file)
        # Metrics:
        self.metrics = GameMetrics()

    def set_starting_game_conditions(self, start_pos: tuple = None,
                                     starts_fixed_position: bool = True,
                                     verbose: bool = False):
        """
        Set starting game conditions. Move for initial position, for example
        """
        if starts_fixed_position:
            if not start_pos:
                ball_pos: list = list(self.features.get_ball_coord())
                starting_corners = get_vertices_around_ball(ball_pos)
                start_pos = random.choice(starting_corners)
            self.actions.move_to_pos(start_pos)
            if verbose:
                print(f"[PLAYER: GAME SET UP] Initial pos= {start_pos}")
        else:
            # Start in current position
            if verbose:
                print(f"[START GAME] Initial pos= RANDOM")
        # Informs the other players that it is ready to start:
        self.game_interface.hfo.say(settings.PLAYER_READY_MSG)

    def exploit_actions(self, state: np.ndarray, verbose: bool = False) -> int:
        q_predict = self.dqn.predict(state)
        max_list = np.where(q_predict == q_predict.max())
        if len(max_list[0]) > 1:
            action = np.random.choice(max_list[0])
        else:
            action = np.argmax(q_predict)
        if verbose:
            print("Q values {} -> {}".format(q_predict, int(action)))
        return int(action)

    def play(self, num_episodes: int, starts_fixed_position: bool = True):
        """
        @param num_episodes: number of episodes to train in this iteration
        @param starts_fixed_position: bool
        @raise ServerDownError
        @return: Game Metrics
        """
        self.metrics.restart()

        for ep in range(num_episodes):
            # Check if server still running:
            try:
                self.game_interface.check_server_is_up()
            except ServerDownError:
                print("!!SERVER DOWN!! Test {}/{}".format(ep, num_episodes))
                metrics_dict = self.metrics.export_to_dict(num_episodes)
                return metrics_dict

            # Update features:
            self.features.update_features(
                observation=self.game_interface.get_observation())

            # Go to origin position:
            self.set_starting_game_conditions(
                starts_fixed_position=starts_fixed_position)

            # metrics:
            touched_ball = False
            passed_ball = False
            while self.game_interface.in_game():
                if self.features.has_ball(): touched_ball = True

                # Update environment features:
                features_array = self.features.get_features()

                # Act:
                act = self.exploit_actions(features_array, verbose=False)
                status, correct_action, passed_ball_succ = \
                    self.actions.execute_action(act)

                # Metrics:
                if passed_ball_succ is True: passed_ball = True
                if correct_action:
                    self.metrics.inc_num_correct_actions()
                else:
                    self.metrics.inc_num_wrong_actions()

            # Update auxiliar variables:
            self.metrics.inc_num_ep()
            if self.game_interface.scored_goal(): self.metrics.inc_num_wins()
            if touched_ball: self.metrics.inc_num_games_touched_ball()
            if passed_ball: self.metrics.inc_num_games_passed_ball()

            # Game Reset
            self.game_interface.reset()

        metrics_dict = self.metrics.export_to_dict(num_episodes)
        return metrics_dict

if __name__ == '__main__':
    parser = argparse.ArgumentParser()
    parser.add_argument('--num_teammates', type=int, default=0)
    parser.add_argument('--num_opponents', type=int, default=0)
    parser.add_argument('--num_episodes', type=int, default=1)
    parser.add_argument('--team_name', type=str, default=None)
    parser.add_argument('--starts_fixed_position', type=str, default="true")
    parser.add_argument('--dir', type=str)
    parser.add_argument('--port', type=int, default=6000)

    # Parse arguments:
    args = parser.parse_args()
    num_team = args.num_teammates
    num_op = args.num_opponents
    num_episodes = args.num_episodes
    team_name = args.team_name
    starts_fixed_position = True if args.starts_fixed_position == "true" \
        else False
    directory = args.dir
    port = args.port

    print(f"[PLASTIC Player: {team_name}] ep={num_episodes}; "
          f"num_t={num_team}; num_op={num_op}; "
          f"starts_fixed_pos={starts_fixed_position};")

    # Model file:
    base_model_file = config.MODEL_FILE_FORMAT.format(team_name=team_name,
                                                      step=step)
    prefix = "new_" if new_model else "re_"
    base_model_file = prefix + base_model_file
    model_file = os.path.join(directory, base_model_file + f".{model_idx}")
    player = Player(num_teammates=num_team, num_opponents=num_op,
                    model_file=model_file)

    # Test Player
    game_metrics: dict = player.play(
        num_episodes=num_episodes,
        starts_fixed_position=starts_fixed_position)

    metrics_file = os.path.join(directory,
                                f"{team_name}_test_metrics_{step}.json")
    # Load metrics data:
    if os.path.isfile(metrics_file):
        with open(metrics_file, "rb") as fp:
            test_metrics = json.load(fp)
    else:
        test_metrics = dict()

    key = "new_models" if new_model else "re_trained_models"
    if test_metrics.get(key):
        test_metrics[key][model_idx] = game_metrics
    else:
        test_metrics[key] = {model_idx: game_metrics}
    with open(metrics_file, 'w+') as fp:
        json.dump(test_metrics, fp)

    print("\n!!!!!!!!! Test End !!!!!!!!!!!!\n\n")


text:
import loopy as lp
import numpy as np

def create_loop_kernel(component_name, domains, instructions, edges, signature):

    domains, assumption_string = create_domain_string(domains, edges)
    globals = create_globals(signature, edges)

    knl = lp.make_function(
        domains,
        instructions,
        globals + ["..."],
        name=component_name,
        assumptions=assumption_string,
        target=lp.CTarget()
    )
    knl = add_instruction_deps(knl)
    if component_name == 'main':
        print(lp.generate_code_v2(knl).device_code())
    return knl

def add_instruction_deps(knl):
    assignees = {}

    for i in knl.instructions:
        assignee = i.assignee_name
        if assignee not in assignees.keys():
            assignees[assignee] = i.id
        else:
            knl = lp.add_dependency(knl, f"id:{i.id}", f"id:{assignees[assignee]}")
            assignees[assignee] = i.id
    return knl

def combine_domains(knl_domains):

    if isinstance(knl_domains, list):
        return _combine_domains_from_list(knl_domains)
    else:
        return _combine_domains_from_kernel(knl_domains)

def _combine_domains_from_list(domain_strs):
    from loopy.kernel.creation import parse_domains
    import islpy as isl
    domains = parse_domains(domain_strs, {})
    result = None
    for dom in domains:
        if result is None:
            result = dom
        else:
            aligned_dom, aligned_result = isl.align_two(
                dom, result, across_dim_types=True)
            result = aligned_result & aligned_dom
    return result

def _combine_domains_from_kernel(knl):
    import islpy as isl
    result = None
    for dom in knl.domains:
        if result is None:
            result = dom
        else:
            aligned_dom, aligned_result = isl.align_two(
                    dom, result, across_dim_types=True)
            result = aligned_result & aligned_dom

    return result
def create_domain_string(domains, edges):
    assumptions = []
    d_str = ["{" + f"[{k}] : {domains[k]}" + "}" for k in domains.keys()]
    for k in domains.keys():
        a = f"{edges[k]['lower']} >= 0"
        if a not in assumptions and not edges[k]['lower'].isdigit():
            assumptions.append(a)
    assumption_str = " and ".join(assumptions)
    return d_str, assumption_str

def create_globals(signature, edges):

    globals = []
    dimension_vars = []
    for s in signature['args']:
        if s in signature['param']:
            g = lp.ValueArg(s, dtype=edges[s]['type'])
            globals.append(g)
        else:
            dims = edges[s]['dimensions']
            for d in dims:
                if d in dimension_vars:
                    continue
                g = lp.ValueArg(d,dtype=np.int32)
                globals.append(g)
                dimension_vars.append(d)
            g = lp.GlobalArg(s, shape=tuple(dims), dtype=edges[s]['type'])
            globals.append(g)

    return globals

def make_component(cname, args, edges):

    dim_vals = []
    for a in args:
        dims = edges[a]['dimensions']
        for d in dims:
            if d in dim_vals:
                continue
            dim_vals.append(d)
    dim_vals = dim_vals + args
    instruction_str = f"{cname}({', '.join(dim_vals)})"
    return instruction_str

def parse_insn(instruction_text):
    import loopy as kc

    inst,_,_ = kc.parse_instructions(instruction_text, {})

    return inst

def add_temporary(stmt):
    return '<>' + stmt

def process_pow(stmt):
    return stmt.replace("^", "**")

def process_if_stmt(prefix, predicate, true_stmt, false_stmt):

    cond_stmt = f"{prefix} = if({predicate}, {true_stmt}, {false_stmt})"
    return cond_stmt

def process_constant(stmt, constant):
    if constant == 'e':
        val = np.e
    elif constant == 'pi':
        val = np.pi
    else:
        val = None
    return stmt.replace(f"{constant}()", str(val))

def process_cdlang_stmt(stmt):

    new_stmt = stmt.replace('][', ',')
    new_stmt = new_stmt.replace(";", "")
    return new_stmt

def make_group_func(lhs, reduce_func, group_expr):
    import re
    from loopy.symbolic import parse
    group_expr = process_cdlang_stmt(group_expr)

    if 'arg' in reduce_func:
        new_arr_expr = reduce_func.replace("][", ",").replace("[", "((").replace("]", "), ", 1)
        indices = re.search('\(\((.*)\)', new_arr_expr).group(1)
        new_arr_expr = new_arr_expr.replace(f"({indices})", f"({indices}),({indices})")
        rhs = f"{new_arr_expr}{group_expr})"
        stmt = f"{lhs}={rhs}"
    else:
        new_arr_expr = reduce_func.replace("][", ",").replace("[", "((").replace("]", "), ", 1)
        rhs = f"{new_arr_expr}{group_expr})"
        stmt = f"{lhs}={rhs}"

    return stmt

def make_lex_schedule(knl):
    import islpy as isl
    from loopy.schedule import RunInstruction
    ctx = knl.isl_context

    for i in knl.instructions:
        insn_id = i.id
        # print(i.depends_on)
    #
    for s in knl.schedule:
        # if isinstance(s, RunInstruction):
        print(s.hash_fields)
        if 'iname' in s.hash_fields:
            print(s.iname)
        # print("\n".join(dir(s)))

def make_domain(prog):

    knl_code = lp.generate_code_v2(prog)

    for id, v in knl_code.implemented_domains.items():
        # print(v[0].compute_schedule())
        print(v[0].to_str())
        print(dir(v[0]))
        # print(f"{id}: {v}")

max_stars_repo_name: he-actlab/cdstack

max_issues_repo_name: he-actlab/cdstack

max_forks_repo_name: he-actlab/cdstack


text:
import torch
import torch.utils.data as data
import torch.nn as nn
import torch.nn.functional as F
import torch.optim as optim
from utils.gtdata import DataAllocate, AudioVisualData
from model.detector import Framework, CAMCluster, ElementLoss_, CAMFeat_, em_cluster, CAMAudio
from utils.logger import AverageMeter, Logger
import sklearn.metrics

import cv2
import librosa
import numpy as np
import os
import pickle
from progress.bar import Bar
import argparse
import time
import shutil
import matplotlib.pyplot as plt
torch.multiprocessing.set_sharing_strategy('file_system')

def parse_args():
    parser = argparse.ArgumentParser(description='AudioVisual Learning')

    parser.add_argument('-d', '--dataset', default='Flickr', type=str)
    parser.add_argument('-j', '--workers', default=0, type=int, metavar='N',
                        help='number of data loading workers (default: 2)')

    parser.add_argument('--epochs', default=15, type=int, metavar='N',
                        help='number of total epochs to run')
    parser.add_argument('--start-epoch', default=0, type=int, metavar='N',
                        help='manual epoch number (useful on restarts)')
    parser.add_argument('--train-batch', default=2, type=int, metavar='N',
                        help='train batchsize')
    parser.add_argument('--val-batch', default=1, type=int, metavar='N',
                        help='test batchsize')
    parser.add_argument('--lr', '--learning-rate', default=0.001, type=float,
                        metavar='LR', help='initial learning rate')
    parser.add_argument('--schedule', type=int, nargs='+', default=[50, 75],
                        help='Decrease learning rate at these epochs.')
    parser.add_argument('--warm-up', dest='wp', default=100, type=int,
                        help='warm up learning rate epoch')
    parser.add_argument('--gamma', type=float, default=0.1, help='LR is multiplied by gamma on schedule.')
    parser.add_argument('--momentum', default=0.9, type=float, metavar='M',
                        help='momentum')
    parser.add_argument('--weight-decay', '--wd', default=5e-4, type=float,
                        metavar='W', help='weight decay (default: 1e-4)')
    parser.add_argument('--mix', default=4, type=int,
                        help='frames per batch for triplet loss')
    parser.add_argument('--rois', default=8, type=int,
                        help='rois per frame for visual stream')
    parser.add_argument('--size', default=3, type=int,
                        help='size for roi-align')
    parser.add_argument('--val-per-epoch','-tp', dest='tp', default=2, type=int,
                        help='number of training epoches between test (default: 30)')
    parser.add_argument('--iter-size', '-is', dest='its', default=4, type=int,
                        help='the forward-backward times within each iteration')

    parser.add_argument('--pretrained', '-pre', dest='pre', default=True, type=bool,
                        help='whether to use pretrained model')
    parser.add_argument('-c', '--checkpoint', default='checkpoint', type=str, metavar='PATH',
                        help='path to save checkpoint (default: checkpoint)')
    parser.add_argument('--resume', default='', type=str, metavar='PATH',
                        help='path to latest checkpoint (default: none)')
    parser.add_argument('-e', '--evaluate', dest='evaluate', action='store_true',
                        help='evaluate model on validation set')

    parser.add_argument('--gpu-id', default='0', type=str,
                        help='id(s) for CUDA_VISIBLE_DEVICES')

    args = parser.parse_args()

    return args

args = parse_args()

os.environ['CUDA_VISIBLE_DEVICES'] = args.gpu_id
use_gpu = torch.cuda.is_available() and int(args.gpu_id) >= 0

with open('utils/gtdata.pkl', 'rb') as f:
    info = pickle.load(f, encoding='bytes')

class Evaluator():

    def __init__(self):
        super(Evaluator, self).__init__()
        self.ciou = []

    def cal_CIOU(self, infer, gtmap, thres=0.01):
        gtmap = gtmap.reshape([256, 256])
        infer = infer.reshape([256, 256])
        infer_map = np.zeros((256, 256))
        infer_map[infer>=thres] = 1
        ciou = np.sum(infer_map*gtmap) / (np.sum(gtmap)+np.sum(infer_map*(gtmap==0)))
        self.ciou.append(ciou)
        return ciou

    def cal_AUC(self):
        results = []
        for i in range(21):
            result = np.sum(np.array(self.ciou)>=0.05*i)
            result = result / len(self.ciou)
            results.append(result)
        x = [0.05*i for i in range(21)]
        auc = sklearn.metrics.auc(x, results)
        print(results)
        return auc

    def final(self):
        ciou = np.mean(np.array(self.ciou)>=0.5)
        return ciou

    def clear(self):
        self.ciou = []

def main():

    f = open('utils/gtest.txt', 'r')
    pool = []
    while True:
        num = f.readline()
        if len(num) == 0:
            break
        pool.append(num.split('.')[0])

    start_epoch = args.start_epoch

    if not os.path.isdir(args.checkpoint):
        os.makedirs(args.checkpoint)

    print('Preparing Dataset %s' % args.dataset)

    valset = AudioVisualData(info, mix=args.mix, pool=pool, training=False)
    valloader = data.DataLoader(valset, batch_size=args.val_batch, shuffle=False, num_workers=args.workers,
                                collate_fn=DataAllocate)

    framework = Framework(args.pre, args.mix, 1, args.rois, args.size)
    cluster = CAMCluster(framework)
    evaluator = Evaluator()
    camfeat = CAMAudio(framework)
    element = ElementLoss_(args.mix, 1)
    cluster_a = np.load('utils/cluster_a.npy')
    cluster_v = np.load('utils/cluster_v.npy')
    cluster_a = torch.FloatTensor(cluster_a.T)
    cluster_v = torch.FloatTensor(cluster_v)

    if use_gpu:
        framework = framework.cuda()
        cluster = cluster.cuda()
        camfeat = camfeat.cuda()
        element = element.cuda()
        cluster_a = cluster_a.cuda()
        cluster_v = cluster_v.cuda()

    args.checkpoint = os.path.join(args.checkpoint, args.dataset)
    if not os.path.exists(args.checkpoint):
        os.mkdir(args.checkpoint)

    if args.resume:
        checkpoint = torch.load(args.resume)
        framework.load_state_dict(checkpoint['state_dict'], strict=True)
        start_epoch = checkpoint['epoch']

    if args.evaluate:
        ciou, auc = test_cam(valloader, framework, camfeat, element, evaluator, cluster_a, cluster_v, start_epoch, use_gpu)
        print(ciou, auc)
        return

def test_avg(dataloader, model, cam, element, evaluator, cluster_a, cluster_v, epoch, use_gpu):
    model.eval()
    model.audio.gru.train()

    data_time = AverageMeter()
    end = time.time()

    bar = Bar('Processing', max=len(dataloader))

    for batch_idx, (audio, visual, roi, gtmap, box) in enumerate(dataloader):
        audio = audio.view(args.val_batch * args.mix, *audio.shape[-2:])
        visual = visual.view(args.val_batch * args.mix, *visual.shape[-3:])
        roi = roi.view(args.val_batch, args.mix, args.rois, 4)
        gt = gtmap.numpy()

        data_time.update(time.time() - end)
        if use_gpu:
            audio = audio.cuda()
            visual = visual.cuda()
            roi = roi.cuda()

        data_time.update(time.time() - end)
        pred_a, feat_a, pred_v, feat_v = model(audio, visual, roi, cluster_a, cluster_v, False)
        feat_v = model.discrim.spa_conv(feat_v)
        feat_v = feat_v.permute([0, 2, 3, 1]).contiguous()

        feat_a = model.discrim.temp_conv(feat_a)
        feat_a = model.discrim.temp_pool(feat_a)

        feat_a = feat_a.view(args.val_batch, 1, 512)
        feat = torch.cat([feat_a.repeat(1, 256, 1),
                          feat_v.view(args.val_batch, 256, 512)], -1)
        feat = feat.view(-1, 1024)
        cams = model.discrim.auto_align(feat)
        cams = cams.view(args.val_batch, 256)
        cams = torch.softmax(-cams*10, -1)
        cams = cams / torch.max(cams, 1)[0].unsqueeze(-1)
        cams = torch.nn.functional.interpolate(cams.view(args.val_batch, 1, 16, 16), (256, 256),
                                               mode='bilinear')
        cams = cams.detach().cpu().numpy()
        for idx in range(args.val_batch):
            cam_visualize(batch_idx * args.val_batch + idx, visual[idx: idx + 1], cams[idx],
                          None, box[idx])
            ciou = evaluator.cal_CIOU(cams[idx], gt[idx], 0.1)
            end = time.time()
        bar.suffix = '({batch}/{size}) Data: {data:.3f}s | CIOU: {ciou:.3f}'.format(
            batch=batch_idx + 1,
            size=len(dataloader),
            data=data_time.val,
            ciou = ciou
        )
        bar.next()

    bar.finish()
    return evaluator.final(), evaluator.cal_AUC()

def test_cam(dataloader, model, cam, element, evaluator, cluster_a, cluster_v, epoch, use_gpu):
    model.eval()
    model.audio.gru.train()

    data_time = AverageMeter()
    end = time.time()

    bar = Bar('Processing', max=len(dataloader))

    for batch_idx, (audio, visual, roi, gtmap, box) in enumerate(dataloader):
        audio = audio.view(args.val_batch * args.mix, *audio.shape[-2:])
        visual = visual.view(args.val_batch * args.mix, *visual.shape[-3:])
        roi = roi.view(args.val_batch, args.mix, args.rois, 4)
        gt = gtmap.numpy()

        data_time.update(time.time() - end)
        if use_gpu:
            audio = audio.cuda()
            visual = visual.cuda()
            roi = roi.cuda()

        data_time.update(time.time() - end)
        pred_a, feat_a, pred_v, feat_v = model(audio, visual, roi, cluster_a, cluster_v, False)
        feat_v = model.discrim.spa_conv(feat_v)
        feat_v = feat_v.permute([0, 2, 3, 1]).contiguous()

        feat_a, ratio, coeff = cam(pred_a, feat_a)
        feat_a = model.discrim.temp_conv(feat_a)
        feat_a = model.avalign.temp_pool(feat_a)

        # coeff = torch.nn.functional.adaptive_max_pool2d(coeff, (feat_a.shape[-2], feat_a.shape[-1]))
        # coeff = coeff.view(*coeff.shape[:2], -1)
        # feat_a = feat_a.view(*feat_a.shape[:2], -1)
        # feat_a = torch.sum(feat_a*coeff, -1) / torch.sum(coeff+1e-10, -1)

        feat_a = feat_a.view(args.val_batch, 7, 1, 512)
        feat = torch.cat([feat_a.repeat(1, 1, 256, 1),
                          feat_v.view(args.val_batch, 1, 256, 512).repeat([1, 7, 1, 1])], -1)
        feat = feat.view(-1, 1024)
        cams = model.discrim.auto_align(feat)
        cams = cams.view(args.val_batch, 7, 256)
        cams = torch.softmax(-cams*10, -1)
        cams = cams / torch.max(cams, -1)[0].unsqueeze(-1)
        cams = torch.nn.functional.interpolate(cams.view(args.val_batch, 7, 16, 16), (256, 256),
                                               mode='bilinear')
        cams = cams.detach().cpu().numpy()
        for idx in range(args.val_batch):
            cam_visualize(batch_idx*args.val_batch+idx, visual[idx: idx+1], cams[idx], None, box[idx])
            index = (ratio[idx].cpu().numpy()>0)
            pred = pred_a[idx].detach().cpu().numpy()
            index[np.argmax(ratio[idx].cpu().numpy())] = 1
            camp = np.sum(cams[idx][index]*pred.reshape(7, 1, 1)[index], 0) / np.sum(pred[index])
            # camp = np.max(cams[idx][index], 0)
            # cam_visualize(batch_idx * args.val_batch + idx, visual[idx: idx + 1], np.expand_dims(camp, 0),
            #               None, box[idx])
            ciou = evaluator.cal_CIOU(camp, gt[idx], 0.1)
            end = time.time()
        bar.suffix = '({batch}/{size}) Data: {data:.3f}s | CIOU: {ciou:.3f}'.format(
            batch=batch_idx + 1,
            size=len(dataloader),
            data=data_time.val,
            ciou = ciou
        )
        bar.next()

    bar.finish()
    return evaluator.final(), evaluator.cal_AUC()

def gt_visualize(batch_idx, visual, gtmap, boxes=None):
    assert visual.shape[0] == gtmap.shape[0]
    segments = visual.shape[0]
    visual = visual.cpu().numpy()
    visual = visual.transpose([0, 2, 3, 1])
    visual = visual * np.array([0.229, 0.224, 0.225])
    visual = visual + np.array([0.485, 0.456, 0.406])
    visual = np.clip(visual, 0.0, 1.0)
    for seg in range(segments):
        actmap = gtmap[seg]
        actmap = np.expand_dims(actmap, -1)
        img = visual[seg]
        if boxes is not None:
            for box in boxes:
                (xmin, ymin, xmax, ymax) = box
                cv2.rectangle(img, (xmin, ymin), (xmax, ymax), (255, 0, 0), 2)
        plt.imsave('./gtvis/'+str(batch_idx)+'_'+str(seg)+'gt.jpg', 0.6*actmap+0.4*img)

def cam_visualize(batch_idx, visual, cam, index=None, boxes=None, thres=None):
    font = cv2.FONT_HERSHEY_SIMPLEX
    segments = cam.shape[0]
    visual = visual.cpu().numpy()
    visual = visual.transpose([0, 2, 3, 1])
    visual = visual * np.array([0.229, 0.224, 0.225])
    visual = visual + np.array([0.485, 0.456, 0.406])
    visual = np.clip(visual, 0.0, 1.0)
    for seg in range(segments):
        actmap = cam[seg]
        if thres is not None:
            actmap[actmap>thres] = 1
        actmap = actmap * 255
        actmap = actmap.astype(np.uint8)
        actmap = cv2.applyColorMap(actmap, cv2.COLORMAP_JET)[:, :, ::-1]
        actmap = actmap / 255.0
        img = visual[0]
        if index is not None:
            img = cv2.putText(img, str(index), (0, 40), font, 1.2, (255, 255, 255), 2)
        # if boxes is not None:
        #     for box in boxes:
        #         (xmin, ymin, xmax, ymax) = box
        #         cv2.rectangle(img, (xmin, ymin), (xmax, ymax), (255, 0, 0), 2)
        plt.imsave('./gtvis/'+str(batch_idx)+'_'+str(seg)+'.jpg', 0.4*actmap+0.6*img)

if __name__ == '__main__':
    main()


text:
import numpy as np
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, ConstantKernel
import logging

class GaussianProcess:
    def __init__(self, numFidelities, dim):
        self.dim = dim
        # Use an anisotropic kernel
        # (independent length scales for each dimension)
        sqrdExp = ConstantKernel() ** 2. * RBF(length_scale=self.dim*[1.])
        numHyperParams = self.dim + 1
        self.models, self.isFit, self.xValues, self.yValues = [], [], [], []
        for _ in range(numFidelities):
            self.models.append(GaussianProcessRegressor(
                                    kernel=sqrdExp,
                                    n_restarts_optimizer=numHyperParams*10))
            self.isFit.append(False)
            self.xValues.append([])
            self.yValues.append([])

    def isValid(self, fidelity):
        return len(self.xValues[fidelity]) >= 2

    def fitModel(self, fidelity):
        if self.isValid(fidelity) and not self.isFit[fidelity]:
            x = np.atleast_2d(self.xValues[fidelity])
            y = np.array(self.yValues[fidelity]).reshape(-1, 1)
            self.models[fidelity].fit(x, y)
            self.isFit[fidelity] = True

    def addSample(self, x, y, fidelity):
        self.xValues[fidelity].append(x)
        self.yValues[fidelity].append(y)
        self.isFit[fidelity] = False

    def getPrediction(self, x, fidelity):
        self.fitModel(fidelity)
        mean, std = self.models[fidelity].predict(x.reshape(-1, self.dim),
                                                  return_std=True)
        return np.array(mean)[:, 0], std

    def plotModel(self, ax, fn, ci):
        assert self.dim == 1
        if not self.isValid(-1):
            return
        import matplotlib.pyplot as plt
        xs = np.linspace(0., 1., 500)
        means, vs = self.getPrediction(xs, -1)
        cs = 1.96 * np.sqrt(vs) # 95% confidence
        lcb, ucb = np.array([ci(x) for x in xs]).T
        ax.set_title('Gaussian Process')
        # ax.plot(xs, means, label='Gaussian Process')
        ax.fill_between(xs, means - cs, means + cs, alpha=.5, color='gray')
        ax.plot(xs, [fn(x) for x in xs], label='f(x)', color='#D95F02')
        ax.plot(xs, lcb, '--', label='LCB', color='#7570B3')
        ax.plot(xs, ucb, '--', label='UCB', color='#1B9E77')
        ax.scatter(self.xValues[-1], self.yValues[-1], label='Samples',
                                                       color='blue')
        ax.legend()
        ax.set_xlim([0., 1.])


text:
from __future__ import division
import argparse
# from . import main, utils
import pkg_resources
import pandas as pd
import ast
import subprocess
import numpy as np

def run():
    parser = argparse.ArgumentParser(
        description='run diffacto for scavager results',
        epilog='''

    Example usage
    -------------
    $ scav2diffacto -S1 sample1_1_proteins.tsv sample1_n_proteins.tsv -S2 sample2_1_proteins.tsv sample2_n_proteins.tsv
    -------------
    ''',
        formatter_class=argparse.ArgumentDefaultsHelpFormatter)

    parser.add_argument('-dif', help='path to Diffacto', required=True)
    parser.add_argument('-S1', nargs='+', help='input files for S1 sample', required=True)
    parser.add_argument('-S2', nargs='+', help='input files for S2 sample', required=True)
    parser.add_argument('-S3', nargs='+', help='input files for S3 sample')
    parser.add_argument('-S4', nargs='+', help='input files for S4 sample')
    parser.add_argument('-peptides', help='name of output peptides file', default='peptides.txt')
    parser.add_argument('-samples', help='name of output samples file', default='sample.txt')
    parser.add_argument('-out', help='name of diffacto output file', default='diffacto_out.txt')
    parser.add_argument('-norm', help='normalization method. Can be average, median, GMM or None', default='None')
    parser.add_argument('-impute_threshold', help='impute_threshold for missing values fraction', default='0.25')
    parser.add_argument('-min_samples', help='minimum number of samples for peptide usage', default='3')
    parser.add_argument('-version', action='version', version='%s' % (pkg_resources.require("scavager")[0], ))
    args = vars(parser.parse_args())

    df_final = False

    allowed_prots = set()
    allowed_peptides = set()

    for sample_num in ['S1', 'S2', 'S3', 'S4']:
        if args[sample_num]:
            for z in args[sample_num]:
                df0 = pd.read_table(z)
                allowed_prots.update(df0['dbname'])

    for sample_num in ['S1', 'S2', 'S3', 'S4']:
        if args[sample_num]:
            for z in args[sample_num]:
                df0 = pd.read_table(z.replace('_proteins.tsv', '_peptides.tsv'))
                allowed_peptides.update(df0['peptide'])

    for sample_num in ['S1', 'S2', 'S3', 'S4']:
        if args[sample_num]:
            for z in args[sample_num]:
                label = z.replace('_proteins.tsv', '')
                df1 = pd.read_table(z.replace('_proteins.tsv', '_PSMs_full.tsv'))
                # df1 = pd.read_table(z.replace('_proteins.tsv', '_PSMs.tsv'))
                df1 = df1[df1['peptide'].apply(lambda z: z in allowed_peptides)]
                # print(df1.shape)
                # print(z.replace('_proteins.tsv', '_PSMs_full.tsv'))
                # print(df1.columns)
                df1['peptide'] = df1.apply(lambda z: z['peptide'] + str(z['assumed_charge']), axis=1)
                df1 = df1.sort_values('MS1Intensity', ascending=False).drop_duplicates(['peptide'])
                df1['peptide'] = df1['peptide'].apply(lambda z: z[:-1])
                df1['MS1Intensity'] = df1.groupby('peptide')['MS1Intensity'].transform(sum)
                df1 = df1.sort_values('q', ascending=True).drop_duplicates(['peptide'])
                df1[label] = df1['MS1Intensity']
                df1[label] = df1[label].replace([0, 0.0], np.nan)
                df1['protein'] = df1['protein'].apply(lambda z: ';'.join([u for u in ast.literal_eval(z) if u in allowed_prots]))
                df1 = df1[df1['protein'].apply(lambda z: z != '')]
                df1 = df1[['peptide', 'protein', label]]
                if df_final is False:
                    df_final = df1
                else:
                    df_final = df_final.reset_index().merge(df1.reset_index(), on='peptide', how='outer')#.set_index('peptide')
                    # df_final = df_final.merge(df1, on='peptide', how='outer')
                    df_final.protein_x.fillna(value=df_final.protein_y, inplace=True)
                    df_final['protein'] = df_final['protein_x']
                    df_final = df_final.drop(columns=['protein_x', 'protein_y', 'index_x', 'index_y'])

    print(df_final.columns)
    df_final = df_final.set_index('peptide')
    df_final['proteins'] = df_final['protein']
    df_final = df_final.drop(columns=['protein'])
    cols = df_final.columns.tolist()
    cols.remove('proteins')
    cols.insert(0, 'proteins')
    df_final = df_final[cols]
    df_final.fillna(value='')
    df_final.to_csv(args['peptides'], sep=',')

    out = open(args['samples'], 'w')
    for sample_num in ['S1', 'S2', 'S3', 'S4']:
        if args[sample_num]:
            for z in args[sample_num]:
                label = z.replace('_proteins.tsv', '')
                out.write(label + '\t' + sample_num + '\n')
    out.close()

    subprocess.call(['python3', args['dif'], '-i', args['peptides'], '-samples', args['samples'], '-out',\
     args['out'], '-normalize', args['norm'], '-impute_threshold', args['impute_threshold'], '-min_samples', args['min_samples']])

if __name__ == '__main__':
    run()

max_stars_repo_name: markmipt/scavager

max_issues_repo_name: markmipt/scavager

max_forks_repo_name: markmipt/scavager


text:
import numpy as np
import torch
import torch.nn as nn
import torch.nn.functional as F

def get_res_net_18(ensemble, **kwargs):

    if ensemble is None:
        return ResNet18(**kwargs)

    elif ensemble == "early_exit":
        return ResNet18EarlyExit(**kwargs)

    elif ensemble == "mc_dropout":
        return ResNet18MCDrop(**kwargs)

    elif ensemble == "deep":
        return ResNet18(**kwargs)

    elif ensemble == "depth":
        return ResNet18Depth(**kwargs)

    else:
        NotImplementedError("ensemble not implemented: '{}'".format(ensemble))

def init_weights(model):

    for module in model.modules():

        if isinstance(module, nn.Conv1d):
            nn.init.kaiming_normal_(module.weight, mode="fan_out", nonlinearity="relu")

        elif isinstance(module, nn.BatchNorm1d):
            nn.init.constant_(module.weight, 1)
            nn.init.constant_(module.bias, 0)

def conv3x3(in_planes, out_planes, stride=1):
    return nn.Conv1d(in_planes, out_planes, kernel_size=3, stride=stride, padding=1, bias=False)

def _conv1x1(in_planes, out_planes, stride=1):
    return nn.Conv1d(in_planes, out_planes, kernel_size=1, stride=stride, bias=False)

class BasicBlock(nn.Module):

    def __init__(self, inplanes, planes, stride=1, downsample=None):
        super().__init__()

        self.conv1 = conv3x3(inplanes, planes, stride)
        self.bn1 = nn.BatchNorm1d(planes)
        self.relu = nn.ReLU(inplace=True)
        self.conv2 = conv3x3(planes, planes)
        self.bn2 = nn.BatchNorm1d(planes)
        self.downsample = downsample
        self.stride = stride

    def forward(self, x):

        identity = x
        out = self.conv1(x)
        out = self.bn1(out)
        out = self.relu(out)
        out = self.conv2(out)
        out = self.bn2(out)

        if self.downsample is not None:
            identity = self.downsample(x)

        out += identity
        out = self.relu(out)

        return out

class ResNet18(nn.Module):

    name = "res_net_18"

    def __init__(self, out_channels, seed=None):
        super().__init__()

        self.out_channels = out_channels
        self.seed = seed

        self.hidden_sizes = [64, 128, 256, 512]
        self.layers = [2, 2, 2, 2]
        self.strides = [1, 2, 2, 2]
        self.inplanes = self.hidden_sizes[0]

        in_block = [nn.Conv1d(1, self.inplanes, kernel_size=7, stride=2, padding=3, bias=False)]
        in_block += [nn.BatchNorm1d(self.inplanes)]
        in_block += [nn.ReLU(inplace=True)]
        in_block += [nn.MaxPool1d(kernel_size=3, stride=2, padding=1)]
        self.in_block = nn.Sequential(*in_block)

        blocks = []
        for h, l, s in zip(self.hidden_sizes, self.layers, self.strides):
            blocks += [self._make_layer(h, l, s)]
        self.blocks = nn.Sequential(*blocks)

        out_block = [nn.AdaptiveAvgPool1d(1)]
        out_block += [nn.Flatten(1)]
        out_block += [nn.Linear(self.hidden_sizes[-1], self.out_channels)]
        self.out_block = nn.Sequential(*out_block)

        if self.seed is not None:
            torch.manual_seed(seed)

        self.apply(init_weights)

    def _make_layer(self, planes, blocks, stride=1):

        downsample = None

        if stride != 1 or self.inplanes != planes:
            downsample = nn.Sequential(_conv1x1(self.inplanes, planes, stride), nn.BatchNorm1d(planes))

        layers = [BasicBlock(self.inplanes, planes, stride, downsample)]
        self.inplanes = planes

        for _ in range(1, blocks):
            layers += [BasicBlock(self.inplanes, planes)]

        return nn.Sequential(*layers)

    def forward(self, x):

        x = self.in_block(x)
        x = self.blocks(x)
        x = self.out_block(x)

        return x

class ExitBlock(nn.Module):

    def __init__(self, in_channels, hidden_sizes, out_channels):
        super().__init__()

        layers = [nn.AdaptiveAvgPool1d(1)]
        layers += [nn.Flatten(1)]
        layers += [nn.Linear(in_channels, hidden_sizes)]
        layers += [nn.ReLU()]
        layers += [nn.Linear(hidden_sizes, out_channels)]
        self.layers = nn.Sequential(*layers)

    def forward(self, x):

        return self.layers(x)

class ResNet18EarlyExit(ResNet18):

    name = "res_net_18_early_exit"

    def __init__(self, *args, exit_after=-1, complexity_factor=1.2, **kwargs):
        self.exit_after = exit_after
        self.complexity_factor = complexity_factor

        super().__init__(*args, **kwargs)

        to_exit = [2, 8, 15, 24, 31, 40, 47, 56]
        hidden_sizes = len(self.hidden_sizes)

        num_hidden = len(self.hidden_sizes)
        exit_hidden_sizes = [int(((self.complexity_factor ** 0.5) ** (num_hidden - idx)) * self.hidden_sizes[-1]) for idx in range(num_hidden)]
        exit_hidden_sizes = [h for pair in zip(exit_hidden_sizes, exit_hidden_sizes) for h in pair]

        if self.exit_after == -1:
            self.exit_after = range(len(to_exit))

        num_exits = len(to_exit)

        if (len(self.exit_after) > num_exits) or not set(self.exit_after).issubset(list(range(num_exits))):
            raise ValueError("valid exit points: {}".format(", ".join(str(n) for n in range(num_exits))))

        self.exit_hidden_sizes = np.array(exit_hidden_sizes)[self.exit_after]

        blocks = []
        for idx, module in enumerate(self.blocks.modules()):
            if idx in to_exit:
                blocks += [module]
        self.blocks = nn.ModuleList(blocks)

        idx = 0
        exit_blocks = []
        for block_idx, block in enumerate(self.blocks):
            if block_idx in self.exit_after:
                in_channels = block.conv1.out_channels
                exit_blocks += [ExitBlock(in_channels, self.exit_hidden_sizes[idx], self.out_channels)]
                idx += 1
        self.exit_blocks = nn.ModuleList(exit_blocks)

        self.apply(init_weights)

    def forward(self, x):

        out = self.in_block(x)

        out_blocks = []
        for block in self.blocks:
            out = block(out)
            out_blocks += [out]

        out_exits = []
        for exit_after, exit_block in zip(self.exit_after, self.exit_blocks):
            out = exit_block(out_blocks[exit_after])
            out_exits += [out]

        out = self.out_block(out_blocks[-1])
        out = torch.stack(out_exits + [out], dim=1)

        return out

class MCDropout(nn.Dropout):

    def forward(self, x):
        return F.dropout(x, self.p, True, self.inplace)

class ResNet18MCDrop(ResNet18EarlyExit):

    name = "res_net_18_mc_drop"

    def __init__(self, *args, drop_after=-1, drop_prob=0.2, **kwargs):
        self.drop_after = drop_after
        self.drop_prob = drop_prob

        super().__init__(*args, exit_after=drop_after, **kwargs)

        self.drop_after = self.exit_after

        self.__delattr__("exit_after")
        self.__delattr__("exit_blocks")

        for block_idx in self.drop_after:
            self.blocks[block_idx].add_module("dropout", MCDropout(self.drop_prob))

    def forward(self, x):

        x = self.in_block(x)
        x = self.blocks(x)
        x = self.out_block(x)

        return x

class ResNet18Depth(ResNet18):

    name = "res_net_18_depth"

    def __init__(self, *args, max_depth=1, **kwargs):
        self.max_depth = max_depth

        super().__init__(*args, **kwargs)

        num_blocks = len(self.hidden_sizes)

        if self.max_depth == -1:
            self.max_depth = len(self.hidden_sizes)

        elif (max_depth > num_blocks) or (max_depth < 1):
            raise ValueError("valid depths: {}".format(", ".join(str(n) for n in range(1, num_blocks + 1))))

        self.blocks = self.blocks[:self.max_depth]

        out_block = [nn.AdaptiveAvgPool1d(1)]
        out_block += [nn.Flatten(1)]
        out_block += [nn.Linear(self.hidden_sizes[self.max_depth - 1], self.out_channels)]
        self.out_block = nn.Sequential(*out_block)

max_stars_repo_name: ajrcampbell/early-exit-ensembles

max_issues_repo_name: ajrcampbell/early-exit-ensembles

max_forks_repo_name: ajrcampbell/early-exit-ensembles


text:
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns

# Function to remove the columns having percentage of null values > 45%.
def remove_nan_columns(df, percent=0.45):
    new_df_shape = df.shape
    colnames = (df.isna().sum() / len(df))
    colnames = list(colnames[colnames.values >= percent].index)
    df.drop(labels=colnames, axis=1, inplace=True)

    print("Number of columns dropped\t: ", len(colnames))
    print("\nOld dataset rows, columns", new_df_shape, "\nNew dataset rows, columns", df.shape)

    return df

# Function to fill all NAN's with the mode of the choosen column from dataframe.
def replace_nan_with_mode(col):
    # Read the column to see what values it have with dropna = False to show NAN's too.
    print(f'Values before NAN replacement:\n {col.value_counts(dropna=False)}')
    # Get the mode of the column (most frequent value).
    print(f'Mode value: {col.mode()[0]}')
    # Fill NAN with mode value.
    col = col.fillna(col.mode()[0])
    '''
    We can see now that value of mode column is increased with the value of NAN's and there are no NAN's.
    All other values remains the same.
    '''
    print(f'Values after NAN replacement:\n {col.value_counts()}')

    return col

# Function to fill all NAN's with the median of the choosen column from dataframe.
def replace_nan_with_median(col):
    # Read the column to see what values it have with dropna = False to show NAN's too.
    print(f'Values before NAN replacement:\n {col.value_counts(dropna=False)}')
    # Get the median of the column (middle value).
    print(f'Median value: {col.median()}')
    # Fill NAN with mode value.
    col = col.fillna(col.median())
    '''
    We can see now that value of median column is increased with the value of NAN's and there are no NAN's.
    All other values remains the same.
    '''
    print(f'Values after NAN replacement:\n {col.value_counts(dropna=False)}')

    return col

# Function for correlations in dataframe.
def get_correlations_greater_than_or_equal_to_05(df):
    # Create correlation matrix.
    corr = round(df.corr(method="pearson"), 2)

    # Select upper triangle of correlation matrix.
    corr = corr.where(np.triu(np.ones(corr.shape), k=1).astype(np.bool))

    # Unstack and reset index.
    corr_new = corr.unstack().reset_index()

    # Set column names.
    corr_new.columns = ["variable1", "variable2", "correlation"]

    # Drop NaN values.
    corr_new.dropna(subset=["correlation"], inplace=True)

    # Convert all negative correlation to their ABS values.
    corr_new["corr_abs"] = corr_new["correlation"].abs()

    # Filter the data - show only highly correlated data greater than or equal to 0.5.
    corr_new = corr_new[corr_new["correlation"] >= 0.5]

    # Show correlations
    return corr_new.sort_values('corr_abs', ascending=False).reset_index(drop=True)

# Function for correlations in dataframe.
def get_correlations_less_than_05(df):
    # Create correlation matrix.
    corr = round(df.corr(method="pearson"), 2)

    # Select upper triangle of correlation matrix.
    corr = corr.where(np.triu(np.ones(corr.shape), k=1).astype(np.bool))

    # Unstack and reset index.
    corr_new = corr.unstack().reset_index()

    # Set column names.
    corr_new.columns = ["variable1", "variable2", "correlation"]

    # Drop NaN values.
    corr_new.dropna(subset=["correlation"], inplace=True)

    # Convert all negative correlation to their ABS values.
    corr_new["corr_abs"] = corr_new["correlation"].abs()

    # Filter the data - show only highly correlated data less than 0.5.
    corr_new = corr_new[corr_new["correlation"] < 0.5]

    # Show correlations
    return corr_new.sort_values('corr_abs', ascending=False).reset_index(drop=True)

# Function to plot boxplot charts.
def plot_boxplot_chart(input_data, input_label, input_title):
    fig = plt.figure(figsize=(7, 5))
    plt.boxplot(input_data)

    plt.ylabel(input_label)
    plt.yscale('log')
    plt.grid()
    plt.title(input_title, pad=20)

    plt.show()

# Function to plot scatterplot charts.
def plot_scatter_chart(input_data_1, input_data_2, input_xlabel, input_ylabel, input_title, input_color):
    plt.figure(figsize=(10, 5))
    plt.scatter(input_data_1, input_data_2, c=input_color)

    plt.title(input_title, pad=20)
    plt.xlabel(input_xlabel)
    plt.ylabel(input_ylabel)

    plt.show()

# Function to plot pie charts.
def plot_pie_chart(input_data, input_labels, input_colors, input_title, input_explode):
    if input_colors != []:
        plt.pie(input_data,
                autopct="%1.2f%%",
                colors=input_colors,
                explode=input_explode,
                shadow=True)
    else:
        plt.pie(input_data,
                autopct="%1.2f%%",
                explode=input_explode,
                shadow=True)
    #     circ = plt.Circle((0,0),.9,color="white")
    #     plt.gca().add_artist(circ)

    plt.title(input_title)
    plt.legend(input_labels, loc="upper right")

    # plt.show()

# Function to plot horizontal bar charts.
def plot_barh_chart(input_data_1, input_data_2, input_color, input_xlabel, input_title):
    fig = plt.figure(figsize=(10, 7))
    plt.barh(input_data_1, input_data_2, color=input_color)

    plt.xlabel(input_xlabel)
    plt.title(input_title, pad=20)

    plt.show()

# Function to plot bar charts.
def plot_bar_chart(input_data_1, input_data_2, input_color, input_ylabel, input_title):
    fig = plt.figure(figsize=(15, 7))
    plt.bar(input_data_1, input_data_2, color=input_color)

    plt.ylabel(input_ylabel)
    plt.tick_params(rotation=90)
    plt.title(input_title, pad=20)

    plt.show()

# Function to plot pie subplots
def plot_two_pie_subplots(data_1, labels_1, colors_1, title_1, data_2, labels_2, colors_2, title_2, explode):
    fig, ax1 = plt.subplots(1, 2, figsize=(15, 10))

    plt.subplot(121)
    plot_pie_chart(data_1,
                   labels_1,
                   colors_1,
                   title_1,
                   explode
                   )

    plt.subplot(122)
    plot_pie_chart(data_2,
                   labels_2,
                   colors_2,
                   title_2,
                   explode
                   )

    plt.show()

def plot_hist_boxplot_subplots(
        hist1_data, hist2_data,
        box_data,
        title, labels, xlabel, ylabel
):
    fig, ax1 = plt.subplots(1, 2, figsize=(15, 10))

    plt.subplot(121)
    plt.hist(hist1_data, bins='fd', density=True, alpha=0.7)
    plt.hist(hist2_data, bins='fd', density=True, alpha=0.7)

    plt.xlabel(xlabel)
    plt.ylabel(ylabel)
    plt.title(title, pad=20)
    plt.legend(labels, loc='upper right')

    plt.subplot(122)
    plt.boxplot(box_data)

    plt.grid()
    plt.title(title, pad=20)
    plt.labels = labels
    plt.xticks([1, 2], labels)
    plt.xlabel(xlabel)
    plt.ylabel(ylabel)

    plt.show()

def plot_hist_boxplot_subplots_without_density(
        hist1_data, hist2_data,
        box_data,
        title, labels, xlabel, ylabel
):
    fig, ax1 = plt.subplots(1, 2, figsize=(15, 10))

    plt.subplot(121)
    plt.hist(hist1_data, bins='fd', alpha=0.7)
    plt.hist(hist2_data, bins='fd', alpha=0.7)

    plt.xlabel(xlabel)
    plt.ylabel(ylabel)
    plt.title(title, pad=20)
    plt.legend(labels, loc='upper right')
    plt.yscale("log")

    plt.subplot(122)
    plt.boxplot(box_data, showfliers=False)

    plt.grid()
    plt.title(title, pad=20)
    plt.labels = labels
    plt.xticks([1, 2], labels)
    plt.xlabel(xlabel)
    plt.ylabel(ylabel)
    plt.yscale("log")

    plt.show()

# Function to boxplot bivariate
def plot_seaborn_boxplot(col1, col2, df1, df2, title1, title2, outliers=False, input_ylim=None):
    plt.figure(figsize=(25, 10))

    plt.subplot(121)
    plt.title(title1)
    sns.boxplot(x=col1, y=col2, data=df1, showfliers=outliers)
    if input_ylim:
        plt.ylim(input_ylim)

    plt.subplot(122)
    plt.title(title2)
    sns.boxplot(x=col1, y=col2, data=df2, showfliers=outliers)
    if input_ylim:
        plt.ylim(input_ylim)

    plt.show()

# Function to plot seaborn heatmaps
def plot_heatmap_subplots(input_data1, input_data2, list_data, group_data, title1, title2):
    plt.figure(figsize=[20, 5])

    plt.subplot(121)
    plt.title(title1)
    res = pd.pivot_table(data=input_data1,
                         values=list_data,
                         index=group_data,
                         aggfunc=np.mean)
    sns.heatmap(res, annot=True, cmap="Spectral", center=0.117)

    plt.subplot(122)
    plt.title(title2)
    res = pd.pivot_table(data=input_data2,
                         values=list_data,
                         index=group_data,
                         aggfunc=np.mean)
    sns.heatmap(res, annot=True, cmap="Spectral", center=0.117)

    plt.show()

max_stars_repo_name: PetkoAndreev/Python-Data-Science

max_issues_repo_name: PetkoAndreev/Python-Data-Science

max_forks_repo_name: PetkoAndreev/Python-Data-Science


text:
#!/usr/bin/env python3
import numpy as np
from typing import Tuple

def prepare_data(x: np.ndarray, y: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
    """ Prepare X and y data. """
    try:
        x = np.array(x, dtype=np.float32)
        y = np.array(y, dtype=np.float32)
        assert x.ndim == y.ndim == 1
        assert x.shape == y.shape
        return x, y
    except Exception as e:
        raise ValueError("x and y must be array-like objects: {0}".format(e))

def r2_score(y: np.ndarray, yhat: np.ndarray) -> float:
    """ Calculate goodness of fit, R^2. """
    rss = np.sum(np.power(y - yhat, 2))
    tss = np.sum(np.power(y - np.mean(y), 2))
    r2 = 1 - (rss / tss)
    return r2

max_stars_repo_name: Dascienz/tensorfit

max_issues_repo_name: Dascienz/tensorfit

max_forks_repo_name: Dascienz/tensorfit


text:
""" Line search using backtracking."""

import numpy as np

from openmdao.core.system import AnalysisError
from openmdao.solvers.solver_base import LineSearch
from openmdao.util.record_util import update_local_meta, create_local_meta

class BackTracking(LineSearch):
    """A line search subsolver using backtracking.

    Options
    -------
    options['atol'] :  float(1e-10)
        Absolute convergence tolerance for line search.
    options['err_on_maxiter'] : bool(False)
        If True, raise an AnalysisError if not converged at maxiter.
    options['iprint'] :  int(0)
        Set to 0 to disable printing, set to 1 to print the residual to stdout
        each iteration, set to 2 to print subiteration residuals as well.
    options['maxiter'] :  int(10)
        Maximum number of line searches.
    options['rtol'] :  float(0.9)
        Relative convergence tolerancee for line search.
    options['solve_subsystems'] :  bool(True)
        Set to True to solve subsystems. You may need this for solvers nested under Newton.
    """

    def __init__(self):
        super(BackTracking, self).__init__()

        opt = self.options
        opt.add_option('atol', 1e-10, lower=0.0,
                       desc='Absolute convergence tolerancee for line search.')
        opt.add_option('rtol', 0.9, lower=0.0,
                       desc='Relative convergence tolerancee for line search.')
        opt.add_option('maxiter', 0, lower=0,
                       desc='Maximum number of line searches.')
        opt.add_option('solve_subsystems', True,
                       desc='Set to True to solve subsystems. You may need this for solvers nested under Newton.')

        self.print_name = 'BK_TKG'

    def solve(self, params, unknowns, resids, system, solver, alpha, fnorm0,
              metadata=None):
        """ Take the gradient calculated by the parent solver and figure out
        how far to go.

        Args
        ----
        params : `VecWrapper`
            `VecWrapper` containing parameters. (p)

        unknowns : `VecWrapper`
            `VecWrapper` containing outputs and states. (u)

        resids : `VecWrapper`
            `VecWrapper` containing residuals. (r)

        system : `System`
            Parent `System` object.

        metadata : dict, optional
            Dictionary containing execution metadata (e.g. iteration coordinate).

        solver : `Solver`
            Parent solver instance.

        alpha : float
            Initial over-relaxation factor as used in parent solver.

        fnorm0 : float
            Initial norm of the residual for relative tolerance check.

        Returns
        --------
        float
            Norm of the final residual
        """

        atol = self.options['atol']
        rtol = self.options['rtol']
        maxiter = self.options['maxiter']
        result = system.dumat[None]
        local_meta = create_local_meta(metadata, system.pathname)

        # Allow different alphas for each value so we can keep moving when we
        # hit a bound.
        alpha = alpha*np.ones(len(unknowns.vec))

        # If our step will violate any upper or lower bounds, then reduce
        # alpha so that we only step to that boundary.
        alpha = unknowns.distance_along_vector_to_limit(alpha, result)

        # Apply step that doesn't violate bounds
        unknowns.vec += alpha*result.vec

        # Metadata update
        update_local_meta(local_meta, (solver.iter_count, 0))

        # Just evaluate the model with the new points
        if solver.options['solve_subsystems']:
            system.children_solve_nonlinear(local_meta)
        system.apply_nonlinear(params, unknowns, resids, local_meta)

        self.recorders.record_iteration(system, local_meta)

        # Initial execution really belongs to our parent driver's iteration,
        # so use its info.
        fnorm = resids.norm()
        if solver.options['iprint'] > 0:
            self.print_norm(solver.print_name, system.pathname, solver.iter_count,
                            fnorm, fnorm0)

        itercount = 0
        ls_alpha = alpha

        # Further backtacking if needed.
        while itercount < maxiter and \
              fnorm > atol and \
              fnorm/fnorm0 > rtol:

            ls_alpha *= 0.5
            unknowns.vec -= ls_alpha*result.vec
            itercount += 1

            # Metadata update
            update_local_meta(local_meta, (solver.iter_count, itercount))

            # Just evaluate the model with the new points
            if self.options['solve_subsystems']:
                system.children_solve_nonlinear(local_meta)
            system.apply_nonlinear(params, unknowns, resids, local_meta)

            solver.recorders.record_iteration(system, local_meta)

            fnorm = resids.norm()
            if self.options['iprint'] > 0:
                self.print_norm(self.print_name, system.pathname, itercount,
                                fnorm, fnorm0, indent=1, solver='LS')

        if itercount >= maxiter and self.options['err_on_maxiter']:
            raise AnalysisError("Solve in '%s': BackTracking failed to converge after %d "
                                "iterations." % (system.pathname, maxiter))

        return fnorm

max_stars_repo_name: colinxs/OpenMDAO

max_issues_repo_name: colinxs/OpenMDAO

max_forks_repo_name: colinxs/OpenMDAO


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

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1v695u-ZH-SzGLFUfiInyZz_aBHgsKGRi

video.ittensive.com/machine-learning/prudential/train.csv.gz
"""

import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.metrics import cohen_kappa_score, confusion_matrix, make_scorer
from sklearn.neighbors import KNeighborsClassifier
from sklearn.linear_model import LogisticRegression
from sklearn import preprocessing
from sklearn.model_selection import GridSearchCV

data = pd.read_csv("https://video.ittensive.com/machine-learning/prudential/train.csv.gz")

data["Product_Info_2_1"] = data["Product_Info_2"].str.slice(0, 1)
data["Product_Info_2_2"] = pd.to_numeric(data["Product_Info_2"].str.slice(1, 2))
data.drop(labels=["Product_Info_2"], axis=1, inplace=True)
for l in data["Product_Info_2_1"].unique():
    data["Product_Info_2_1" + l] = data["Product_Info_2_1"].isin([l]).astype("int8")
data.drop(labels=["Product_Info_2_1"], axis=1, inplace=True)
data.fillna(value=-1, inplace=True)

columns_groups = ["Insurance_History", "InsurеdInfo", "Medical_Keyword",
                  "Family_Hist", "Medical_History", "Product_Info"]
columns = ["Wt", "Ht", "Ins_Age", "BMI"]
for cg in columns_groups:
    columns.extend(data.columns[data.columns.str.startswith(cg)])

scaler = preprocessing.StandardScaler()
data_transformed = pd.DataFrame(scaler.fit_transform(pd.DataFrame(data,
                                            columns=data.columns)))
columns_transformed = data_transformed.columns
data_transformed["Response"] = data["Response"]

def reduce_mem_usage (df):
    start_mem = df.memory_usage().sum() / 1024**2
    for col in df.columns:
        col_type = df[col].dtypes
        if str(col_type)[:5] == "float":
            c_min = df[col].min()
            c_max = df[col].max()
            if c_min > np.finfo("f2").min and c_max < np.finfo("f2").max:
                df[col] = df[col].astype(np.float16)
            elif c_min > np.finfo("f4").min and c_max < np.finfo("f4").max:
                df[col] = df[col].astype(np.float32)
            else:
                df[col] = df[col].astype(np.float64)
        elif str(col_type)[:3] == "int":
            c_min = df[col].min()
            c_max = df[col].max()
            if c_min > np.iinfo("i1").min and c_max < np.iinfo("i1").max:
                df[col] = df[col].astype(np.int8)
            elif c_min > np.iinfo("i2").min and c_max < np.iinfo("i2").max:
                df[col] = df[col].astype(np.int16)
            elif c_min > np.iinfo("i4").min and c_max < np.iinfo("i4").max:
                df[col] = df[col].astype(np.int32)
            elif c_min > np.iinfo("i8").min and c_max < np.iinfo("i8").max:
                df[col] = df[col].astype(np.int64)
        else:
            df[col] = df[col].astype("category")
    end_mem = df.memory_usage().sum() / 1024**2
    print('Потребление памяти меньше на', round(start_mem - end_mem, 2), 'Мб (минус', round(100 * (start_mem - end_mem) / start_mem, 1), '%)')
    return df

data_transformed = reduce_mem_usage(data_transformed)

data_train, data_test = train_test_split(data_transformed,
                                         test_size=0.2)

"""KNN"""

model_knn = KNeighborsClassifier(n_neighbors=100)

y = data_train['Response']
x = data_train.drop('Response',axis=1)

model_knn.fit(x,y)

test_data = data_test.drop('Response',axis=1)

test_data['target'] = model_knn.predict(test_data)

cohen_kappa_score(test_data['target'],data_test['Response'])

"""Multiple Log Regs
### Логистическая регрессия
В обучающих данных пометим все классы, кроме 6 и 8, как 0 - и проведем обучение по такому набору данных.

Затем в оставшихся данных (в которых класс не равен 6 или 8) заменим все классы, кроме 7 и 1, на 0 - и снова проведем обучение. И т.д. Получим иерархию классификаторов:
8/6/нет -> 7/1/нет -> 2/5/нет -> 4/3
"""

def regression_model(df,columns):
  x = pd.DataFrame(df,columns=columns)
  y = df['Response']
  model = LogisticRegression(max_iter=10000)
  model.fit(x,y)
  return model
def logistic_regression(columns,df_train):
  model = regression_model(df_train,columns)
  logr_model = GridSearchCV(model, {},cv=5,n_jobs=2,scoring=make_scorer(cohen_kappa_score))
  x = pd.DataFrame(df_train,columns=columns)
  y = df_train['Response']
  logr_model.fit(x,y)
  return logr_model.best_score_

def find_opt_params(data_train):
  kappa_score_opt = 0
  columns_opt = []
  for col in columns_transformed:
    kappa_score = logistic_regression([col],data_train)
    if kappa_score > kappa_score_opt:
      columns_opt = [col]
      kapp_score = kappa_score_opt
  for col in columns_transformed:
    if col not in columns_opt:
      columns_opt.append(col)
      kappa_score = logistic_regression(columns_opt,data_train)
      if kappa_score < kappa_score_opt:
        columns_opt.pop()
      else:
        kappa_score_opt = kappa_score
  return kappa_score_opt,columns_opt

data_transformed = data_transformed.iloc[:1000]
data_train, data_test = train_test_split(data_transformed,
                                         test_size=0.2)

responses = [[6,8],[2,5],[3,4],[1,7]]
logr_models = [{}] * len(responses)
i = 0
data_train_current = data_train.copy()

for response in responses:
  state_train = data_train_current.copy()
  if response != [3,4]:
    state_train['Response'] = state_train['Response'].apply(lambda x: 0 if x not in response else x)
  kappa_score,columns_opt = find_opt_params(state_train)
  print(i,kappa_score,columns_opt)
  logr_models[i] = {
      'model': regression_model(state_train,columns_opt),
      'columns': columns_opt
  }
  if response != [3,4]:
    data_train_current = data_train_current[~data_train_current['Response'].isin(response)]
  i += 1

logr_models[1]

def logr_hierarchy(x):
  for response in range(len(responses)):
    if x['target' + str(response)] > 0:
      x['target'] = x['target' + str(response)]
      break
  return x

for response in range(len(responses)):
  model = logr_models[response]['model']
  columns = logr_models[response]['columns']
  x = pd.DataFrame(data_test,columns=columns)
  data_test['target' + str(response)] = model.predict(x)

data_test = data_test.apply(logr_hierarchy,axis=1,result_type='expand')

data_test.head()

cohen_kappa_score(data_test['target'],data_test['Response'])