rpirena/multiplicities.py

#! /usr/bin/python

import sys, getopt, fileinput
import numpy as np
import matplotlib.pyplot as plt
options,files = getopt.getopt(sys.argv[1:], "V:vqm:D:o:", ["verbose="])

max_multiplicity = 10
max_dtime = 65535
verbose = 1
out = sys.stdout

# 0: errors
# 1: warnings
# 7: parsed items
# 8: valid input lines
# 9: all input

def logger(priority, message):
    if priority <= verbose:
        sys.stderr.write(message)

for o,v in options:
    if o=="-V" or o=="--verbose":
        verbose = int(v)
    if o=="-v":
        verbose += 1
    if o=="-q":
        verbose -= 1
    if o=="-m":
        max_multiplicity = int(v, 0)
    if o=="-D":
        max_dtime = int(v, 0)
        if max_dtime<0:
            raise ValueError("max_dtime must be positive %d" % max_dtime)
        if max_dtime>0xffff:
            logger(1, "max_dtime = 0xffff")
            max_dtime = 0xffff
    if o=="-o":
        out = file(v, "w")  
        
def analyse_burst(burst, Time, Pressure):
    ll = [l for b,l in burst]
    burst = [b for b,l in burst]
    if not burst:
        return
    logger(4, "burst of length %d at %d\n" % (len(burst), Time))
    #for i in range(1, len(burst)-1):
        #for j in range(2, len(burst)+1-i):
            #dti = sum(burst[i:i+j])
            #if dti <= max_dtime:
                #out.write("D %d %.1f %d %d\n" % (Time, Pressure, dti, j))
    l = len(burst)
    eburst = l*[max_dtime] + burst + l*[max_dtime]
    for i in range(l):
        dtis=[]
        for j in range(1, min(max_multiplicity, l)):
            dti = min([sum(eburst[i+l-k+1:i+l-k+1+j]) for k in range(j+1)])
            dtis.append(dti)
            if dti <= max_dtime:
                out.write("DT %d %.1f %d %d\n" %(Time, Pressure, dti, j+1))
        out.write(("EIDT "+" ".join(ll[i][1:])+len(dtis)*" %d"+"\n") % tuple(dtis))

burst = []
for line in fileinput.input(files):
   logger(8 if line[0] in "HEP" else 9, line)
   if line[:2]=="H ":
        Time = float(line.split()[1])
        logger(7, "Time: %d\n" % Time)
   if line[:3]=="EI ":
        ll = line.split()
        dtime = int(ll[4])
        if dtime>=max_dtime:
            analyse_burst(burst, Time, 0.0)
            burst = []
        burst.append((dtime,ll))


       
#plt.figure()
#plt.plot(hist(dti), mult)
#plt.show()


""" np.savetxt("dti.txt", dti)
dti = np.loadtxt("dti.txt", unpack=True)
 """