def showSAXSProfiles(exp_data, model_data):
    #Read experimental data
    #Read model data
    #Plot data
    from matplotlib import pyplot;
    from pylab import genfromtxt;
    mat0 = genfromtxt(exp_data);
    mat1 = genfromtxt(model_data);
    pyplot.plot(mat0[:,0], mat0[:,1], label = "Experimental");
    pyplot.plot(mat1[:,0], mat1[:,1], label = "Model");
    pyplot.legend();
    pyplot.show();

def parseResults(files):
    ''' Reads all of the results files and puts them into a list with the
    results. Returns field, dither, fiber, and redshift.

    '''

    r = []
    for f in files:
        print f
        cluster, field, dither = f.split('_')
        data = pyl.genfromtxt(f, delimiter='\t', names=True, dtype=None)
        try:
            for fiber, z, Q in zip(data['Fiber'], data['Redshift'],
                                   data['Quality']):
                if Q == 0:
                    r.append((field, 'D' + str(dither.rstrip('.results')),
                              fiber, z))
        except TypeError:
            fiber = int(data['Fiber'])
            z = float(data['Redshift'])
            Q = int(data['Quality'])
            if Q == 0:
                r.append(
                    (field, 'D' + str(dither.rstrip('.results')), fiber, z))

    print len(r), 'objects read'
    return r

 def draw_plot(self):
     global LAST_CALL, LAST_MODIFIED_DATE, IMAGE_BUFFER
     data_filename = "stats.txt"
     try:
         mtime = os.path.getmtime(data_filename)
     except OSError:
         mtime = 0
     modified_date = datetime.fromtimestamp(mtime)
     if LAST_CALL == self.path and modified_date == LAST_MODIFIED_DATE:
         IMAGE_BUFFER.seek(0)
         return IMAGE_BUFFER
     LAST_CALL = self.path
     LAST_MODIFIED_DATE = modified_date
     data = pylab.genfromtxt(data_filename , delimiter=',', dtype=int)
     y_data = data[:, 0]
     x_data = data[:, 1]
     if self.op == 'game':
         y_data = y_data[-self.game_count:]
         x_data = x_data[-self.game_count:]
     pylab.plot(x_data, y_data, '-')
     # pylab.show()
     IMAGE_BUFFER = io.BytesIO()
     pylab.savefig(IMAGE_BUFFER, format='png')
     IMAGE_BUFFER.seek(0)
     # pylab.legend()
     # pylab.title("Title of Plot")
     # pylab.xlabel("X Axis Label")
     # pylab.ylabel("Y Axis Label")
     pylab.close()
     return IMAGE_BUFFER

def main():

    is_transparent = False 
    
    f = open("pi_data.txt","r")
    
    # this is a little different than normal becase of the complex data for the floquet stability
    # multipliers. When we use the "dtype" option we get a single array of tuples so slicing is a
    # little more awkward has to look like data[#][#] to get a single value NOT data[#,#].
    data = pl.genfromtxt(f,comments="e",dtype="complex,complex,float")
   
    eigs1 = pl.array([])
    eigs2 = pl.array([])
    A = pl.array([])
    
    for i,j in enumerate(data):
        eigs1 = pl.append(eigs1,j[0])
        eigs2 = pl.append(eigs2,j[1])
        A = pl.append(A,j[2])

    fig1, ax1 = pl.subplots(2,2,sharex=True)
    ax1[0,0].plot(A,[k.real for k in eigs1],color = "Black")
    ax1[1,0].plot(A,[k.imag for k in eigs1],color = "Black")
    ax1[0,1].plot(A,[k.real for k in eigs2],color = "Black")
    ax1[1,1].plot(A,[k.imag for k in eigs2],color = "Black")

    ax1[0,0].set_ylabel("Re[$\lambda_1$]",fontsize=25)
    ax1[1,0].set_ylabel("Im[$\lambda_1$]",fontsize=25)
    ax1[0,1].set_ylabel("Re[$\lambda_2$]",fontsize=25)
    ax1[1,1].set_ylabel("Im[$\lambda_2$]",fontsize=25)
    ax1[1,0].set_xlabel("$A$",fontsize=25)
    ax1[1,1].set_xlabel("$A$",fontsize=25)
    fig1.tight_layout()
    fig1.savefig("paper_A_vs_eigs.png",dpi=300,transparent=is_transparent)
    os.system("open paper_A_vs_eigs.png")

def readSparseSystem (filename):
    """

    Convert's iSAM library's output (when printing sparse 
    matrices) to a scipy sparse matrix

    Returns
    -------
    
    a sparse COO matrix, the Cholesky factor R of the 
    information matrix
    
    """
    f = open (filename, 'r')
    line = f.readline ()
    f.close ()
    dimStr = re.search ('[0-9]+x[0-9]+', line).group (0)
    dim = int (dimStr.split ('x')[0])

    data = pl.genfromtxt (filename)
    data = data[1:,:]
    rows = data[:,0].astype (int)
    cols = data[:,1].astype (int)
    vals = data[:,2]

    R = scipy.sparse.coo_matrix ((vals, (rows, cols)), shape=(dim,dim))
    return R

def main():
    f = open("final_position.txt","r")
    data = pl.genfromtxt(f,comments = "L")
    
    # need to get every other
    x = pl.array([])
    y = pl.array([])
    for i,j in enumerate(data[:-7,2]):
        if i%4 == 0:
            x = pl.append(x,data[i,4])
            y = pl.append(y,j)
    
    print(x)
    print(y)
    fit = np.polyfit(x,y,2)

    print(fit)
    
    #fited = fit[0]+fit[1]*x + fit[2]*x**2
    fited = np.poly1d(fit)
    print(fited)

    #pl.plot(pl.append(x,[.262,.264,.266]),fited(pl.append(x,[.262,.264,.266])),color="black")
    pl.scatter(x,y,color = "black")
    pl.xlabel("$A$",fontsize="30")
    pl.ylabel("$x$",fontsize="30")
    pl.savefig("fin_pts.png",transparent=True,dpi=300)
    
    os.system("open fin_pts.png")

def main():

    parser = argparse.ArgumentParser()
    parser.add_argument("--dir", action="store", dest="dir", type=str, required=True)
    inargs = parser.parse_args()

    os.chdir(os.path.expanduser("~/Data/EC/4DBlock/Old/" + inargs.dir))

    os.mkdir("Last")

    b_num = get_b_num()

    os.system("cp info.txt Last/")

    all = os.listdir(".")
    for i, j in enumerate(all):
        if "poindat" in j:
            curfile = open(j, "r")
            to_file = open("Last/" + j, "w")
            to_file.write(curfile.readline())
            to_file.write(curfile.readline())
            cur_dat = pl.genfromtxt(curfile, comments="Z")
            for l in range(b_num):
                to_file.write(
                    str(cur_dat[-b_num + l, 0])
                    + " "
                    + str(cur_dat[-b_num + l, 1])
                    + " "
                    + str(cur_dat[-b_num + l, 2])
                    + " "
                    + str(cur_dat[-b_num + l, 3])
                )
                to_file.write("\n")
            to_file.close()
            curfile.close()

def sheat_vs_tempature(ancl):

    if 'ssheat_data.txt' not in os.listdir('.'):
        print('Need specific heat data')
        os.system('say Need specific heat data')
    if 'temp_granular_sliced.txt' not in os.listdir('.'):
        print('Need granular tempature data')
        os.system('say Need granular tempature data')

    tempature_file = open('temp_granular_sliced.txt','r')
    sheat_file = open('ssheat_data.txt','r')

    # first line is labels
    tempature_labels = tempature_file.readline()
    tempature_plotting_data = pl.genfromtxt(tempature_file)
    tempature_arr = tempature_plotting_data[:,1]

    # first line is labels
    sheat_labels = sheat_file.readline()
    sheat_plotting_data = pl.genfromtxt(sheat_file)
    #first column sweep variables
    var_arr = sheat_plotting_data[:,0]
    # evergy_stuff is next coulumn
    delta_E_sqrd = sheat_plotting_data[:,1]
    s_heat_arr = sheat_plotting_data[:,2]

    fig = pl.figure()
    ax = fig.add_subplot(111)
    # form of errorbar(x,y,xerr=xerr_arr,yerr=yerr_arr)
    pl.scatter(tempature_arr,s_heat_arr,c='k')
    #pl.errorbar(var_arr,averages_2,yerr=std_arr,c='b',ls='none',fmt='o')
    ax.set_xlabel(r'T_g',fontsize=30)
    ax.set_ylabel('Specific heat per particle',fontsize=20)
    fig.tight_layout()
    fig.savefig('T_vs_s_heat.png',dpi=300)
    pl.close(fig)


    print('\a')
    os.system('say finnished plotting tempature against specific heat')

def spatio_temporal(ancl):

    os.mkdir('SpatioTemporalVels')

    print('RIGHT NOW THIS IS ONLY FOR VX!!!!!!!')

    p_arr = pl.arange(0,ancl.N)
    

    # How many cycles do we want to look at?
    how_many = 10

    var_arr = pl.array([])
    for i,j in enumerate(os.listdir('.')):
        if 'poindat.txt' not in j:
            continue
        print('working on file ' + j)
        poin_num = int(j[:j.find('p')])
        cur_file = open(j,'r')
        cur_sweep_var = float(cur_file.readline().split()[-1])
        cur_data=pl.genfromtxt(cur_file)
        cur_file.close()

        var_arr = pl.append(var_arr,cur_sweep_var)
        
        count = 0
        grid = cur_data[-int(how_many*2.0*pl.pi/ancl.dt):,:ancl.N]

        # in 1D because particles never cross eachother we can order them in the images to mathch
        # their physical order.
        grid_ordered = pl.zeros(pl.shape(grid))
        # can just use the initial conditions to figure out where each is
        init_x = cur_data[0,ancl.N:2*ancl.N]
        sorted_x = sorted(init_x)
        for a,alpha in enumerate(sorted_x):
            for b,beta in enumerate(init_x):
                if alpha == beta:
                    grid_ordered[:,a]=grid[:,b]
        
    
        print('shape of grid_ordered: ' + str(pl.shape(grid_ordered)))
        
        fig = pl.figure()
        ax = fig.add_subplot(111)
        # form of errorbar(x,y,xerr=xerr_arr,yerr=yerr_arr)
        ax.imshow(grid_ordered,interpolation="nearest", aspect='auto')
        ax.set_xlabel('Particle',fontsize=30)
        #ax.set_aspect('equal')
        ax.set_ylabel(r'$ t $',fontsize=30)
        fig.tight_layout()
        fig.savefig('SpatioTemporalVels/%(number)04d.png'%{'number':poin_num})
        pl.close(fig)

    def __init__(self, xml_tree):
        self.root = xml_tree.getroot()
        self.misc = self.root.find('Misc')

        self.duration = self.misc.find('Duration')
        self.final_t = float(self.duration.attrib['Time'])

        self.pressure = self.misc.find('Pressure')
        self.pressure_avg = float(self.pressure.attrib['Avg'])
        self.pressure_tensor = self.pressure.find('Tensor')
        stream = io.BytesIO()
        stream.write(self.pressure_tensor.text)
        stream.seek(0)
        self.pressure_tensor = pylab.genfromtxt(stream)

def get_fit(which):
    f = open("final_position.txt","r")
    data = pl.genfromtxt(f,comments = "L")

    if which=="x":
        datnum = 2
    if which=="vx":
        datnum = 0

    x = pl.array([])
    y = pl.array([])
    for i,j in enumerate(data[:-7,datnum]):
        if i%2 == 0:
            x = pl.append(x,data[i,4])
            y = pl.append(y,j)

    fit = pl.polyfit(x,y,2)

    fitted = pl.poly1d(fit)
    return fitted 

def loadPlayByPlay(csvfile, vbose=0):
    skeys = ['game_id','type','playerName','posTeam','awayTeam','homeTeam']
    ikeys = ['seas','igame_id','dwn','ytg','yfog','yds']
    fkeys = []

    lines = [l.strip() for l in open(csvfile).readlines()]
    hd = lines[0]
    ks = hd.split(',')
    dt = []
    for k in ks:
        # postgres copy to file makes headers lower-case; this is a kludge
        if k=='playername':
            k = 'playerName'
        elif k=='posteam':
            k = 'posTeam'
        elif k=='away_team':
            k = 'awayTeam'
        elif k=='awayteam':
            k = 'awayTeam'
        elif k=='home_team':
            k = 'homeTeam'
        elif k=='hometeam':
            k = 'homeTeam'

        if k in skeys:
            tmp = (k, 'S16')
        elif k in ikeys:
            tmp = (k, 'i4')
        else:
            tmp = (k, 'f8')

        if vbose>=1:
            print k, tmp

        dt.append(tmp)
    dt = pylab.dtype(dt)
    data = pylab.genfromtxt(csvfile, dtype=dt, delimiter=',', skip_header=1)
    return data

def readCsv(ifile):
    skeys = ['date', 'homeTeam', 'awayTeam', 'game_id','player','posteam','oldstate','newstate']
    ikeys = ['seas','igame_id','yds']
    fkeys = []

    dt = []
    lines = [l.strip() for l in open(ifile).readlines()]
    hd = lines[0]
    ks = hd.split(',')
    for k in ks:
        if k in skeys:
            tmp = (k, 'S64')
        elif k in ikeys:
            tmp = (k, 'i4')
        elif k in fkeys:
            tmp = (k, 'f4')
        else:
            tmp = (k, 'f8')
        dt.append(tmp)

    dt = pylab.dtype(dt)
    data = pylab.genfromtxt(ifile, dtype=dt, skip_header=1, delimiter=',')
    return data

def main(cluster):
    cluster = cluster + '_r_mosaic.fits'
    f = pyl.figure(1, figsize=(8, 8))

    gc = aplpy.FITSFigure(cluster, north=True, figure=f)
    gc.show_grayscale(stretch='arcsinh')
    gc.set_theme('publication')

    gc.set_tick_labels_format(xformat='hh:mm:ss', yformat='dd:mm:ss')
    #gc.set_tick_labels_size('small')

    data = pyl.genfromtxt('./../analysis_all/redshifts/' +\
            cluster.split('_')[0]+'_redshifts.csv', delimiter=',', names=True,
            dtype=None)

    try:
        # filter out the specz's
        x = pyl.isnan(data['Specz'])
        # draw the specz's
        gc.show_markers(data['ra'][~x], data['dec'][~x], edgecolor='#ffbf00',
                        facecolor='none', marker='D', s=200)

    except ValueError:
        print 'no Speczs found'

    # draw observed but not redshifted
    x = data['Q'] == 2
    gc.show_markers(data['ra'][x], data['dec'][x], edgecolor='#a60628',
                    facecolor='none', marker='s', s=150)

    # draw redshifted
    x = (data['Q'] == 0) | (data['Q'] == 1)
    gc.show_markers(data['ra'][x], data['dec'][x], edgecolor='#188487',
                    facecolor='none', marker='o', s=150)

    pyl.tight_layout()
    pyl.show()

def crunchZfile(f,aCol,sCol,bCol,normFactor):
    '''
    Takes a zAveraged... data file generated from the crunchData
    function of this library and produces the arithmetic mean
    as well as the standard error from all seeds.  The error
    is done through the propagation of errors as:
    e = sqrt{ \sum_k (c_k e_k)^2 } where e_k are the individual
    seed's standard errors and c_k are the weighting coefficients
    obeying \sum_k c_k = 1.
    '''
    avgs,stds,bins = pl.genfromtxt(f, usecols=(aCol, sCol, bCol),
            unpack=True, delimiter=',')

    # get rid of any items which are not numbers..
    # this is some beautiful Python juju.
    bins = bins[pl.logical_not(pl.isnan(bins))]
    stds = stds[pl.logical_not(pl.isnan(stds))]
    avgs = avgs[pl.logical_not(pl.isnan(avgs))]

    # normalize data.
    stds *= normFactor
    avgs *= normFactor

    weights = bins/pl.sum(bins)

    avgs *= weights
    stds *= weights  # over-estimates error bars

    stds *= stds

    avg = pl.sum(avgs)
    stdErr = pl.sum(stds)

    stdErr = stdErr**0.5

    return avg, stdErr

source.add( composition )

# run simulation
sim.setShowProgress(True)
sim.run(source, 20000, True)


# ### (Optional) Plotting

# In[2]:

get_ipython().magic(u'matplotlib inline')
import pylab as pl
# load events
output.close()
d = pl.genfromtxt('events.txt', names=True)

# observed quantities
Z = pl.array([chargeNumber(id) for id in d['ID'].astype(int)])  # element
A = pl.array([massNumber(id) for id in d['ID'].astype(int)])  # atomic mass number
lE = pl.log10(d['E']) + 18  # energy in log10(E/eV))

lEbins = pl.arange(18, 20.51, 0.1)  # logarithmic bins
lEcens = (lEbins[1:] + lEbins[:-1]) / 2  # logarithmic bin centers
dE = 10**lEbins[1:] - 10**lEbins[:-1]  # bin widths

# identify mass groups
idx1 = A == 1
idx2 = (A > 1) * (A <= 7)
idx3 = (A > 7) * (A <= 28)
idx4 = (A > 28)

def make_int_c_sqrd_plot(v):

    qq,dt,beta,A,cycles,N,x_num_cell,y_num_cell,order,sweep_str,Dim  = of.get_system_info()
    
    y_lbl = r'$\int C^2(\tau)$'
    x_lbl = sweep_str

    sweep_var_arr = pl.array([])
    int_c_sqrd_arr = pl.array([])
    # loop over all files
    for i,j in enumerate(os.listdir('.')):
        if 'poindat.txt' not in j:
            continue
        
        work_file = open(j,'r')
        sweep_var_arr = pl.append(sweep_var_arr,float(work_file.readline().split()[-1]))
        data = pl.genfromtxt(work_file)
        work_file.close()

        average_out = pl.array([])

        # average over every particle in the simulation
        for a in range(N):
            
            if v == 'x':
                input_arr = data[:,Dim*N+a]
            if v == 'vx':
                input_arr = data[:,a]
            # Will only get into these if asking for 2D stuff anyway
            if v == 'y':
                if Dim==1:
                    print('No y in 1D')
                    quit()
                input_arr = data[:,(Dim+1)*N+a]
            if v == 'vy':
                if Dim==1:
                    print('No y in 1D')
                    quit()
                input_arr = data[:,(Dim-1)*N+a]

            print('shape of input_arr: ' + str(pl.shape(input_arr)))

            # lets try this for t from 0 to 100 cycles
            #t_arr = pl.arange(0,int(100.0*2.0*pl.pi),int(2.0*pl.pi))
            tau_arr = pl.arange(1,int(cycles*2.0*pl.pi),1)
            output_arr = pl.array([])
            for i,j in enumerate(tau_arr):
                cur = acf(j,input_arr)
                #print('current run of acf is (should be one number): ' +str(cur))
                
                output_arr = pl.append(output_arr,cur)

            # for average acf plot
            if a == 0:
                average_out = pl.append(average_out,output_arr)
            else:
                average_out += output_arr

        average_out = average_out/a
        print('shape of average_out (should be 1D array): ' + str(pl.shape(average_out)))

        # Romberg Integration integrate averages**2 to get int c^2 number
        #average_out = average_out[:257]
        #print('shape of :2**13+1 averag out: ' +str(len(average_out)))
        #int_c_sqrd = scint.romb(average_out**2,show=True)

        # simpson integration
        int_c_sqrd = scint.simps(average_out, x=None, dx=1, axis=-1, even='avg')
        int_c_sqrd_arr = pl.append(int_c_sqrd_arr, int_c_sqrd)
        print('int_c_sqrd_arr (should be number): ' + str(int_c_sqrd_arr))

         
            
    fig = pl.figure()
    ax = fig.add_subplot(111)
    #ax.set_xlim(x_rng)
    #ax.set_ylim(y_rng)
    ax.set_xlabel(x_lbl,fontsize=30)
    ax.set_ylabel(y_lbl,fontsize=30)
    ax.scatter(sweep_var_arr,int_c_sqrd_arr,c='k')
    fig.tight_layout()
    fig.savefig('int_c_sqrd_vs_sweep.png')
    pl.close(fig)

    print('\a')
    print('\a')
    os.system('open int_c_sqrd_vs_sweep.png')

    gc.show_grayscale(stretch='arcsinh', pmin=1, pmax=99.9)
    gc.set_tick_labels_format(xformat='hh:mm:ss', yformat='dd:mm:ss')
    gc.set_theme('publication')
    gc.set_tick_labels_size('small')

    # now for the axis lables
    if not i % 3 == 0:
        gc.axis_labels.hide_y()
    if i < 6:
        gc.axis_labels.hide_x()

    ax = fig.axes[-1]
    ax.set_title(cluster)

    data = pyl.genfromtxt('./../analysis_all/redshifts/' +\
            cluster.split('_')[0]+'_redshifts.csv', delimiter=',', names=True,
            dtype=None)

    try:
        # filter out the specz's
        x = pyl.isnan(data['Specz'])
        # draw the specz's
        gc.show_markers(data['ra'][~x], data['dec'][~x], edgecolor='#ffbf00',
                        facecolor='none', marker='D', s=50)

    except ValueError:
        print 'no Speczs found'

    # draw observed but not redshifted
    x = data['Q'] == 2
    gc.show_markers(data['ra'][x], data['dec'][x], edgecolor='#a60628',

import pylab
import time
import os
from matplotlib import dates
import datetime


os.environ["DISPLAY"] = ":0.0"

plotdailylogfilename = '/home/pi/Desktop/powersystem-logfiles/' + time.strftime("%Y-%m-%d") + '.log'
#dt=pylab.dtype([('f0','datetime64[ms]'),('f1',str,(2)),('f2',float,(5))])
dt=pylab.dtype([('f0','datetime64[ms]'),('s1',str,(15)),('f2',float,(1)),('f3',float,(1)),('f4',float,(1)),('f5',float,(1)),('f6',float,(1)),('f7',float,(1)),('f8',float,(1)),('f9',float,(1)),('f10',float,(1)),('f11',float,(1)),('s12',str,(15)),('f13',float,(1)),('f14',float,(1)),('f15',float,(1)),('f16',float,(1)),('f17',float,(1)),('f19',float,(1)),('f20',float,(1)),('f21',float,(1)),('f22',float,(1)),('f23',float,(1)),('f24',float,(1)),('f25',float,(1)),('f26',float,(1)),('f27',float,(1)),('f28',float,(1)),('f29',float,(1)),('f30',float,(1)),('f31',float,(1)),('f32',float,(1)),('f33',float,(1)),('f34',float,(1))])
data = pylab.genfromtxt(plotdailylogfilename,dtype=dt,delimiter='\t')

pylab.plot( data['f0'], data['f2'], label="Output")
pylab.plot( data['f0'], data['f3'], label="Battery")
pylab.legend()
pylab.title("Charge Data for " + time.strftime("%Y-%m-%d"))
pylab.xlabel("Time")
pylab.ylabel("Voltage: \n")
pylab.savefig("/var/www/html/voltagegraph.png")
pylab.clf()

pylab.plot( data['f0'], data['f21'], label="Amp")
pylab.plot( data['f0'], data['f22'], label="Watts")
pylab.legend()
pylab.title("Charge Data for " + time.strftime("%Y-%m-%d"))
pylab.xlabel("Time")
pylab.ylabel("Charge")
pylab.legend(loc='lower left')
pylab.savefig("/var/www/html/whcgraph.png")

#http://matplotlib.org/users/pyplot_tutorial.html
from matplotlib import pyplot;
from pylab import genfromtxt;  
mat0 = genfromtxt("data0.txt");
mat1 = genfromtxt("data1.txt");
pyplot.plot(mat0[:,0], mat0[:,1], label = "data0");
pyplot.plot(mat1[:,0], mat1[:,1], label = "data1");
pyplot.legend();
pyplot.show();