text:
'''
Created on Jun 12, 2017

@author: vital
'''

#!/usr/bin/python
import sys
import numpy as np
from numpy import float_, absolute as abs, random as ran

import time
import ssp_Hector_Fit3D_tools as ssp
import pyfits as pyf
from pyfits import getheader as ghead, getdata as gdata, writeto as wfits
from scipy.interpolate.interpolate import interp1d
import ssp_Hector_Fit3D_my as my
import os.path as pt
import matplotlib
import os

#Example command
#
'''
cd workspace/dazer/bin/lib/ssp_functions/

Pyhon command to run original file
python ssp_Hector_Fit3D.py NGC5947.spec_5.txt ssp_lib.fits,ssp_lib.3.fits auto_ssp.NGC5947.cen.only.out mask_elines.txt auto_ssp_V500_several_Hb.config 1 -1 40 3850 6800 emission_lines.txt 0.02 0.001 0.015 0.025 2 0.5 1 9 0.5 0.1 0.0 1.6
'''

def sycall(comand):
    from subprocess import call
    line=comand.split(" ")
    fcomand=[]
    fcomand.extend(line)
    linp=''
    nx=len(fcomand)
    for i in range(0, nx):
        linp=linp+fcomand[i]+" "
    print linp
    call(fcomand)

sys.argv=filter(None,sys.argv)
ran.seed(None)
vel_light=299792.458
red_elines=0.0
sec_ini=ssp.print_time()
time1=ssp.get_time()
if len(sys.argv) < 7:
    print "USE: auto_ssp.py SPEC1.txt SSP_SFH.fits,SSP_KIN.fits OUTFILE MASK_LIST CONFIG_FILE PLOT [min max] [wmin wmax] [redshift_elines_to_mask] [input_redshift delta_redshift min_redshift max_redshift] [input_sigma delta_sigma min_sigma max_sigma] [input_Av delta_Av min_Av max_Av]"
    print "CONFIG_FILE:"
    print "redshift delta_redshift min_redshift max_redshift"
    print "sigma delta_sigma min_sigma max_sigma"
    print "Av delta_Av min_Av max_Av [Same range for all]"
    print "N_SYSTEMS"
    print "(1) START_W END_W MASK_FILE CONFIG_FILE NPOLY MASK_FILE_POLY"
    print "..."
    print "(N) START_W END_W MASK_FILE CONFIG_FILE NPOLY MASK_FILE_POLY"
    print "MIN_DELTA_CHISQ MAX_NITER CUT_MEDIAN_FLUX"
    print "start_w_peak end_w_peak"
    print "wavelength_to_norm width_AA new_back_templates.fits"
    inline_params = ['ssp_Hector_Fit3D.py', 'NGC5947.spec_5.txt','ssp_lib.fits,ssp_lib.3.fits','auto_ssp.NGC5947.cen.only.out','mask_elines.txt','auto_ssp_V500_several_Hb.config' ,'1', '-1', '40', '3850', '6800', 'emission_lines.txt', '0.02', '0.001', '0.015', '0.025', '2', '0.5', '1', '9', '0.5', '0.1', '0.0', '1.6']
    sys.argv = inline_params
unc_file=sys.argv[1]
print 'This file', unc_file, os.getcwd()

clean_file="clean_"+sys.argv[1]
junk_back_list=sys.argv[2]
data=junk_back_list.split(',')
if len(data) == 2:
    back_list=data[0]
    back_list2=data[1]
else:
    back_list=junk_back_list
    back_list2=junk_back_list
outfile=sys.argv[3]
out_elines="elines_"+outfile
out_single="single_"+outfile
out_fit="fit_"+outfile
out_coeffs_file="coeffs_"+outfile
out_fit="output."+outfile+".fits"
out_ps=outfile
#######################################
# Clean previous results
#######################################
call="rm -rf "+outfile
sycall(call)
call="rm -rf "+out_elines
sycall(call)
call="rm -rf "+out_single
sycall(call)
call="rm -rf "+out_fit
sycall(call)
D_SYS_VEL=100
mask_list=sys.argv[4]
config_file=sys.argv[5]
plot=int(sys.argv[6])
if plot == 2:
    matplotlib.use('Agg')
    import matplotlib.pyplot as plt
    import matplotlib.colors as colors
    import matplotlib.cm as cmx
    dev_plot=outfile+".pdf"
    dev_plot_single="single_"+outfile+".pdf"
else:
    if plot == 0:
        matplotlib.use('Agg')
    dev_plot="null"
    dev_plot_single="null"
    import matplotlib.pyplot as plt
    import matplotlib.colors as colors
    import matplotlib.cm as cmx
smooth=1
MIN_CHISQ=1e12
out_file="junk.junk"
factor=1
box=1
deft=0
if len(sys.argv) == 9:
    min=float_(sys.argv[7])
    max=float_(sys.argv[8])
    deft=1
if len(sys.argv) == 11:
    min=float_(sys.argv[7])
    max=float_(sys.argv[8])
    min_wave=sys.argv[9]
    max_wave=sys.argv[10]
    deft=2
if len(sys.argv) == 12:
    min=float_(sys.argv[7])
    max=float_(sys.argv[8])
    min_wave=sys.argv[9]
    max_wave=sys.argv[10]
    elines_mask=sys.argv[11]
    deft=2
input_redshift=0
if len(sys.argv) == 13:
    min=float_(sys.argv[7])
    max=float_(sys.argv[8])
    min_wave=sys.argv[9]
    max_wave=sys.argv[10]
    elines_mask=sys.argv[11]
    input_redshift=float_(sys.argv[12])
    redshift=input_redshift
    deft=2
f=open(config_file,'r')
line=f.readline()
data=line.split(" ")
data=filter(None,data)
redshift=float_(data[0])
d_redshift=float_(data[1])
min_redshift=float_(data[2])
max_redshift=float_(data[3])
DV=float_(data[4])
RV=float_(data[5])
DS=float_(data[6])
RS=float_(data[7])
MIN_W=float_(data[8])
MAX_W=float_(data[9])
if len(sys.argv) == 16:
    min=float_(sys.argv[7])
    max=float_(sys.argv[8])
    min_wave=sys.argv[9]
    max_wave=sys.argv[10]
    elines_mask=sys.argv[11]
    input_redshift=float_(sys.argv[12])
    input_d_redshift=float_(sys.argv[13])
    input_min_redshift=float_(sys.argv[14])
    input_max_redshift=float_(sys.argv[15])
    redshift=input_redshift
    d_redshift=input_d_redshift
    min_redshift=input_min_redshift
    max_redshift=input_max_redshift
    deft=2
line=f.readline()
data=line.split(" ")
data=filter(None,data)
sigma=data[0]
d_sigma=data[1]
min_sigma=data[2]
max_sigma=data[3]
if len(sys.argv) == 20:
    min=float_(sys.argv[7])
    max=float_(sys.argv[8])
    min_wave=sys.argv[9]
    max_wave=sys.argv[10]
    elines_mask=sys.argv[11]
    input_redshift=float_(sys.argv[12])
    input_d_redshift=float_(sys.argv[13])
    input_min_redshift=float_(sys.argv[14])
    input_max_redshift=float_(sys.argv[15])
    sigma=float_(sys.argv[16])
    d_sigma=float_(sys.argv[17])
    min_sigma=float_(sys.argv[18])
    max_sigma=float_(sys.argv[19])
    redshift=input_redshift
    d_redshift=input_d_redshift
    min_redshift=input_min_redshift
    max_redshift=input_max_redshift
    deft=2
line=f.readline()
data=line.split(" ")
data=filter(None,data)
Av_IN=data[0]
d_Av_IN=data[1]
min_Av=data[2]
max_Av=data[3]
if len(sys.argv) == 24:
    min=float_(sys.argv[7])
    max=float_(sys.argv[8])
    min_wave=sys.argv[9]
    max_wave=sys.argv[10]
    elines_mask=sys.argv[11]
    input_redshift=float_(sys.argv[12])
    input_d_redshift=float_(sys.argv[13])
    input_min_redshift=float_(sys.argv[14])
    input_max_redshift=float_(sys.argv[15])
    sigma=float_(sys.argv[16])
    d_sigma=float_(sys.argv[17])
    min_sigma=float_(sys.argv[18])
    max_sigma=float_(sys.argv[19])
    Av_IN=float_(sys.argv[20])
    d_Av_IN=float_(sys.argv[21])
    min_Av=float_(sys.argv[22])
    max_Av=float_(sys.argv[23])
    redshift=input_redshift
    d_redshift=input_d_redshift
    min_redshift=input_min_redshift
    max_redshift=input_max_redshift
    deft=2
data=min_wave.split(',')
if len(data) == 2:
    min_wave=float_(data[0])
    min_wave2=float_(data[1])
else:
    min_wave=float_(min_wave)
    min_wave2=min_wave
data=max_wave.split(',')
if len(data) == 2:
    max_wave=float_(data[0])
    max_wave2=float_(data[1])
else:
    max_wave=float_(max_wave)
    max_wave2=max_wave
REDSHIFT=redshift
Av_ini=Av_IN
if d_redshift !=0:
    fit_redshift=1
else:
    fit_redshift=0
print "FIT_RED "+str(fit_redshift)+" "+str(d_redshift)+" "+str(len(sys.argv))
line=f.readline()
data=line.split(" ")
data=filter(None,data)
ns=int(data[0])
start_w_min=1e12
end_w_max=-1e12
start_w_E=[]
end_w_E=[]
mask_E=[]
config_E=[]
n_line_E=[]
npoly_E=[]
mask_poly_E=[]
nmin_E=[]
nmax_E=[]
config_line_E=[]
for i in range(0, ns):
    line=f.readline()
    data=line.split(" ")
    data=filter(None,data)
    start_w_e=float_(data[0])
    end_w_e=float_(data[1])
    mask_e=data[2]
    config_e=data[3]
    npoly_e=int(data[4])
    mask_poly_e=data[5]
    nmin_e=float_(data[6])
    nmax_e=float_(data[7])
    start_w_E.extend([start_w_e])
    end_w_E.extend([end_w_e])
    mask_E.extend([mask_e])
    config_E.extend([config_e])
#
# We read all the information
#
    n_line=0
    linef=""
    f2=open(config_e,'r')
    for line in f2:
        linef=linef+line+";"
        n_line=n_line+1
    config_line_E.extend([linef])
    f2.close()
    n_line_E.extend([n_line])
    npoly_E.extend([npoly_e])
    mask_poly_E.extend([mask_poly_e])
    nmin_E.extend([nmin_e])
    nmax_E.extend([nmax_e])
    if start_w_e < start_w_min:
        start_w_min=start_w_e
    if end_w_e > end_w_max:
        end_w_max=end_w_e
line=f.readline()
data=line.split(" ")
data=filter(None,data)
MIN_DELTA_CHISQ=float_(data[0])
MAX_NITER=int(data[1])
CUT_MEDIAN_FLUX=float_(data[2])
ABS_MIN=0.5*CUT_MEDIAN_FLUX
line=f.readline()
data=line.split(" ")
data=filter(None,data)
start_w_peak=float_(data[0])
end_w_peak=float_(data[1])
line=f.readline()
data=line.split(" ")
data=filter(None,data)
if len(data) == 3:
    wave_norm=data[0]
    w_wave_norm=data[1]
    new_back_file=data[2]
else:
    wave_norm=[]
    w_wave_norm=[]
    new_back_file=[]
f.close()
#
# SFH template
#
[pdl_flux_c_ini,hdr]=gdata(back_list, 0, header=True)
[nf,n_c]=pdl_flux_c_ini.shape
coeffs=np.zeros([nf,3])
crpix=hdr["CRPIX1"]
cdelt=hdr["CDELT1"]
crval=hdr["CRVAL1"]
n_mc=30
#
# Kinematics template
#
[pdl_flux_c_ini2,hdr2]=gdata(back_list2, 0, header=True)
[nf2,n_c2]=pdl_flux_c_ini2.shape
coeffs2=np.zeros([nf2,3])
crpix2=hdr2["CRPIX1"]
cdelt2=hdr2["CDELT1"]
crval2=hdr2["CRVAL1"]
Av=np.zeros(nf)
d_Av=np.zeros(nf)
for i in range(0, nf):
    Av[i]=Av_IN
    d_Av[i]=d_Av_IN
if mask_list == "none":
    nmask=0
else:
    f=open(mask_list,'r')
    start_mask=[]
    end_mask=[]
    for line in f:
        data=line.split(" ")
        data=filter(None,data)
        if len(data) != 0 and data[0] != "\n":
            start_mask.extend([float_(data[0])])
            end_mask.extend([float_(data[1])])
    nmask=len(start_mask)
    f.close()
    n_mask_org=nmask
if elines_mask == "none":
    nmask_e=0
    nline=0
else:
    f=open(elines_mask,'r')
    nmask_e=0
    nline=0
    w_eline=[]
    start_mask_e=[]
    end_mask_e=[]
    for line in f:
        data=line.split(" ")
        data=filter(None,data)
        if data[0] != "#":
            w_eline.extend([float_(data[0])])
            start_mask_e.extend([w_eline[nline]*(1+input_redshift)-4*sigma])
            end_mask_e.extend([w_eline[nline]*(1+input_redshift)+4*sigma])
            nmask_e=nmask_e+1
            nline=nline+1
    f.close()
#
#  We read the input spectrum
#
n_unc=0
y_min=1e12
y_max=-1e12
f=open(unc_file,'r')
i_scale=0
FLUX=0
have_error=0
index_unc=[]
wave_unc=[]
flux_unc=[]
flux_unc_org=[]
flux_unc_input=[]
e_flux_unc=[]
color_unc=[]
masked=[]
masked2=[]
masked_Av=[]
flux_masked=[]
flux_masked2=[]
e_flux_unc_kin=[]
wave_scale=0
for line in f:
    data=line.split(' ')
    data=filter(None,data)
    if data[0] != "#":
        index_unc.extend([float_(data[0])])
        wave_unc.extend([float_(data[1])])
        flux_unc.extend([float_(data[2])])
        flux_unc_org.extend([float_(data[2])])
        flux_unc_input.extend([float_(data[2])])
        if len(data) > 2:
        # Variance Column!
            e_flux_unc.extend([np.sqrt(abs(float_(data[3])))])
            color_unc.extend([1])#$data[4];
            have_error=1
        else:
            e_flux_unc.extend([np.sqrt(abs(float_(data[2])))/10])
            color_unc.extend([1])
        if np.isnan(flux_unc[n_unc]):
            flux_unc[n_unc]=flux_unc[n_unc-1]
        if flux_unc[n_unc] < y_min:
            y_min=flux_unc[n_unc]
        if flux_unc[n_unc] > y_max:
            y_max=flux_unc[n_unc]
        if n_unc > 0:
            if wave_unc[n_unc-1] <= wave_scale and wave_unc[n_unc] > wave_scale:
                i_scale=n_unc
        masked.extend([1])
        masked2.extend([1])
        masked_Av.extend([1])
        if flux_unc[n_unc] == 0:
            masked[n_unc]=0
            masked2[n_unc]=0
        w_test=wave_unc[n_unc-1]
        for j in range(0, nmask):
            if w_test > start_mask[j] and w_test < end_mask[j]:
                masked[n_unc]=0
                masked2[n_unc]=0
                masked_Av[n_unc]=0
        if deft == 2:
            if w_test < min_wave:
                masked[n_unc]=0
                masked_Av[n_unc]=0
            if w_test > max_wave:
                masked[n_unc]=0
                masked_Av[n_unc]=0
            if w_test < min_wave2:
                masked2[n_unc]=0
            if w_test > max_wave2:
                masked2[n_unc]=0
        for j in range(0, nmask_e):
            if w_test > start_mask_e[j] and w_test < end_mask_e[j]:
                masked2[n_unc]=0
                masked_Av[n_unc]=0
        flux_masked.extend([flux_unc[n_unc]*masked[n_unc]])
        flux_masked2.extend([flux_unc[n_unc]*masked2[n_unc]])
        if wave_unc[n_unc] > min_wave and wave_unc[n_unc] < max_wave:
            FLUX=FLUX+flux_masked[n_unc]
        e_flux_unc_kin.extend([e_flux_unc[n_unc]])
        n_unc=n_unc+1
f.close()

sigma_e=np.median(e_flux_unc)
print "SIGMA_E = "+str(sigma_e)
for i in range(0, n_unc):
    if e_flux_unc[i] > 1.5*sigma_e:
        e_flux_unc[i]=1.5*sigma_e
        e_flux_unc_kin[i]=1.5*sigma_e
if deft == 2:
    y_min=min
    y_max=max
else:
    min_wave=np.amin(wave_unc)
    max_wave=np.amax(wave_unc)
if deft == 1:
    y_min=min
    y_max=max
median_flux=np.median(flux_masked)
dpix_unc=wave_unc[1]-wave_unc[0]
max=3*median_flux
pdl_output=np.zeros([6,n_unc])
#
# We create a kernel
#
med_flux=np.median(flux_unc)
chi_sq_min_now=1e12
min_chi_sq=chi_sq_min_now
print '-----The redshift', redshift
print '-----The sigma', sigma
print '-----The Av', Av.shape
print '-----The crval2', crval2
print '-----The cdelt2', cdelt2
print '-----The crpix2', crpix2
print '-----back_list2', back_list2
print '-----nf2', nf2
print '-----n_c2', n_c2
print '-----pdl_flux_c_ini2', pdl_flux_c_ini2.shape
#print '-----hdr2', hdr2
print '-----wave_unc', wave_unc
print '-----masked_Av', len(masked_Av)
print '-----e_flux_unc', e_flux_unc
print '-----flux_unc', flux_unc
print '-----n_mc', n_mc
print '-----chi_sq_min_now', chi_sq_min_now
print '-----min_chi_sq', min_chi_sq

ssp_dat=ssp.fit_ssp_lin_no_zero(redshift,sigma,Av,crval2,cdelt2,crpix2,nf2,n_c2,pdl_flux_c_ini2,hdr2,wave_unc,masked_Av,e_flux_unc,flux_unc,n_mc,chi_sq_min_now,min_chi_sq,plot)
min_chi_sq=ssp_dat[0]
# print "CUT = "+str(med_flux)+" "+str(CUT_MEDIAN_FLUX)
# print str(redshift)+","+str(sigma)

print 'Aqui acabamos'
sys.exit(0)

if med_flux < ABS_MIN:
# WHAT TO DO???
# We print all!!!
    sys.exit(0)
if med_flux > CUT_MEDIAN_FLUX:
    if MIN_W == 0:
        MIN_W = min_wave
    if MAX_W == 0:
        MAX_W=max_wave

################
# REDSHIFT DETERMINATION
    my_plot=2
    K=0
    nr=0
    chi_r=[]
    red_r=[]
    if d_redshift > 0:
        min_chi_sq=1e30
        RED=min_redshift
        while RED < max_redshift:
            ssp_dat1=ssp.fit_ssp_lin_no_zero_no_cont(RED,sigma,Av,crval2,cdelt2,crpix2,nf2,n_c2,pdl_flux_c_ini2,hdr2,wave_unc,masked2,e_flux_unc_kin,flux_masked2,n_mc,chi_sq_min_now,min_chi_sq,my_plot)
            chi_now=ssp_dat1[0]
            chi_r.extend([chi_now])
            red_r.extend([RED])
#            print RED,chi_now,d_redshift
            if nr > 1  and chi_r[nr-1] < min_chi_sq and chi_r[nr-1] > 0:
                redshift=red_r[nr-1]
                min_chi_sq=chi_r[nr-1]
                K=nr-1
            nr=nr+1
            RED=RED+d_redshift
        #
        # TWO
        #
        e_redshift=d_redshift
        nr=0
        chi_r=[]
        red_r=[]
        RED=redshift-1.5*d_redshift
        max_redshift=redshift+1.5*d_redshift
        d_redshift=0.1*d_redshift
        while RED < max_redshift:
            ssp_dat2=ssp.fit_ssp_lin_no_zero_no_cont(RED,sigma,Av,crval2,cdelt2,crpix2,nf2,n_c2,pdl_flux_c_ini2,hdr2,wave_unc,masked2,e_flux_unc,flux_masked2,n_mc,chi_sq_min_now,min_chi_sq,my_plot)
            chi_now=ssp_dat2[0]
            chi_r.extend([chi_now])
            red_r.extend([RED])
            if nr > 1 and chi_r[nr-1] < chi_r[nr-2] and chi_r[nr-1] < chi_r[nr] and chi_r[nr-1] <= min_chi_sq:
                a=red_r[nr-2]
                b=red_r[nr-1]
                c=red_r[nr]
                fa=chi_r[nr-2]
                fb=chi_r[nr-1]
                fc=chi_r[nr]
                den=(fc-2*fb+fa)
                redshift=c-(b-a)*((fc-fb)/den+0.5)
                slope=abs(0.5*(fc-fb)/(c-b))+abs(0.5*(fa-fb)/(a-b))
                if slope > 0:
                    e_redshift=0.01*redshift/slope
                else:
                    e_redshift=0.01*redshift
                min_chi_sq=chi_r[nr-1]
                K=nr-1
            nr=nr+1
            a_rnd=ran.rand(2)
            RED=RED+d_redshift*(a_rnd[0])
        fit_redshift=0
        d_redshift=0
    else:
        fit_redshift=0
        e_redshift=0
    print "REDSHIFT = "+str(redshift)+" +- "+str(e_redshift)
    #sys.exit()
    REDSHIFT=redshift
    # sigma DETERMINATION
    K=0
    nr=0
    chi_s=[]
    sigma_s=[]
    print "D_SIGMA = "+str(d_sigma)
    if d_sigma > 0:
        min_chi_sq = 1e30
        SIGMA=min_sigma
        while SIGMA < max_sigma:
            ssp_dat3=ssp.fit_ssp_lin_no_zero_no_cont(redshift,SIGMA,Av,crval2,cdelt2,crpix2,nf2,n_c2,pdl_flux_c_ini2,hdr2,wave_unc,masked2,e_flux_unc,flux_masked2,n_mc,chi_sq_min_now,min_chi_sq,my_plot)
            chi_now=ssp_dat3[0]
            chi_s.extend([chi_now])
            sigma_s.extend([SIGMA])
            if chi_s[nr-1] < min_chi_sq:
                sigma=sigma_s[nr-1]
                min_chi_sq=chi_s[nr-1]
                K=nr
            nr=nr+1
            SIGMA=SIGMA+d_sigma
        SIGMA=sigma-1.5*d_sigma
        max_sigma=sigma+1.5*d_sigma
        d_sigma=0.33*d_sigma
        #
        #TWO
        #
        nr=0
        chi_s=[]
        sigma_s=[]
        e_sigma=d_sigma
        while SIGMA < max_sigma:
            ssp_dat4=ssp.fit_ssp_lin_no_zero_no_cont(redshift,SIGMA,Av,crval2,cdelt2,crpix2,nf2,n_c2,pdl_flux_c_ini2,hdr2,wave_unc,masked2,e_flux_unc,flux_masked2,n_mc,chi_sq_min_now,min_chi_sq,my_plot)
            chi_now=ssp_dat4[0]
            chi_s.extend([chi_now])
            sigma_s.extend([SIGMA])
            if nr > 1 and chi_s[nr-1] < chi_s[nr-2] and chi_s[nr-1] < chi_s[nr] and chi_s[nr-1] <= min_chi_sq:
                a=sigma_s[nr-2]
                b=sigma_s[nr-1]
                c=sigma_s[nr]
                fa=chi_s[nr-2]
                fb=chi_s[nr-1]
                fc=chi_s[nr]
                den=(fc-2*fb+fa)
                sigma=c-(b-a)*((fc-fb)/den+0.5)
                min_chi_sq=chi_s[nr-1]
                K=nr
                SIGMA=max_sigma
            nr=nr+1
            a_rnd=ran.rand(2)
            SIGMA=SIGMA+d_sigma*(a_rnd[0])
        slope=(chi_s[nr-1]-min_chi_sq)/(sigma_s[nr-1]-sigma)
        if slope > 0:
            e_sigma=sigma/slope/10.
        else:
            e_sigma=sigma/10.
        fit_sigma=0
        d_sigma=0
    else:
        fit_sigma=0
        e_sigma=0
    sigma=abs(sigma)
    e_sigma=abs(e_sigma)
    print "SIGMA = "+str(sigma)+"+-"+str(e_sigma)
else:
#
# Below the cut!
#
    for i in range(0, nf):
        Av[i]=0
        d_Av[i]=0
# Av DETERMINATION
K=0
nr=0
chi_Av=[]
Av_s=[]
Av_p_chi=[]
print "D_Av = "+str(d_Av_IN)
nr_min=0
if d_Av_IN > 0:
    min_chi_sq = 1e30
    Av_NOW=min_Av
    while Av_NOW < max_Av:
        for i in range(0, nf):
            Av[i]=Av_NOW
        #
        # Not allow negative coeffs!!!!
        #
        ssp_dat5=ssp.fit_ssp_lin_no_zero(redshift,sigma,Av,crval2,cdelt2,crpix2,nf2,n_c2,pdl_flux_c_ini2,hdr2,wave_unc,masked_Av,e_flux_unc,flux_masked,n_mc,chi_sq_min_now,min_chi_sq,my_plot)
        chi_now=ssp_dat5[0]
        chi_Av.extend([chi_now])
        Av_s.extend([Av_NOW])
        if chi_now > 0:
            Av_p_chi.extend([Av_NOW/(chi_now)])
        if K == 0 and chi_Av[nr] < min_chi_sq:
            Av_F=Av_s[nr]
            nr_min=nr
            min_chi_sq=chi_now
        if nr > 1 and chi_Av[nr-1] < chi_Av[nr-2] and chi_Av[nr-1] < chi_Av[nr] and chi_Av[nr-1] <= min_chi_sq:
            a=Av_s[nr-2]
            b=Av_s[nr-1]
            c=Av_s[nr]
            fa=chi_Av[nr-2]
            fb=chi_Av[nr-1]
            fc=chi_Av[nr]
            den=(fc-2*fb+fa)
            Av_F=c-(b-a)*((fc-fb)/den+0.5)
            min_chi_sq=chi_Av[nr-1]
            K=nr
        nr=nr+1
        a_rnd=ran.rand(2)
        Av_NOW=Av_NOW+d_Av_IN*(a_rnd[0])
    if Av_s[nr-1] != Av_F:
        slope=(chi_Av[nr-1]-min_chi_sq)/(Av_s[nr-1]-Av_F)
        if slope > 0 :
            e_Av=abs(Av_F/slope/3.)
        else:
            e_Av=d_Av_IN
    else:
        e_Av=d_Av_IN
    fit_Av=0
    d_Av_NOW=0
else:
    fit_Av=0
if d_Av_IN == 0:
    Av_F=Av_IN
if e_Av == 0:
    e_Av=d_Av_IN
print "AV = "+str(Av_F)+" +- "+str(e_Av)
for i in range(0, nf):
    Av[i]=Av_F
fit_redshift=0
redshift_abs=redshift
delta_chi=10
NITER=0
niter_tmp_max=10
chi_sq_min_now=1e12
min_chi_sq_limit=min_chi_sq
n_mc=10
pdl_rat_master=np.ones(n_unc+1)
[min_chi_sq,pdl_age_mod,pdl_met_mod,pdl_ml,pdl_Av,coeffs,coeffs_N,coeffs_NM,pdl_model_spec_min,pdl_res]=ssp.fit_ssp_lin_no_zero(redshift,sigma,Av,crval2,cdelt2,crpix2,nf2,n_c2,pdl_flux_c_ini2,hdr2,wave_unc,masked_Av,e_flux_unc,flux_unc,n_mc,chi_sq_min_now,min_chi_sq,plot)
#
# We substract the continuum!
#
pdl_mod_SSP=pdl_model_spec_min
pdl_res_SSP=pdl_res
nx=n_unc
i0_now=int(0.4*nx)
i1_now=int(0.6*nx)
stats_res=np.std(pdl_res[i0_now:i1_now])+np.mean(pdl_res[i0_now:i1_now])
stats_mod=np.mean(pdl_model_spec_min[i0_now:i1_now])
SN=0
if stats_res > 0:
    SN=stats_mod/stats_res
print "Signal-to-Noise = "+str(SN)
old=1
if old == 1 and SN > 10:
    pdl_model_spec_min[np.where(pdl_model_spec_min == 0)[0]]=1.
    pdl_rat=pdl_res/pdl_model_spec_min+1
    rat=pdl_rat
    med_rat=my.median_filter(int(5*2.354*sigma),rat)
    pdl_med_rat=np.array(med_rat)
    n_unc_1=n_unc-1
    pdl_wave_unc=wave_unc[0]+(wave_unc[1]-wave_unc[0])*np.arange(0,n_unc_1)
    med_rat=my.median_filter(int(7*2.354*sigma),rat)
    med_sigma=int(1.5*sigma)
    if med_sigma < 3:
        med_sigma=3
    med_rat_box=my.median_box(med_sigma,med_rat)
    med_wave_box=my.median_box(med_sigma,wave_unc)
    y_rat = interp1d(med_wave_box, med_rat_box,bounds_error=False,fill_value=0.)(wave_unc)
    if plot > 0:
        out_ps_now="junk2"
        title="ratio"
        ssp.plot_results_min_max(2,wave_unc,[flux_unc,pdl_model_spec_min,pdl_res,pdl_rat,y_rat],out_ps_now,title,-0.2,1.5)
    i0_now=int(0.4*n_unc)
    i1_now=int(0.6*n_unc)
    stats_rat0=np.mean(y_rat[i0_now:i1_now])
    stats_rat1=np.std(y_rat[i0_now:i1_now])+stats_rat0
    if stats_rat0 > 0 and stats_rat1 > 0.02:
        for i in range(0, n_unc):
            val=y_rat[i]
            if val > 0:
                flux_unc[i]=flux_unc[i]/val
print "Deriving SFH...."
#my pdl_mod_JOINT;
#my pdl_res_JOINT;#=$pdl_res_SSP; #Modificacion 12 de Marzo de 2015 (en caso de no entrar en el ciclo while)
#my pdl_no_gas;
#my age_min;
#my met_min;
#my Av_min;
#my age_min_mass;
#my met_min_mass;
#my Av_min_mass;
coeffs_cat=np.zeros([nf+1,n_mc])
while MIN_CHISQ > MIN_DELTA_CHISQ and NITER < MAX_NITER:
    if NITER == 1:
        MIN_CHISQ=1e12
######################################################################
# Fitting the emission lines
######################################################################
    a_fixed=np.zeros([1,9])
    a_type_fixed=[]
    n_mod_fixed=0
    if ns > 0:
        ks=0
        SYS_VEL=vel_light*REDSHIFT
        REN=[]
        e_REN=[]
        sycall(call)
        for ist in range(0,ns):
            if red_elines > 0:
                SYS_VEL=vel_light*red_elines
                if ist == 0:
                    SYS_VEL_MAX=vel_light*red_elines+D_SYS_VEL
                    SYS_VEL_MIN=vel_light*red_elines-D_SYS_VEL
                else:
                    SYS_VEL_MAX=vel_light*red_elines+D_SYS_VEL
                    SYS_VEL_MIN=vel_light*red_elines-D_SYS_VEL
            else:
                SYS_VEL=vel_light*REDSHIFT
                if ist == 0:
                    SYS_VEL_MAX=vel_light*REDSHIFT+D_SYS_VEL
                    SYS_VEL_MIN=vel_light*REDSHIFT-D_SYS_VEL
                else:
                    SYS_VEL_MAX=vel_light*REDSHIFT+D_SYS_VEL
                    SYS_VEL_MIN=vel_light*REDSHIFT-D_SYS_VEL
            start_w_e=start_w_E[ist]
            end_w_e=end_w_E[ist]
            mask_e=mask_E[ist]
            config_e=config_E[ist]
            npoly_e=npoly_E[ist]
            mask_poly_e=mask_poly_E[ist]
            nmin_e=nmin_E[ist]
            nmax_e=nmax_E[ist]
            print "CONF="+config_e
            wave_elines=[]
            flux_elines=[]
            flux_e_elines=[]
            masked_elines=[]
            n_e=0
            for i in range(0, n_unc):
                if wave_unc[i] > start_w_e and wave_unc[i] < end_w_e:
                    wave_elines.extend([wave_unc[i]])
                    flux_elines.extend([flux_unc_org[i]-pdl_mod_SSP[i]])
                    flux_e_elines.extend([abs(e_flux_unc[i])])
                    masked_elines.extend([1])
                    n_e=n_e+1
            pdl_wave_elines=np.array(wave_elines)
            pdl_flux_elines=np.array(flux_elines)
            pdl_flux_e_elines=np.array(flux_e_elines)
            pdl_masked_elines=np.array(masked_elines)
            stats0=np.mean(pdl_flux_elines)
            stats4=np.amax(pdl_flux_elines)
            y_max=stats4-stats0
            deft=1
            data=filter(None, config_line_E[ist].split(';')[0].split(" "))
            #print float_(filter(None, config_line_E[0].split(';')[4].split(" ")))[2]
            junk=data[0]
            n_mod=int(data[1])
            chi_goal=float_(data[2])
            d_chi_goal=float_(data[3])
            n_line=n_line_E[ist]
            i_mod=1
            typef=[]
            a=np.zeros([n_mod,9])
            ia=np.zeros([n_mod,9])
            ea=np.zeros([n_mod,9])
            a0=np.zeros([n_mod,9])
            a1=np.zeros([n_mod,9])
            link=np.zeros([n_mod,9])
            for ii in range(0, n_mod):
                cnf=filter(None, config_line_E[ist].split(';')[i_mod].split(" "))
                i_mod=i_mod+1
                typef.extend(cnf)
                for j in range(0, 9):
                    data=config_line_E[ist].split(';')[i_mod].replace('\t',' ')
                    data=filter(None, data.split(' '))
                    i_mod=i_mod+1
                    a[ii][j]=float_(data[0])
                    ia[ii][j]=float_(data[1])
                    ea[ii][j]=0
                    a0[ii][j]=float_(data[2])
                    a1[ii][j]=float_(data[3])
                    link[ii][j]=float_(data[4])
                if deft == 1:
                    a1_max=2*y_max*(a[ii][2]*((2*3.1416)**0.5))
                    a0_min=0.01*1.2*y_max*(a[ii][2]*((2*3.1416)**0.5))
                    if a1[ii][1] > a1_max:
                        a1[ii][1]=a1_max
                        a0[ii][1]=a0_min
            #
            # Force vicitiny in the velocity
            #
            a[0][3]=SYS_VEL
            ia[0][3]=1
            a0[0][3]=SYS_VEL_MIN
            a1[0][3]=SYS_VEL_MAX
            i_ter=0
            n_loops=5
            n_mc_e=30
            breakt=0
            scale_ini=0.15;
            deft=0;
            pdl_model=np.zeros(n_e)
            pdl_model_cont=np.zeros(n_e)
            pdl_model_tmp=np.zeros(n_e)
            pdl_model_cont_tmp=np.zeros(n_e)
            a_out_now=ssp.copy_a(n_mod,a)
            a_out_tmp=ssp.copy_a(n_mod,a)
            chi_sq_now=1e12
            a_results=np.zeros([1, n_mod, 9])
            nnk=0
            max_time=5
            time=ssp.get_seconds()
            d_time=ssp.get_seconds()-time
            i_loops=0
            ran.seed(None)
            while i_ter < n_loops and breakt == 0:
                chi_iter=chi_sq_now
                chi_single=0
                [chi_sq_now,pdl_a,pdl_model_tmp,pdl_model_cont_tmp]=ssp.fit_elines_grad_rnd_new(pdl_wave_elines,pdl_flux_elines,pdl_flux_e_elines,n_mod,chi_goal,d_chi_goal,typef,a_out_tmp,ia,a0,a1,link,n_mc_e,pdl_masked_elines,deft,scale_ini)#,max_time)
                a_out_now=ssp.copy_a_pdl(n_mod,pdl_a)
                #print chi_sq_now, pdl_a[:,1],a_out_tmp[:,1]
                if chi_sq_now < chi_iter:
#####################################################
                # Close to a result, narrow the range
                    for i in range(0, n_mod):
                        for j in range(0, 9):
                            if typef[i] == "eline\n":
                                if ia[i][j] == 1:
                                    if link[i][j] == -1:
                                        delta_now=abs(a1[i][j]-a0[i][j])/(2.)
                                        a0_tmp=a0[i][j]
                                        a1_tmp=a1[i][j]
                                        if j != 3:
                                            a0_tmp=a_out_now[i][j]-delta_now
                                            a1_tmp=a_out_now[i][j]+delta_now
                                        else:
                                            a0_tmp=a_out_now[i][j]-0.5*delta_now
                                            a1_tmp=a_out_now[i][j]+0.5*delta_now
                                        if a0_tmp < a0[i][j]:
                                            a0_tmp=a0[i][j]
                                        if a1_tmp > a1[i][j]:
                                            a1_tmp=a1[i][j]
                                        a0[i][j]=a0_tmp
                                        a1[i][j]=a1_tmp
####################################################
                    a_out_tmp=ssp.copy_a(n_mod,a_out_now)
                    a_results=ssp.copy_a_results(n_mod,nnk,a_out_now,a_results)
                    pdl_model=pdl_model_tmp
                    pdl_model_cont=pdl_model_cont_tmp
                    nnk=nnk+1
                    i_ter=i_ter+1
                else:
                    rnd_a=ran.rand(10);
                    a_out_now=ssp.copy_a(n_mod,a_out_now)
                    i_loops=i_loops+1
                    if i_loops > 5*n_loops:
                        breakt=1
            out_ps_now="fit_"+outfile+"."+str(start_w_e)+"_"+str(end_w_e)
            title="["+str(start_w_e)+","+str(end_w_e)+"]"
            if pdl_model.shape[0] == len(pdl_wave_elines):
                pdl_model=np.transpose(pdl_model)
            ssp.plot_results(plot,pdl_wave_elines,[pdl_flux_elines,pdl_model[0,:],(pdl_flux_elines-pdl_model[0,:])],out_ps_now,title)
            print "----------------------------------------";
            a_final=ssp.mean_a_results_last(n_mod,nnk,a_results,ia)
            #
            # Background noise
            #
            pdl_res_now=pdl_flux_elines-pdl_model
            stats_back1=np.mean(pdl_res_now)+np.std(pdl_res_now)
            a_final=ssp.add_back_noise(n_mod,a_final,typef,chi_sq_now,stats_back1)
            ssp.print_a_final(n_mod,a_final,typef,chi_sq_now)
            out_fit_spectra=out_elines
            ssp.print_a_final_file_add(n_mod,a_final,typef,chi_sq_now,out_fit_spectra)
            [n_mod_fixed,junk_a_fixed,junk_a_type_fixed]=ssp.add_a_results_elines(n_mod,a_final,typef,n_mod_fixed,a_fixed,a_type_fixed)
            a_fixed=junk_a_fixed
            a_type_fixed=junk_a_type_fixed
            nmin_e=int(0.1*n_unc)
            nmax_e=int(0.9*n_unc)
###############################
            # Low order polynomical!
            out_fit_now=out_fit+"."+str(start_w_e)+"_"+str(end_w_e)+".pdf"
            box=int(sigma*6)
            print "DONE FIT ELINES CONFIG "+str(ist)
        #
        # We create a FIXED model of the emission lines
        #
        pdl_model_elines=np.zeros(n_unc)
        pdl_model_cont=np.zeros(n_unc)
        pdl_wave_elines=np.array(wave_unc)
        NN=len(pdl_wave_elines)
        NN1=len(pdl_model_elines)
        for i in range(0, n_mod_fixed):
            pdl_tmp=ssp.create_single_model(pdl_wave_elines,i,a_type_fixed,a_fixed)
            NN2=len(pdl_tmp[0,:])
            pdl_model_elines=pdl_model_elines+pdl_tmp[0,:]
        #
        # We remove the gas before a new iteration
        #
        for i in range(0, n_unc):
            flux_unc[i]=flux_unc_org[i]-pdl_model_elines[i]
        pdl_mod_JOINT=pdl_mod_SSP+pdl_model_elines
        pdl_res_JOINT=pdl_res_SSP-pdl_model_elines
        pdl_no_gas=np.array(flux_unc)
        #############################################################
        # We rescale!
        ##############################################################
        y_rat=np.ones(nx+1)
        jy=0
        if SN > 10:
            pdl_mod_JOINT[np.where(pdl_mod_JOINT == 0)[0]]=1.
            pdl_rat=pdl_res_JOINT/pdl_mod_JOINT+1
            rat=pdl_rat
            n_unc_1=n_unc-1
            pdl_wave_unc=wave_unc[0]+(wave_unc[1]-wave_unc[0])*np.arange(0, n_unc_1)
            med_rat=my.median_filter(int(5*2.354*sigma),rat);
            med_sigma=int(1.5*sigma)
            if med_sigma < 3:
                med_sigma=3
            med_rat_box=my.median_box(med_sigma,med_rat)
            med_wave_box=my.median_box(med_sigma,wave_unc)
            y_rat = interp1d(med_wave_box, med_rat_box,bounds_error=False,fill_value=0.)(wave_unc)
            i0_now=int(0.4*nx)
            i1_now=int(0.6*nx)
            stats_rat0=np.mean(y_rat[i0_now:i1_now])
            stats_rat1=np.std(y_rat[i0_now:i1_now])+stats_rat0
            if plot > 1:
                out_ps_now="junk3"
                title="ratio = "+str(stats_rat0)+", rms="+str(stats_rat1)
                print title
                ssp.plot_results_min_max(2,wave_unc,[flux_unc,pdl_model_spec_min,pdl_res,pdl_rat,y_rat],out_ps_now,title,-0.2,1.5)
            if stats_rat0 > 0 and stats_rat1 > 0.02:
                if jy == 0:
                    # Continuum shape correction on/off
                    pdl_rat_master=y_rat
                    pdl_rat_master[np.where(pdl_rat_master == 0)[0]]=1.
                    y_rat=pdl_rat_master
                else:
                    y_rat=pdl_rat_master
        for i in range(0, n_unc):
            val=y_rat[i]
            if val > 0:
                flux_unc[i]=flux_unc[i]/val
                flux_unc_org[i]=flux_unc_org[i]/val
        ##############################################################
        # End re-scale
        ##############################################################
        ML=0
        if med_flux > CUT_MEDIAN_FLUX:
            n_mc=20
            [min_chi_sq,pdl_age_mod,pdl_met_mod,pdl_ml,pdl_Av,coeffs,coeffs_N,coeffs_NM,pdl_mod_SSP,pdl_res_SSP,coeffs_N_input,e_coeffs_N_input]=ssp.fit_ssp_lin_MC(redshift,sigma,Av,crval,cdelt,crpix,nf,n_c,pdl_flux_c_ini,hdr,wave_unc,masked,e_flux_unc,flux_unc,n_mc,chi_sq_min_now,MIN_CHISQ,plot)
            smooth_ratiot=ssp.smooth_ratio(flux_unc,pdl_mod_SSP,int(sigma))
            pdl_mod_SSP_no_cor=np.copy(pdl_mod_SSP)
            pdl_mod_SSP=pdl_mod_SSP*smooth_ratiot
            f1=open(out_coeffs_file, "w")
            f1.write("# ID   AGE     MET    COEFF   Min.Coeff  log(M/L)   AV   N.Coeff   Err.Coeff\n")
            print "------------------------------------------------------------------------------"
            print "ID   AGE     MET    COEFF   Min.Coeff  log(M/L)   AV   N.Coeff   Err.Coeff"
            print "------------------------------------------------------------------------------"
            age_mod=pdl_age_mod
            met_mod=pdl_met_mod
            Av_mod=pdl_Av
            ml=pdl_ml
            a_coeffs=coeffs[:,0]
            a_e_coeffs=coeffs[:,1]
            a_min_coeffs=coeffs[:,2]
            a_coeffs_N=coeffs_N
            a_e_coeffs_N=a_e_coeffs
            l_age_min=0
            l_met_min=0
            l_Av_min=0
            l_age_min_mass=0
            l_met_min_mass=0
            l_Av_min_mass=0
            e_l_age_min=0
            e_l_met_min=0
            e_l_Av_min=0
            e_l_age_min_mass=0
            e_l_met_min_mass=0
            e_l_Av_min_mass=0
            for k in range(0, nf):
                if a_coeffs[k] > 0:
                    a_e_coeffs_N[k]=a_e_coeffs[k]*(a_coeffs_N[k]/a_coeffs[k])
                else:
                    a_e_coeffs_N[k]=0
                l_age_min=l_age_min+a_coeffs[k]*np.log10(age_mod[k])
                l_met_min=l_met_min+a_coeffs[k]*np.log10(met_mod[k])
                l_Av_min=l_Av_min+a_coeffs[k]*np.log10(Av_mod[k])
                l_age_min_mass=l_age_min_mass+ml[k]*a_coeffs_N[k]*np.log10(age_mod[k])
                l_met_min_mass=l_met_min_mass+ml[k]*a_coeffs_N[k]*np.log10(met_mod[k])
                l_Av_min_mass=l_Av_min_mass+ml[k]*a_coeffs_N[k]*np.log10(Av_mod[k])
                e_l_age_min=e_l_age_min+a_e_coeffs[k]*np.log10(age_mod[k])
                e_l_met_min=e_l_met_min+a_e_coeffs[k]*np.log10(met_mod[k])
                e_l_Av_min=e_l_Av_min+a_e_coeffs[k]*np.log10(Av_mod[k])
                e_l_age_min_mass=e_l_age_min_mass+ml[k]*a_e_coeffs_N[k]*np.log10(age_mod[k])
                e_l_met_min_mass=e_l_met_min_mass+ml[k]*a_e_coeffs_N[k]*np.log10(met_mod[k])
                e_l_Av_min_mass=e_l_Av_min_mass+ml[k]*a_e_coeffs_N[k]*np.log10(Av_mod[k])
                ML=ML+ml[k]*a_coeffs_N[k]
                C_ini=coeffs_N_input[k]
                e_C_ini=e_coeffs_N_input[k]
                f1.write(("%2d" % k)+" "+("%7.4f" % age_mod[k])+" "+("%7.4f" % met_mod[k])+" "+("%7.4f" % a_coeffs_N[k])+" "+("%7.4f" % a_min_coeffs[k])+" "+("%4.4f" % np.log10(ml[k]))+" "+("%4.2f" % Av_mod[k])+" "+("%7.4f" % a_coeffs[k])+" "+("%7.4f" % a_e_coeffs[k])+"\n")
                if a_coeffs[k] > 1e-5:
                    print ("%2d" % k)+" "+("%7.4f" % age_mod[k])+" "+("%7.4f" % met_mod[k])+" "+("%7.4f" % a_coeffs_N[k])+" "+("%7.4f" % a_min_coeffs[k])+" "+("%4.4f" % np.log10(ml[k]))+" "+("%4.2f" % Av_mod[k])+" "+("%7.4f" % a_coeffs[k])+" "+("%7.4f" % a_e_coeffs[k])+" "+("%7.4f" % C_ini)+" "+("%7.4f" % e_C_ini)
            print "------------------------------------------------------------------------------"
            f1.close
            age_min=10**(l_age_min)
            met_min=10**(l_met_min)
            Av_min=10**(l_Av_min)
            age_min_mass=10**(l_age_min_mass/ML)
            met_min_mass=10**(l_met_min_mass/ML)
            Av_min_mass=10**(l_Av_min_mass/ML)
            e_age_min=abs(0.43*e_l_age_min*age_min)
            e_met_min=abs(0.43*e_l_met_min*met_min)
            e_Av_min=abs(0.43*e_l_Av_min*Av_min)
            e_age_min_mass=abs(0.43*e_l_age_min*age_min_mass)
            e_met_min_mass=abs(0.43*e_l_met_min*met_min_mass)
            e_Av_min_mass=abs(0.43*e_l_Av_min*Av_min_mass)
            if min_chi_sq > 0:
                delta_chi=abs((chi_sq_min_now-min_chi_sq)/min_chi_sq)
        wpeak=6562
        Fpeak=-1e12
        pdl_mod_JOINT=pdl_mod_SSP+pdl_model_elines
        pdl_res_JOINT=pdl_res_SSP-pdl_model_elines
        pdl_no_gas=np.array(flux_unc)
        # Copy output!
        pdl_output[0,:]=np.array(flux_unc_org)
        pdl_output[1,:]=pdl_mod_SSP
        pdl_output[2,:]=pdl_mod_JOINT
        pdl_res_SSP=np.array(flux_unc_org)-pdl_mod_SSP
        pdl_res_SSP_no_cor=np.array(flux_unc_input)-pdl_mod_SSP_no_cor
        pdl_output[3,:]=pdl_res_SSP_no_cor
        pdl_tmp=np.array(flux_unc_org)
        nx_1=n_unc#-1
        if len(pdl_rat_master)-len(pdl_mod_JOINT)==1:
            pdl_res_JOINT=pdl_tmp/(pdl_rat_master[0:nx_1])-pdl_mod_JOINT
        else:
            pdl_res_JOINT=pdl_tmp/(pdl_rat_master)-pdl_mod_JOINT
        pdl_output[4,:]=pdl_res_JOINT
        pdl_output[5,:]=np.array(flux_unc_org)-(pdl_res_SSP-pdl_res_JOINT)
        title="X="+str(chi_sq_now)+" T="+str(age_min)+" ("+str(age_min_mass)+") Z="+str(met_min)+" ("+str(met_min_mass)+") Av="+str(Av_min)+" z="+str(redshift)+" sigma="+str(sigma)
        ssp.plot_results(plot,pdl_wave_elines,pdl_output,out_ps,title)
    print "I.Iter = "+str(NITER)+" DONE"
    NITER=NITER+1
# Write output file
h=pyf.PrimaryHDU().header
h["NAXIS"]=2
h["NAXIS1"]=n_unc
h["NAXIS2"]=6
h["COMMENT"]="OUTPUT auto_ssp_elines_rnd.pl FITs"
h["CRVAL1"]=wave_unc[0]
h["CDELT1"]=wave_unc[1]-wave_unc[0];
h["CRPIX1"]=1
if pt.exists(out_fit) == False:
    wfits(out_fit,pdl_output,h)
else:
    sycall("rm "+out_fit)
    wfits(out_fit,pdl_output,h)
################################
print "--------------------------------------------------------------"
pdl_masked=np.array(masked)
pdl_chi_now=((pdl_masked*pdl_res_JOINT)**2)/((np.array(e_flux_unc))**2)
pdl_chi_now[np.isnan(pdl_chi_now)]=0
chi_joint=np.sum(pdl_chi_now)
chi_joint=(chi_joint/(n_unc-n_mod_fixed-nf-1))**0.5
rms=np.std(pdl_masked*pdl_res_JOINT)
j1=int(0.4*n_unc)
j2=int(0.6*n_unc)
rms=np.std(pdl_res_JOINT[j1:j2])
pdl_flux_unc_now=np.array(flux_unc)
med_flux=np.median(pdl_flux_unc_now[j1:j2])
title="X="+str(chi_joint)+" T="+str(age_min)+" ("+str(age_min_mass)+") Z="+str(met_min)+" ("+str(met_min_mass)+") Av="+str(Av_min)+" z="+str(redshift)+" sigma="+str(sigma)
ssp.plot_results(plot,wave_unc,pdl_output,out_ps,title)
MASS=ML*med_flux
lML=np.log10(ML)
print "MSP CHISQ="+str(chi_joint)+" AGE="+str(age_min)+"+-"+str(e_age_min)+" MET="+str(met_min)+"+-"+str(e_met_min)+" AV="+str(Av_min)+"+-"+str(e_Av_min)+" REDSHIFT="+str(redshift)+"+-"+str(e_redshift)+" SIGMA_DISP="+str(sigma)+"+-"+str(e_sigma)+" RMS="+str(rms)+" MED_FLUX="+str(med_flux)+" AGE_mass="+str(age_min_mass)+"+-"+str(e_age_min_mass)+" MET_mass="+str(met_min_mass)+"+-"+str(e_met_min_mass)+" MASS="+str(MASS)+" log_M/L="+str(lML)
j1=int(0.4*n_unc)
j2=int(0.6*n_unc)
wave_norm=(wave_unc[j1]+wave_unc[j2])/2.
f=open(outfile, "w")
f.write("# (1) MIN_CHISQ\n")
f.write("# (2) LW Age (Gyr)\n")
f.write("# (3) LW Age error\n")
f.write("# (4) LW metallicity\n")
f.write("# (5) LW metallicity error\n")
f.write("# (6) Av\n")
f.write("# (7) AV error\n")
f.write("# (8) redshift \n")
f.write("# (9) redshift error\n")
f.write("# (10) velocity dispersion sigma, in AA\n")
f.write("# (11) velocity dispersion error\n")
f.write("# (12) median_FLUX\n")
f.write("# (13) redshift_ssp\n")
f.write("# (14) med_flux \n")
f.write("# (15) StdDev_residual \n")
f.write("# (16) MW Age (Gyr)\n")
f.write("# (17) MW Age error\n")
f.write("# (18) MW metallicity\n")
f.write("# (19) MW metallicity error\n")
f.write("# (20) Systemic Velocity km/s \n")
f.write("# (21) Log10 Average Mass-to-Light Ratio \n")
f.write("# SSP_SFH $back_list \n")
f.write("# SSP_KIN $back_list2 \n")
f.write("# WAVE_NORM $wave_norm AA\n")
if chi_joint == 0:
    chi_joint=1
f.write(str(chi_joint)+","+str(age_min)+","+str(e_age_min)+","+str(met_min)+","+str(e_met_min)+","+str(Av_min)+","+str(e_Av_min)+","+str(redshift)+","+str(e_redshift)+","+str(sigma)+","+str(e_sigma)+","+str(FLUX)+","+str(redshift_abs)+","+str(med_flux)+","+str(rms)+","+str(age_min_mass)+","+str(e_age_min_mass)+","+str(met_min_mass)+","+str(e_met_min_mass)+","+str(SYS_VEL)+","+str(lML)+"\n")
f.close
sec_end=ssp.print_time()
sec_total=sec_end-sec_ini
print "# SECONDS = "+str(sec_total)
#
# Write the output!
#
#

max_stars_repo_name: Delosari/dazer

max_issues_repo_name: Delosari/dazer

max_forks_repo_name: Delosari/dazer

reason: import numpy,from numpy,from scipy


text:
import socket
import struct
import traceback
import logging
import time
import numpy as np
import queue
import random,threading,time
from translate import pose_predict

s = socket.socket(socket.AF_INET, socket.SOCK_STREAM)
port = 60000
s.connect(("192.168.50.14", port))
print('socket created ')

print('socket listensing ... ')

# def kinect_transform_to_model(poses):
#     new_pose = np.zeros((poses.shape[0],poses.shape[1]+3))

#     # base
#     new_pose[:,0:3] = poses[:,0:3]
#     # left leg
#     new_pose[:,3:12] = poses[:,54:63]
#     new_pose[:,12:15] = poses[:,63:66]
#     new_pose[:,15:18] = poses[:,63:66]

#     # right leg
#     new_pose[:,18:27] = poses[:,66:75]
#     new_pose[:,27:30] = poses[:,75:78]
#     new_pose[:,30:33] = poses[:,75:78]

#     #spine
#     new_pose[:,33:42] = poses[:,3:12]

#     # head
#     new_pose[:,42:48] = poses[:,78:84]

#     # left arm
#     new_pose[:,48:69] = poses[:,12:33]
#     new_pose[:,69:72] = poses[:,27:30]

#     # right arm
#     new_pose[:,72:93] = poses[:,33:54]
#     new_pose[:,93:96] = poses[:,48:51]

#     return new_pose

# def model_transform_to_kinect(poses):
#     new_pose = np.zeros((poses.shape[0],poses.shape[1]-3))

#     #base
#     new_pose[:,0:3] = poses[:,0:3]

#     #spine
#     new_pose[:,3:12] = poses[:,33:42]

#     # left arm
#     new_pose[:,12:33] = poses[:,48:69]
#     # new_pose[:,27:30] = poses[:,69:72]

#     # right arm
#     new_pose[:,33:54] = poses[:,72:93]
#     # new_pose[:,48:51] = poses[:,93:96]

#     # left leg
#     new_pose[:,54:66] = poses[:,3:15]
#     # new_pose[:,63:66] = poses[:,12:15]

#     # right leg
#     new_pose[:,66:75] = poses[:,18:27]
#     # new_pose[:,75:78] = poses[:,27:30]

#     # head
#     new_pose[:,78:84] = poses[:,42:48]

def reciveing():

    try:
        bytes_received = s.recv(388)
        pose_received = bytes_received[0:384]
        array_received = np.frombuffer(pose_received, dtype=np.float32) #converting into float array
        return array_received
        # print(array_received.shape)
        # input_pose = np.vstack((input_pose,array_received))
        # input_pose = transform(input_pose)

    except Exception as e:
        logging.error(traceback.format_exc())
        print("error")

def sending(result):
    print("sending result")
    try:
        bytes_to_send = struct.pack('%sf' % len(result), *result) #converting float to byte
        s.sendall(bytes_to_send) #sending back
    except Exception as e:
        print("sending result error!")

class Producer(threading.Thread):
    def __init__(self,t_name,data):
        threading.Thread.__init__(self,name=t_name)
        self.data =data

    def run(self):
        while 1 :
            input_pose = reciveing()
            if(len(self.data)==102):
                self.data[0:101] = self.data[1:102]
                self.data[101] = input_pose
            else:
                self.data.append(input_pose)

"""
q = []

model = pose_predict()
producer = Producer("get.",q)
producer.start()

while 1 :
    if(len(q)==102):
        poses = q
        poses = np.array(poses)
        poses = kinect_transform_to_model(poses)

        result = model.sample(poses)
        result = model_transform_to_kinect(result)
        result = model(input_pose)

        sending(result)

    else:
        continue

"""
import csv
import pandas as pd

save = open("test.csv", mode="w", newline="")
writer = csv.writer(save)

input_pose = np.zeros((100,99))

data = pd.read_csv('IRL_0_pose_4.csv', encoding = 'big5',header = None, dtype = np.float32)
data = data.to_numpy()
data = data[:,0:-1]

input_pose[:,0:96] = data[0:100,:]

header = ["ground"]
writer.writerow(header)
for idx,i in enumerate(data[100:110,:]):
    if(idx%2 ==0):
        writer.writerow(i)

model = pose_predict()

result = model.sample(input_pose)
header = ["predict"]
writer.writerow(header)
for i in result:
    writer.writerow(i)

ext: py

reason: import numpy


text:
import cv2
import numpy as np

fileName = ['9','8','7','6','5','4','3','2','1','0']

img = cv2.imread('/home/spline/document-scanner/box.jpg')

gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
gray = cv2.bilateralFilter(gray, 11, 17, 17)

kernel = np.ones((5,5),np.uint8)
erosion = cv2.erode(gray,kernel,iterations = 2)
kernel = np.ones((4,4),np.uint8)
dilation = cv2.dilate(erosion,kernel,iterations = 2)

edged = cv2.Canny(dilation, 30, 200)

_, contours, hierarchy = cv2.findContours(edged, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)

rects = [cv2.boundingRect(cnt) for cnt in contours]
rects = sorted(rects,key=lambda  x:x[1],reverse=True)

i = -1
j = 1
y_old = 5000
x_old = 5000
for rect in rects:
    x,y,w,h = rect
    area = w * h

    if area > 47000 and area < 70000:

        if (y_old - y) > 200:
            i += 1
            y_old = y

        if abs(x_old - x) > 300:
            x_old = x
            x,y,w,h = rect

            out = img[y+10:y+h-10,x+10:x+w-10]
            cv2.imwrite('cropped\\' + fileName[i] + '_' + str(j) + '.jpg', out)

            j+=1

max_stars_repo_name: rbsuchetas/object-tracker

max_issues_repo_name: rbsuchetas/object-tracker

max_forks_repo_name: rbsuchetas/object-tracker


text:
'''
Ensemble script for ImageNet
Copyright (c) Junming CHEN, 2018
'''
from __future__ import print_function

import argparse
import os
import shutil
import time
import random

import torch
import torch.nn as nn
import torch.nn.parallel
import torch.backends.cudnn as cudnn
import torch.optim as optim
import torch.utils.data as data
import torchvision.transforms as transforms
import torchvision.datasets as datasets
import torchvision.models as models
import models.imagenet as customized_models
import hashlib
import numpy as np
import pdb
from itertools import combinations
import matplotlib.pyplot as plt
import heapq
import json

from utils import Bar, Logger, AverageMeter, accuracy, mkdir_p, savefig

############################################################
MAX_ENS_NUM = 8
TEMs = [0.1, 0.3, 1.0, 3.0, 10.0, 100.0]
BestT = 100.0
Path = 'Ensemble/imagenet'
############################################################

def main():

    Models = ['vgg19_bn', 'resnet152', 'densenet161', 'densenet121', 'densenet201', 'resnet101', 'densenet169', 'resnet50']

    m = nn.Softmax()

    Targets = torch.load(Path +'/{}/Targets of {}.pt'.format(Models[0], Models[0]), map_location=lambda storage, loc: storage)

    OutputsList = []
    for model in Models:
        Logits = torch.load(Path +'/{}/Logit Outputs of {}.pt'.format(model, model), map_location=lambda storage, loc: storage)
        OutputsList.append(Logits)
    print('\nOutputlist')
    LogitTestError = test(OutputsList, Targets)
    f = open(Path +'/LogitTestError.json', 'w')
    JJ = json.dumps(LogitTestError)
    f.write(JJ)
    f.close()

        # var = torch.autograd.Variable(Logits, volatile=True)
        # SoftmaxOutList.append(m(var))

    # Temperature
    SoftmaxOutList_Tem_Dict = {}
    for T in TEMs:
        SoftmaxOutList_Tem_Dict[T] = []
    # for model in Models:
    #     TempOutsDict = torch.load(Path +'/{}/Softmax_Outs_tem_Dict of {}.pt'.format(model, model), map_location=lambda storage, loc: storage)
    #     for T in TEMs:
    #         SoftmaxOutList_Tem_Dict[T].append(TempOutsDict[T])

    SoftmaxTestError_Tem_Dict = {}
    for T in TEMs:
        print('\nSoftmaxOutList_Tem_Dict: temperature: {}'.format(T))
        SoftmaxTestError_Tem_Dict[T] = test(SoftmaxOutList_Tem_Dict[T], Targets)
    f = open(Path +'/SoftmaxTestError_Tem_Dict.json', 'w')
    JJ = json.dumps(SoftmaxTestError_Tem_Dict)
    f.write(JJ)
    f.close()

    print('\nSoftmaxTestError_GM_BestT')
    SoftmaxTestError_GM_BestT = test(SoftmaxOutList_Tem_Dict[BestT], Targets, 'GM')
    f = open(Path +'/SoftmaxTestError_GM_BestT.json', 'w')
    JJ = json.dumps(SoftmaxTestError_GM_BestT)
    f.write(JJ)
    f.close()

    SoftmaxOutList = []
    for model in Models:
        SoftmaxOuts = torch.load(Path +'/{}/Softmax Outputs of {}.pt'.format(model, model), map_location=lambda storage, loc: storage)
        SoftmaxOutList.append(SoftmaxOuts)

    print('\nSoftmaxOutList')
    SoftmaxTestError = test(SoftmaxOutList, Targets)
    f = open(Path +'/SoftmaxTestError.json', 'w')
    JJ = json.dumps(SoftmaxTestError)
    f.write(JJ)
    f.close()

    print('\nSoftmaxTestError_RMS')
    SoftmaxTestError_RMS = test(SoftmaxOutList, Targets, 'RMS')
    f = open(Path +'/SoftmaxTestError_RMS.json', 'w')
    JJ = json.dumps(SoftmaxTestError_RMS)
    f.write(JJ)
    f.close()

    print('\nSoftmaxTestError_GM')
    SoftmaxTestError_GM = test(SoftmaxOutList, Targets, 'GM')
    f = open(Path +'/SoftmaxTestError_GM.json', 'w')
    JJ = json.dumps(SoftmaxTestError_GM)
    f.write(JJ)
    f.close()

    print('\nSoftmaxTestError_HM')
    SoftmaxTestError_HM = test(SoftmaxOutList, Targets, 'HM')
    f = open(Path +'/SoftmaxTestError_HM.json', 'w')
    JJ = json.dumps(SoftmaxTestError_HM)
    f.write(JJ)
    f.close()

        # for T in TEMs:
        #     SoftmaxOutList_Tem_Dict[T].append(m(var/T))

    logger.close()

    print('Done')

def RMS(combination):
    ms = 0
    for array in combination:
        ms = ms + array**2
    ms = ms / len(combination)
    rms = np.sqrt(ms)
    return rms

def GM(combination):
    gm = 1.0
    for array in combination:
        gm = gm * array
    gm = gm ** (1.0/len(combination))
    return gm

def HM(combination):
    hm = 0
    for array in combination:
        hm = hm + (1.0/array)
    hm = float(len(combination)) / hm
    return hm

def EnsCombin(OutputsList, n, method = 'AM'):
    EnsCombination = []
    for combination in list(combinations(OutputsList, n)):
        if method == 'AM':
            Sum = 0
            for array in combination:
                Sum += array
            ensemble = Sum / float(n)
        elif method == 'RMS':
            ensemble = RMS(combination)
        elif method == 'GM':
            ensemble = GM(combination)
        elif method == 'HM':
            ensemble = HM(combination)
        else:
            print('Wrong: No such method')
        EnsCombination.append(ensemble)
    return EnsCombination

def test(OutputsList, Targets, method = 'AM'):
    tt = Targets.cpu().numpy().reshape(-1,)

    TestError_top1 = []
    TestError_top5 = []

    for i in range(MAX_ENS_NUM):
        print(i)
        m = 1 # Softmax Only
        TElist_top1 = []
        Aver_top1 = 100
        Best_top1 = 100
        Worst_top1 = 0
        # WrongQuantity_top1 = len(Targets)
        TElist_top5 = []
        Aver_top5 = 100
        Best_top5 = 100
        Worst_top5 = 0
        # WrongQuantity_top5 = len(Targets)
        for ensemble in EnsCombin(OutputsList, i+1, method):
            WrongQuantity_top1 = len(Targets)
            WrongQuantity_top5 = len(Targets)
            oo = ensemble.cpu().numpy()
            pred_top1 = []
            for row, target in zip(oo, tt):
                index_top1 = np.argwhere(row == max(row)).reshape(-1,)
                if len(index_top1) != 1:
                    for idx in index_top1:
                        if idx == target:
                            print(len(index_top1),' ',row[idx],' ',method, file = fff)
                        #     index_top1 = np.array([target])
                        #     break
                        # else:
                        #     index_top1 = np.array([idx])
                    # pdb.set_trace()
                    index_top1 = np.array([index_top1[0]])
                # pred_top1.append(index_top1)

                index_top5 = index_top1
                if index_top1[0] == target:
                    WrongQuantity_top1 -= 1;
                    WrongQuantity_top5 -= 1;
                else:
                    for value in heapq.nlargest(5, row):
                        idx = np.argwhere(row == value).reshape(-1,)
                        if len(idx) != 1:
                            for ID in idx:
                                if ID == target:
                                    print(len(idx),' ',row[ID],' ',method, file = ff)
                            idx = idx[0]
                        if idx == target:
                            # index_top5 = idx
                            WrongQuantity_top5 -= 1;
                            break
                # pred_top5.append(index_top5)

            # pred_top1 = np.asarray(pred_top1).reshape(-1,)
            # unequal = (pred_top1 != tt)
            # unique, counts = np.unique(unequal, return_counts=True)
            # WrongQuantity = dict(zip(unique, counts))[True]
            # TestErr = WrongQuantity / float(len(unequal))
            TestErr = WrongQuantity_top1 / float(len(Targets))
            TElist_top1.append(TestErr)
            if TestErr < Best_top1:
                Best_top1 = TestErr
            if TestErr > Worst_top1:
                Worst_top1 = TestErr

            # pred_top5 = np.asarray(pred_top5).reshape(-1,)
            # unequal = (pred_top5 != tt)
            # unique, counts = np.unique(unequal, return_counts=True)
            # WrongQuantity = dict(zip(unique, counts))[True]
            # TestErr = WrongQuantity / float(len(unequal))
            TestErr = WrongQuantity_top5 / float(len(Targets))
            TElist_top5.append(TestErr)
            if TestErr < Best_top5:
                Best_top5 = TestErr
            if TestErr > Worst_top5:
                Worst_top5 = TestErr

        Aver_top1 = np.mean(TElist_top1)
        Aver_top5 = np.mean(TElist_top5)

        TestError_top1.append([Best_top1, Aver_top1, Worst_top1])
        TestError_top5.append([Best_top5, Aver_top5, Worst_top5])

        TestError = {'top1': TestError_top1, 'top5': TestError_top5}
    return TestError

# def save_checkpoint(state, is_best, checkpoint='checkpoint', filename='checkpoint.pth.tar'):
#     filepath = os.path.join(checkpoint, filename)
#     torch.save(state, filepath)
#     if is_best:
#         shutil.copyfile(filepath, os.path.join(checkpoint, 'model_best.pth.tar'))

if __name__ == '__main__':
    main()

max_stars_repo_name: JeremyCJM/Ensemble

max_issues_repo_name: JeremyCJM/Ensemble

max_forks_repo_name: JeremyCJM/Ensemble


text:
import json

import numpy as np

def get_data_folder_name():
    return "data/"

def get_json_file_name(endpoint, num_pages):
    return get_data_folder_name() + f"{endpoint}_{num_pages}.json"

def get_snapshot_fname():
    return get_data_folder_name() + "IStoreService.json"

def load_json(fname):
    with open(fname, "r", encoding="utf8") as f:
        data = json.load(f)

    return data

def save_json(data, fname):
    with open(fname, "w", encoding="utf8") as f:
        json.dump(data, f)

    return

def save_npy(data, fname):
    # NB: data type is uint32 (max value ~ 4 billion), which is far greater than max appID (~ 2 million)

    np.save(
        fname,
        np.asarray(data, dtype=np.uint32),
        allow_pickle=False,
        fix_imports=False,
    )

    return

def save_app_ids(app_ids, fname_root="app_ids"):
    save_json(data=app_ids, fname=get_data_folder_name() + f"{fname_root}.json")
    save_npy(data=app_ids, fname=get_data_folder_name() + f"{fname_root}.npy")

    return

max_stars_repo_name: woctezuma/steam-popular-appids

max_issues_repo_name: woctezuma/steam-popular-appids

max_forks_repo_name: woctezuma/steam-popular-appids


text:
import ast
import networkx as nx
import matplotlib.pyplot as plt
from typing import Any, Optional, NoReturn, Tuple
from collections import defaultdict

class ASTVisualizer(ast.NodeVisitor):
    def __init__(self):
        self.stack = []
        self.ast = nx.Graph()
        self.label_dict = {}
        self.class_count = defaultdict(int)

    def extend_ast(self, node: ast.AST, node_label: str):
        class_name = node.__class__.__name__
        node_name = f"{class_name}_{self.class_count[class_name]}"
        self.class_count[class_name] += 1
        parent_name = None
        if self.stack:
            parent_name = self.stack[-1]

        self.stack.append(node_name)
        self.ast.add_node(node_name)
        self.label_dict[node_name] = node_label

        if parent_name:
            self.ast.add_edge(parent_name, node_name)

        super(self.__class__, self).generic_visit(node)
        self.stack.pop()

    def visit_FunctionDef(self, node: ast.FunctionDef) -> Any:
        label = f"{node.__class__.__name__}\nname: {node.name}"
        self.extend_ast(node, label)

    def visit_arg(self, node: ast.arg) -> Any:
        label = f"{node.__class__.__name__}\nname: {node.arg}"
        self.extend_ast(node, label)

    def visit_Name(self, node: ast.Name) -> Any:
        label = f"{node.__class__.__name__}\nid: {node.id}"
        self.extend_ast(node, label)

    def visit_Constant(self, node: ast.Constant) -> Any:
        label = f"{node.__class__.__name__}\nvalue: {node.value}"
        self.extend_ast(node, label)

    def generic_visit(self, node: ast.AST) -> Any:
        label = node.__class__.__name__
        self.extend_ast(node, label)

    def visualize(self, save_fn: Optional[str] = None, figsize: Optional[Tuple[int, int]] = None) -> NoReturn:
        plt.figure(1, figsize=figsize)
        nx.draw(self.ast, labels=self.label_dict, with_labels=True)
        if save_fn:
            plt.savefig(save_fn)
        plt.show()

def visualize_ast(input_fn: str, output_fn: Optional[str] = None, figsize: Optional[Tuple[int, int]] = None):
    with open(input_fn, 'r') as fin:
        code = fin.read()

    root = ast.parse(code)
    print(ast.dump(root,  indent=4))
    ast_vis = ASTVisualizer()
    ast_vis.visit(root)
    ast_vis.visualize(save_fn=output_fn, figsize=figsize)

if __name__ == "__main__":
    visualize_ast("fibonacci.py", "artifacts/fibonacci_ast.png", figsize=(15, 15))

reason: import networkx


text:
import h5py
import numpy as np

def print_group(f, group_name, indent=''):
    view = f[group_name]
    for k in view.keys():
        if isinstance(view[k], h5py.Group):
            print_group(view, k, indent='  ')
        else:
            print(f"{indent}{group_name}/{k}: {type(view[k])}")

with h5py.File('data_dump.h5', 'r') as f:
    print(f.keys())
    for iteration_group in f.keys():
        print_group(f, iteration_group)

max_stars_repo_name: domischi/SpringBox

max_issues_repo_name: domischi/SpringBox

max_forks_repo_name: domischi/SpringBox


text:
import networkx as nx
import pytest

from qupid.matching import hopcroft_karp_matching

@pytest.fixture
def simple_graph():
    G = nx.Graph()

    G.add_node("1")
    G.add_node("2")
    G.add_node("3")

    G.add_node("A")
    G.add_node("B")
    G.add_node("C")

    G.add_edge("1", "A")
    G.add_edge("1", "C")

    G.add_edge("2", "B")
    G.add_edge("2", "C")

    G.add_edge("3", "A")
    G.add_edge("3", "B")

    return G

def test_hk(simple_graph):
    all_matches = set()
    for i in range(20):
        M = hopcroft_karp_matching(simple_graph)
        match = frozenset(tuple((x, M[x]) for x in M))
        all_matches.add(match)
    print(all_matches)

    exp_matches = {
        frozenset((("1", "A"), ("2", "C"), ("3", "B"))),
        frozenset((("1", "C"), ("2", "B"), ("3", "A"))),
    }
    assert all_matches == exp_matches

max_stars_repo_name: gibsramen/matchlock

max_issues_repo_name: gibsramen/matchlock

max_forks_repo_name: gibsramen/matchlock


text:
import numpy as np
from easydict import EasyDict as edict

config = edict()

config.num_label = 10
config.img_width = 280
config.img_height = 32
config.seq_length = config.img_width//8
config.image_path = 'images'
config.num_classes = 5990
config.num_hidden = 100
config.num_lstm_layer = 2
config.use_lstm = False
config.no4x1pooling=False
config.to_gray = True

#config.PIXEL_MEANS = np.array([103.939, 116.779, 123.68])

# network related params
#config.FIXED_PARAMS = ['conv1', 'conv2']

default = edict()
default.network = 'simplenet' # simplenet, resnet, mobilenet, densenet
default.dataset = 'chinese' #chinese or vgg
#default.network = 'resnet'
default.pretrained = ''
default.pretrained_epoch = 0
default.lr = 0.0001
default.lr_step = '5,10'
# default dataset
default.dataset_path = ''
# default training
default.frequent = 20
default.batch_size = 128
default.kvstore = 'device'
# default e2e
default.prefix = 'model/e2e'
default.epoch = 20

network = edict()

dataset = edict()

dataset.chinese = edict()
dataset.chinese.dataset_path = '/gpu/data2/jiaguo/ocr_cn_dataset'
dataset.chinese.prefix = 'model/chinese'

dataset.vgg = edict()
dataset.vgg.dataset_path = '/gpu/data2/jiaguo/ocr_dataset'
dataset.vgg.prefix = 'model/vgg'
dataset.vgg.num_label = 23
dataset.vgg.img_width = 96
dataset.vgg.img_height = 32
dataset.vgg.seq_length = dataset.vgg.img_width//8
dataset.vgg.num_classes = 62
dataset.vgg.no4x1pooling=False
dataset.vgg.to_gray = False

def generate_config(_network, _dataset):
    if _network in network:
      for k, v in network[_network].items():
          if k in config:
              config[k] = v
          elif k in default:
              default[k] = v
    if _dataset in dataset:
      for k, v in dataset[_dataset].items():
          if k in config:
              config[k] = v
          elif k in default:
              default[k] = v

max_stars_repo_name: deepinsight/insightocr

max_issues_repo_name: gregbugaj/insightocr

max_forks_repo_name: gregbugaj/insightocr


text:
import pathlib

import numpy as np
from scipy import sparse
import pandas as pd
from sklearn.preprocessing import normalize
from sklearn.utils.extmath import safe_sparse_dot
from nilearn import image

from neuroquery.img_utils import get_masker
from neuroquery import tokenization, smoothed_regression, ridge

_MAX_SIMILAR_DOCS_RETURNED = 100

class NeuroQueryModel:
    """Text -> brain map encoder.

    It encodes text into statistical maps of the brain and also provides a list
    of related terms.

    It can be initialized with a fitted regression model
    (`neuroquery.smoothed_regression.SmoothedRegression`) or loaded using
    `from_data_dir`. Most users will probably load a pre-trained model with
    `from_data_dir`.

    Parameters
    ----------
    vectorizer : `neuroquery.tokenization.TextVectorizer`
        An object that transforms text into TFIDF features.

    smoothed_regression : `neuroquery.smoothed_regression.SmoothedRegression`
        A reduced-rank regression that combines feature smoothing, projection,
        and linear regression. The input features must correspond to the
        outputs of `vectorizer`.

    mask_img : Nifti1Image
        Mask of the regression targets. The non-zero voxels correspond to the
        dependent variables.

    corpus_info : dict, optional (default=None)
        Data required to report which studies are most relevant for a query.
        Must contain:
            - "metadata": pandas DataFrame, each row describing a study
            - "tfidf": scipy sparse matrix or numpy array, TFIDF features for
              the documents. Rows must correspond to the same studies as in
              "metadata", and columns to the terms in the vectorizer's
              vocabulary.
        If corpus_info is not available the model will not report most similar
        studies.
    """

    @classmethod
    def from_data_dir(cls, model_dir):
        """Load a pre-trained TextToBrain model.

        Parameters
        ----------
        model_dir : str
            path to a directory containing the serialized trained model.
            The directory must be organized as the one returned by
            `neuroquery.datasets.fetch_neuroquery_model`, except that
            `corpus_metadata.csv` and `corpus_tfidf.npz` are optional.
        """
        model_dir = pathlib.Path(model_dir)
        vectorizer = tokenization.TextVectorizer.from_vocabulary_file(
            str(model_dir / "vocabulary.csv"),
            voc_mapping="auto",
            add_unigrams=False,
        )
        regression = smoothed_regression.SmoothedRegression.from_data_dir(
            str(model_dir)
        )
        mask_img = image.load_img(str(model_dir / "mask_img.nii.gz"))
        corpus_tfidf = model_dir / "corpus_tfidf.npz"
        corpus_metadata = model_dir / "corpus_metadata.csv"
        if corpus_tfidf.is_file() and corpus_metadata.is_file():
            corpus_info = {}
            corpus_info["tfidf"] = sparse.load_npz(str(corpus_tfidf))
            corpus_info["metadata"] = pd.read_csv(
                str(corpus_metadata), encoding="utf-8"
            )
        else:
            corpus_info = None
        return cls(vectorizer, regression, mask_img, corpus_info=corpus_info)

    def to_data_dir(self, model_dir):
        """Save the model so it can later be loaded with `from_data_dir`."""
        model_dir = pathlib.Path(model_dir)
        model_dir.mkdir(parents=True, exist_ok=True)
        self.vectorizer.to_vocabulary_file(str(model_dir / "vocabulary.csv"))
        self.smoothed_regression.to_data_dir(model_dir)
        self.get_masker().mask_img_.to_filename(
            str(model_dir / "mask_img.nii.gz")
        )
        if self.corpus_info is not None:
            sparse.save_npz(
                str(model_dir / "corpus_tfidf.npz"),
                sparse.csr_matrix(self.corpus_info["tfidf"]),
            )
            self.corpus_info["metadata"].to_csv(
                str(model_dir / "corpus_metadata.csv"), index=False
            )

    def __init__(
        self, vectorizer, smoothed_regression, mask_img, corpus_info=None
    ):
        self.vectorizer = vectorizer
        self.smoothed_regression = smoothed_regression
        self.mask_img = mask_img
        self.corpus_info = corpus_info

    def full_vocabulary(self):
        """All the terms recognized by the model."""
        return self.vectorizer.get_vocabulary()

    def _supervised_features(self):
        if not hasattr(
            self.smoothed_regression.regression_, "selected_features_"
        ):
            return np.arange(
                self.smoothed_regression.regression_.coef_.shape[1]
            )
        return self.smoothed_regression.regression_.selected_features_

    def supervised_vocabulary(self):
        """Terms selected as features for the supervised regression."""
        return np.asarray(self.full_vocabulary())[self._supervised_features()]

    def document_frequencies(self):
        if self.corpus_info is None:
            return None
        if not hasattr(self, "document_frequencies_"):
            document_frequencies = (self.corpus_info["tfidf"] > 0).sum(axis=0)
            document_frequencies = np.asarray(document_frequencies).ravel()
            self.document_frequencies_ = pd.Series(
                document_frequencies, index=self.full_vocabulary()
            )
        return self.document_frequencies_

    def _similar_words(self, tfidf, vocabulary=None):
        if vocabulary is None:
            vocabulary = self.full_vocabulary()
        if sparse.issparse(tfidf):
            tfidf = tfidf.A.squeeze()
        similar = pd.Series(tfidf, index=vocabulary).sort_values(
            ascending=False
        )
        return similar[similar > 0]

    def similar_documents(self, tfidf):
        if self.corpus_info is None:
            return None
        similarities = safe_sparse_dot(
            tfidf, self.corpus_info["tfidf"].T, dense_output=True
        ).ravel()
        order = np.argsort(similarities)[::-1]
        order = order[similarities[order] > 0][:_MAX_SIMILAR_DOCS_RETURNED]
        ordered_simil = similarities[order]
        similar_docs = (
            self.corpus_info["metadata"].iloc[order].reset_index(drop=True)
        )
        similar_docs["similarity"] = ordered_simil
        return similar_docs

    def _beta_norms(self):
        return np.linalg.norm(
            self.smoothed_regression.regression_.coef_, axis=0
        )

    def get_masker(self):
        if not hasattr(self, "masker_"):
            self.masker_ = get_masker(self.mask_img)
        return self.masker_

    def _supervised_vocabulary_set(self):
        if not hasattr(self, "supervised_vocabulary_set_"):
            self.supervised_vocabulary_set_ = set(self.supervised_vocabulary())
        return self.supervised_vocabulary_set_

    def transform(self, documents):
        """Transform a set of documents

        Parameters
        ----------
        documents : list or array of str
            the text snippets to transform

        Returns
        -------
        list of dict, each containing:
            - "brain_map": a nifti image of the most relevant brain regions.
            - "raw_tfidf": the vectorized documents.
            - "smoothed_tfidf": the tfidf after semantic smoothing.
            - "z_map" is an alias for "brain_map" for backwards compatibility

        """
        raw_tfidf = self.vectorizer.transform(documents)
        raw_tfidf = normalize(raw_tfidf, copy=False)
        self.smoothed_regression.regression_.intercept_ = 0.0
        brain_maps = self.smoothed_regression.transform_to_brain_maps(
            raw_tfidf
        )
        masker = self.get_masker()
        brain_maps_unmasked = list(map(masker.inverse_transform, brain_maps))
        smoothed_tfidf = self.smoothed_regression.smoothing_.transform(
            raw_tfidf
        )
        smoothed_tfidf = normalize(smoothed_tfidf, copy=False)
        return {
            "brain_map": brain_maps_unmasked,
            "z_map": brain_maps_unmasked,
            "raw_tfidf": raw_tfidf,
            "smoothed_tfidf": smoothed_tfidf,
        }

    def __call__(self, document):
        """Transform a document

        Parameters
        ----------
        document : str
            the text to transform

        Returns
        -------
        dict containing:
            - "brain_map": a nifti image of the most relevant brain regions.
            - "similar_words": pandas DataFrame containing related terms.
                - "similarity" is how much the term is related.
                - "weight_in_brain_map" is the contribution of the term in the
                  predicted "brain_map".
                - "weight_in_query" is the TFIDF of the term in `document`.
            - "similar_documents": if no corpus_info was provided, this is
              `None`. Otherwise it is a DataFrame containing information about
              the most relevant studies.
            - "highlighted_text": an XML document showing which terms were
              recognized in the provided text.
            - "smoothed_tfidf": the tfidf after semantic smoothing.
            - "raw_tfidf": the vectorized documents.
            - "z_map" is an alias for "brain_map" for backwards compatibility
        """
        self.vectorizer.tokenizer.keep_pos = True
        result = self.transform([document])
        result = {k: v[0] for k, v in result.items()}
        similar_words = pd.DataFrame(
            {
                "similarity": self._similar_words(result["smoothed_tfidf"]),
                "weight_in_query": self._similar_words(result["raw_tfidf"]),
                "weight_in_brain_map": self._similar_words(
                    result["smoothed_tfidf"][self._supervised_features()]
                    * self._beta_norms(),
                    self.supervised_vocabulary(),
                ),
            },
            columns=["similarity", "weight_in_brain_map", "weight_in_query"],
        )
        similar_words.fillna(0.0, inplace=True)
        similar_words.sort_values(
            by="weight_in_brain_map", ascending=False, inplace=True
        )
        doc_freq = self.document_frequencies()
        if doc_freq is not None:
            similar_words["n_documents"] = doc_freq.loc[similar_words.index]
            similar_words = similar_words.loc[
                :,
                [
                    "similarity",
                    "weight_in_brain_map",
                    "weight_in_query",
                    "n_documents",
                ],
            ]
        result["similar_words"] = similar_words
        result["similar_documents"] = self.similar_documents(
            result["smoothed_tfidf"]
        )
        self._supervised_vocabulary_set()
        result[
            "highlighted_text"
        ] = self.vectorizer.tokenizer.highlighted_text(
            lambda w: {
                "in_model": (
                    "true" if w in self.supervised_vocabulary_set_ else "false"
                )
            }
        )
        return result

class SimpleEncoder:
    """Basic text to brain map encoder"""

    @classmethod
    def from_data_dir(cls, model_dir):
        model_dir = pathlib.Path(model_dir)
        vectorizer = tokenization.TextVectorizer.from_vocabulary_file(
            str(model_dir / "vocabulary.csv"),
            voc_mapping="auto",
            add_unigrams=False,
        )
        regression = ridge.FittedLinearModel.from_data_dir(model_dir)
        mask_img = image.load_img(str(model_dir / "mask_img.nii.gz"))
        return cls(vectorizer, regression, mask_img)

    def to_data_dir(self, model_dir):
        """Save the model so it can later be loaded with `from_data_dir`."""
        model_dir = pathlib.Path(model_dir)
        model_dir.mkdir(parents=True, exist_ok=True)
        self.vectorizer.to_vocabulary_file(str(model_dir / "vocabulary.csv"))
        self.regression.to_data_dir(model_dir)
        self.get_masker().mask_img_.to_filename(
            str(model_dir / "mask_img.nii.gz")
        )

    def __init__(self, vectorizer, regression, mask_img):
        self.vectorizer = vectorizer
        self.regression = regression
        self.mask_img = mask_img

    def get_masker(self):
        if not hasattr(self, "masker_"):
            self.masker_ = get_masker(self.mask_img)
        return self.masker_

    def __call__(self, document):
        self.vectorizer.tokenizer.keep_pos = True
        self.regression.intercept_ = 0.0
        result = {}
        tfidf = self.vectorizer.transform([document])
        masked_map = self.regression.predict(tfidf).squeeze()
        result["brain_map"] = self.get_masker().inverse_transform(masked_map)
        result[
            "highlighted_text"
        ] = self.vectorizer.tokenizer.highlighted_text()
        return result

    def full_vocabulary(self):
        """All the terms recognized by the model."""
        return self.vectorizer.get_vocabulary()

max_stars_repo_name: RaphaelMeudec/neuroquery

max_issues_repo_name: RaphaelMeudec/neuroquery

max_forks_repo_name: RaphaelMeudec/neuroquery

reason: import numpy,from scipy


text:
from numpy import mean
from numpy import std
from numpy import percentile

def std_dev(dataset):

    data_mean, data_std = mean(dataset), std(dataset)
    # Identify outliers
    cut_off = data_std * 3
    lower, upper = data_mean - cut_off, data_mean + cut_off

    #outliers = [x for x in data if x < lower or x > upper]
    outliers_removed = [x for x in dataset if x > lower and x < upper]

    return outliers_removed

def iqr(dataset):

    # calculate interquartile range
    q25, q75 = percentile(dataset, 25), percentile(dataset, 75)
    iqr = q75 - q25

    # calculate the outlier cutoff
    cut_off = iqr * 1.5
    lower, upper = q25 - cut_off, q75 + cut_off

    # remove outliers
    outliers_removed = [x for x in dataset if x > lower and x < upper]

    return outliers_removed

max_stars_repo_name: jseidl/pi-weather-station

max_issues_repo_name: jseidl/pi-weather-station

max_forks_repo_name: jseidl/pi-weather-station

reason: from numpy


text:
# %%
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
import mifs
# %%
dataset = pd.read_csv('data/parkinsons.csv')
X = dataset.iloc[:, 1:-1]   # input feactures
y = dataset.iloc[:, -1]     # output
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.25, random_state=0)
# %%
for i in range(1, 22):
    feat_selector = mifs.MutualInformationFeatureSelector('MRMR', k=i)
    feat_selector.fit(X_train, y_train)
    # call transform() on X to filter it down to selected features
    X_filtered = feat_selector.transform(X_train.values)
    # Create list of features
    feature_name = X_train.columns[feat_selector.ranking_]
    print(feature_name)


text:
import dataset
from dataset.Inshop import Inshop_Dataset
import torch
from collections import defaultdict
from combine_sampler import CombineSampler, CombineSamplerAdvanced, CombineSamplerSuperclass, CombineSamplerSuperclass2
import numpy as np
import pickle
import matplotlib.pyplot as plt
from  matplotlib.lines import Line2D
import similarity
import sys

def create_loaders(data_root, num_classes, is_extracted, num_workers, num_classes_iter, num_elements_class, size_batch, dict_superclass=None, iterations_for_epoch=None):
    Dataset = dataset.Birds(
        root=data_root,
        labels=list(range(0, num_classes)),
        is_extracted=is_extracted,
        transform=dataset.utils.make_transform())

    Dataset_2 = dataset.Birds(
        root=data_root,
        labels=list(range(num_classes, 2*num_classes)),
        is_extracted=is_extracted,
        transform=dataset.utils.make_transform())

    ddict = defaultdict(list)
    for idx, label in enumerate(Dataset.ys):
        ddict[label].append(idx)

    list_of_indices_for_each_class = []
    for key in ddict:
        list_of_indices_for_each_class.append(ddict[key])

    ddict_2 = defaultdict(list)
    for idx, label in enumerate(Dataset_2.ys):
        ddict_2[label].append(idx)

    list_of_indices_for_each_class_2 = []
    for key in ddict_2:
        list_of_indices_for_each_class_2.append(ddict_2[key])

    dl_ev_2= torch.utils.data.DataLoader(
        Dataset_2,
        batch_size=size_batch,
        shuffle=False,
        sampler=CombineSampler(list_of_indices_for_each_class_2, num_classes_iter, num_elements_class),
        num_workers=num_workers,
        drop_last=True,
        pin_memory=True
    )

    dl_tr = torch.utils.data.DataLoader(
        Dataset,
        batch_size=size_batch,
        shuffle=False,
        sampler=CombineSampler(list_of_indices_for_each_class, num_classes_iter, num_elements_class),
        num_workers=num_workers,
        drop_last=True,
        pin_memory=True
    )

    dl_tr_superclass=None
    if(dict_superclass is not None):
        dl_tr_superclass = torch.utils.data.DataLoader(
            Dataset,
            batch_size=size_batch,
            shuffle=False,
            sampler=CombineSamplerSuperclass(list_of_indices_for_each_class, num_classes_iter, num_elements_class, dict_superclass, iterations_for_epoch),
            num_workers=num_workers,
            drop_last=True,
            pin_memory=True
        )

    if data_root == 'Stanford':
        class_end = 2 * num_classes - 2
    else:
        class_end = 2 * num_classes

    dl_ev = torch.utils.data.DataLoader(
        dataset.Birds(
            root=data_root,
            labels=list(range(num_classes, class_end)),
            is_extracted=is_extracted,
            transform=dataset.utils.make_transform(is_train=False)
        ),
            batch_size=50,
            shuffle=False,
            num_workers=1,
            pin_memory=True
    )

    dl_finetune = torch.utils.data.DataLoader(
        dataset.Birds(
            root=data_root,
            labels=list(range(num_classes)),
            is_extracted=is_extracted,
            transform=dataset.utils.make_transform(is_train=False)
        ),
        batch_size=size_batch,
        shuffle=True,
        num_workers=num_workers,
        pin_memory=True
    )

    dl_train_evaluate = torch.utils.data.DataLoader(
        dataset.Birds(
            root=data_root,
            labels=list(range(num_classes)),
            is_extracted=is_extracted,
            transform=dataset.utils.make_transform(is_train=False)
        ),
        batch_size=150,
        shuffle=False,
        num_workers=1,
        pin_memory=True
    )

    return dl_tr, dl_ev, dl_ev_2, dl_finetune, dl_train_evaluate, dl_tr_superclass

def create_loaders_Inshop(data_root, num_classes, is_extracted, num_workers, num_classes_iter, num_elements_class, size_batch, dict_class_distances=None, iterations_for_epoch=None):
    Dataset = Inshop_Dataset(
        root = data_root,
        mode = 'train',
        transform = dataset.utils.make_transform())

    ddict = defaultdict(list)
    for idx, label in enumerate(Dataset.ys):
        ddict[label].append(idx)

    list_of_indices_for_each_class = []
    for key in ddict:
        list_of_indices_for_each_class.append(ddict[key])

    dl_tr = torch.utils.data.DataLoader(
        Dataset,
        batch_size = 180,
        sampler=CombineSampler(list_of_indices_for_each_class, num_classes_iter, num_elements_class),
        shuffle = False,
        num_workers = num_workers,
        drop_last = True,
        pin_memory = True
    )

    dl_tr_ev = torch.utils.data.DataLoader(
    Dataset,
    batch_size=150,
    shuffle=False,
    num_workers=1,
    pin_memory=True
    )

    query_dataset = Inshop_Dataset(
            root = data_root,
            mode = 'query',
            transform = dataset.utils.make_transform(is_train=False
    ))

    dl_query = torch.utils.data.DataLoader(
        query_dataset,
        batch_size = 180,
        shuffle = False,
        num_workers = num_workers,
        pin_memory = True
    )

    gallery_dataset = Inshop_Dataset(
            root = data_root,
            mode = 'gallery',
            transform = dataset.utils.make_transform(
                is_train = False
    ))

    dl_gallery = torch.utils.data.DataLoader(
        gallery_dataset,
        batch_size = 180,
        shuffle = False,
        num_workers = num_workers,
        pin_memory = True
    )
    return dl_tr, dl_query, dl_gallery,dl_tr_ev

def create_loaders_finetune(data_root, num_classes, is_extracted, num_workers, size_batch):

    if data_root == 'Stanford':
        class_end = 2 * num_classes - 2
    else:
        class_end = 2 * num_classes

    dl_ev = torch.utils.data.DataLoader(
        dataset.Birds(
            root=data_root,
            labels=list(range(num_classes, class_end)),
            is_extracted=is_extracted,
            transform=dataset.utils.make_transform(is_train=False)
        ),
            batch_size=150,
            shuffle=False,
            num_workers=1,
            pin_memory=True
    )

    dl_finetune = torch.utils.data.DataLoader(
        dataset.Birds(
            root=data_root,
            labels=list(range(num_classes)),
            is_extracted=is_extracted,
            transform=dataset.utils.make_transform(is_train=False)
        ),
        batch_size=size_batch,
        shuffle=True,
        num_workers=num_workers,
        pin_memory=True
    )

    return dl_ev, dl_finetune

def create_loaders_new(data_root, num_classes, is_extracted, dict_class_distances, num_classes_iter=0, num_elements_class=0, iterations_for_epoch=200):
    Dataset = dataset.Birds(
       root=data_root,
       labels=list(range(0, num_classes)),
       is_extracted=is_extracted,
       transform=dataset.utils.make_transform())
    ddict = defaultdict(list)
    for idx, label in enumerate(Dataset.ys):
       ddict[label].append(idx)
    list_of_indices_for_each_class = []
    for key in ddict:
       list_of_indices_for_each_class.append(ddict[key])
    # with open('list_of_indices_for_each_class_stanford.pickle', 'wb') as handle:
    #     pickle.dump(list_of_indices_for_each_class, handle, protocol=pickle.HIGHEST_PROTOCOL)
    list_of_indices_for_each_class = pickle.load(open('list_of_indices_for_each_class_stanford.pickle', 'rb'))
    dl_tr = torch.utils.data.DataLoader(
       Dataset,
       batch_size=num_classes_iter * num_elements_class,
       shuffle=False,
       sampler=CombineSamplerSuperclass2(list_of_indices_for_each_class, num_classes_iter, num_elements_class, dict_class_distances, iterations_for_epoch=iterations_for_epoch),
       num_workers=8,
       drop_last=True,
       pin_memory=True
    )
    return dl_tr

def get_labeled_and_unlabeled_points(labels, num_points_per_class, num_classes=100):
    labs, L, U = [], [], []
    labs_buffer = np.zeros(num_classes)
    num_points = labels.shape[0]
    for i in range(num_points):
        if labs_buffer[labels[i]] == num_points_per_class:
            U.append(i)
        else:
            L.append(i)
            labs.append(labels[i])
            labs_buffer[labels[i]] += 1
    return labs, L, U

def get_W_gt(labels):
    n=labels.shape[0]
    W = torch.zeros(n, n)
    for i, yi in enumerate(labels):
        for j, yj in enumerate(labels[i+1:],i+1):
              W[i, j] = W[j, i] = yi==yj
    W = W.cuda()
    W.requires_grad=False
    return W

def debug_info(gtg, model):
    for name, param in gtg.named_parameters():
        if param.requires_grad:
            if param.grad is not None:
                # print(name, param.grad.data.norm(2))
                print(name, torch.mean(param.grad.data))
    for name, param in model.named_parameters():
        if param.requires_grad:
            if param.grad is not None:
                # print(name, param.grad.data.norm(2))
                print(name, torch.mean(param.grad.data))
    print("\n\n\n")

def plot_grad_flow(named_parameters,name):
    '''Plots the gradients flowing through different layers in the net during training.
    Can be used for checking for possible gradient vanishing / exploding problems.

    Usage: Plug this function in Trainer class after loss.backwards() as
    "plot_grad_flow(self.model.named_parameters())" to visualize the gradient flow'''
    ave_grads = []
    max_grads= []
    layers = []
    for n, p in named_parameters:
        if(p.requires_grad) and ("bias" not in n):
            layers.append(n)
            ave_grads.append(p.grad.abs().mean())
            max_grads.append(p.grad.abs().max())
    plt.bar(np.arange(len(max_grads)), max_grads, alpha=0.1, lw=1, color="c")
    plt.bar(np.arange(len(max_grads)), ave_grads, alpha=0.1, lw=1, color="b")
    plt.hlines(0, 0, len(ave_grads)+1, lw=2, color="k" )
    plt.xticks(range(0,len(ave_grads), 1), layers, rotation="vertical")
    plt.xlim(left=0, right=len(ave_grads))
    plt.ylim(bottom = -0.001, top=0.02) # zoom in on the lower gradient regions
    plt.xlabel("Layers")
    plt.ylabel("average gradient")
    plt.title("Gradient flow")
    plt.grid(True)
    #plt.tight_layout()
    plt.legend([Line2D([0], [0], color="c", lw=4),
                Line2D([0], [0], color="b", lw=4),
                Line2D([0], [0], color="k", lw=4)], ['max-gradient', 'mean-gradient', 'zero-gradient'])
    plt.savefig('save_graphics/'+str(name))

def plot_grad_flow_v2(named_parameters,name, e):
    plt.figure(int(e))
    ave_grads = []
    layers = []
    for n, p in named_parameters:
        if(p.requires_grad) and ("bias" not in n):
            layers.append(n)
            ave_grads.append(p.grad.abs().mean())
    plt.plot(ave_grads, alpha=0.25, color="r")
    plt.hlines(0, 0, len(ave_grads)+1, linewidth=1, color="k" )
    plt.xticks(range(0,len(ave_grads), 1), layers, rotation="vertical")
    plt.xlim(xmin=0, xmax=len(ave_grads))
    plt.xlabel("Layers")
    plt.ylabel("average gradient per layer")
    plt.title("Gradient flow")
    #plt.tight_layout()
    plt.grid(True)
    plt.savefig('save_graphics/'+str(name)+'_epoch:'+str(e))

def plot_losses(loss_class, loss_group, path_to_save='save_graphics/', file_name='model'):
    plt.figure(100)
    plt.plot(range(len(loss_class)),loss_class, c="r", alpha=1, label="classification loss")
    plt.plot(range(len(loss_group)),loss_group, c="b", alpha=1, label="group loss")
    plt.grid(True)
    plt.title('plot of losses')
    plt.xlabel("average loss per epoch")
    plt.ylabel("loss")
    plt.legend(loc='upper right')
    plt.savefig(path_to_save+file_name+'_loss_curve')
    plt.clf()
    return None

def calc_anchors(model, dataloader, net_type):
    fc7s, L = [], []
    with torch.no_grad():
        for X, Y in dataloader:
            _, fc7 = model(X.cuda())
            fc7s.append(fc7.cpu())
            L.append(Y)
    fc7, Y = torch.cat(fc7s), torch.cat(L)
    print(fc7.shape)
    print(Y.shape)
    return torch.squeeze(fc7), torch.squeeze(Y)