#  Generates summary reports on description txt files in present folder
#  There are three output files, A, B and C.
#  A produces a report with titles and column descriptions along with event data
#  B only produces event data, no titles; otherwise, same as A
#  C produces a report with titles and amplitude stats; otherwise same as A

# Things to do
# Re-examine magnitude philosophy on segments and for first 30 seconds.
# print DST year via dictionary entry

# sentuser txt (track data) contents on each line
# 1. Frame count relative to the trigger time.  Negative values are before trigger.  Positive values are after trigger.
# 2. Number of pixels that are above threshold.
# 3. Amplitude.  Basically, the sum of all pixels above threshold.
# 4. X coordinate of the event centroid, in pixels.
# 5. Y coordinate of the event centroid, in pixels.
#
# There is a header line, which is followed by an empty line

# Note that some of the functions here are not used here  They are
# part of a class, and as such are perfectly acceptable Python's reuseability
# philosophy. They produce statistics that may be needed in the future,
# particularly the segment calcs.  See the CollectTrkData class.

import sys, os, glob
import string
from numpy import *
from datetime import datetime, timedelta
import time
from scipy import stats as stats # scoreatpercentile

# DST for several years from
# Naval Obs. <http://aa.usno.navy.mil/faq/docs/daylight_time.php>
DST_dict = { # West coast, 8 hours from Greenwich for PST (Pacific States)
         2007:("2007/03/11 02:00:00", "2007/11/04 02:00:00"), 
         2008:("2008/03/09 02:00:00", "2008/11/02 02:00:00"),
         2009:("2009/03/08 02:00:00", "2009/11/01 02:00:00"),
         2010:("2010/03/14 02:00:00", "2010/11/07 02:00:00"),
         2011:("2011/03/13 02:00:00", "2011/11/06 02:00:00"),
         2012:("2012/03/11 02:00:00", "2012/11/04 02:00:00")}

def adjust_DST(DT_stamp, spring, fall):
    # yyyy/mm/dd hh:mm:ss in, adjusted out
    format = '%Y/%m/%d %H:%M:%S'
    dt      = datetime(*(time.strptime(DT_stamp, format)[0:6])) # get six tuple
    dspring = datetime(*(time.strptime(spring, format)[0:6]))
    dfall   = datetime(*(time.strptime(fall, format)[0:6]))
    if ((dt >= dspring) and (dt <= dfall)):
        #print "   mid-year adjustment"
        adj = timedelta(seconds = 3600) 
        dt = dt + adj       
    else:
        pass
        #print "   no adjustment"
    format = '%Y/%m/%d %H:%M:%S'
    return dt.strftime(format)

def adjust_ftime(afilename, sec):
    # yyyymmdd_hhmmss & seconds in, new yyyy/mm/dd hh:mm:ss with new time out
    td_stamp = afilename
    format = '%Y%m%d_%H%M%S'
    d = datetime(*(time.strptime(td_stamp, format)[0:6])) # get six tuple
    diff = timedelta(seconds = sec)
    d = d + diff
    format = '%Y/%m/%d %H:%M:%S'
    return d.strftime(format)

global output_fileA

class Track_data:
    def __init__(self, line):
        self.aline = line
        return

    def get_data(self, line):
        self.aline = line
        self.frame_count = int(self.aline[0:4])
        self.pix_above_thold = int(self.aline[5:11])
        self.amplitude = float(self.aline[14:19])
        self.cntrd_x = float(self.aline[20:27])
        self.cntrd_y = float(self.aline[28:35])
        data = [self.frame_count, self.pix_above_thold, \
                self.amplitude, self.cntrd_x, self.cntrd_y]
        return data

    def print_me(self):
        print self.amplitude, self.cntrd_x,self.cntrd_y
        return
        
class CollectTrkData:
    # Collect track data into a list of 3-tuples, xya. Misc track
    # calculations.
    def __init__(self):
        self.xya = []   # (x , y, amplitude)
        self.ntrk_pts = 0
        return

    def collection(self, amp, x, y):
        self.xya.append((x, y, amp)) # list of 3-tuples
        self.ntrk_pts += 1
        return
    
    def get_collection(self):
        return self.xya

    def distance(self):
        # total distance from start to end of track
        x0 = self.xya[0][0]
        y0 = self.xya[0][1]
        xn = self.xya[self.ntrk_pts-1][0]
        yn = self.xya[self.ntrk_pts-1][1]
        #print "wtw nframes: ", self.ntrk_pts
        #print "wtw-xy-0: ", x0,y0, "xy-n: ",xn,yn
        dist = sqrt((xn-x0)**2+(yn-y0)**2)
        return dist

    def angle(self):       
        x0 = self.xya[0][0]
        y0 = self.xya[0][1]
        xn = self.xya[self.ntrk_pts-1][0]
        yn = self.xya[self.ntrk_pts-1][1]
        x = xn - x0
        y = yn - y0
        a = (arctan2((xn-x0), (yn-y0))*57.3) -90.0 + 360.0
        angle = fmod(a,360.0)     
        return angle
    
    def print_sums(self):
        print self.total_amp, self.total_cntrd_x, \
              self.total_cntrd_y, self.ntrk_pts

    def print_xya(self, begin, end):
        for n in range(begin, end):
            print n, self.xya[n][0], self.xya[n][1]

    def get_amplitude_stats(self): # Simple track statistics
        abc = array(self.xya)[:,2]
        per_tile25 = stats.scoreatpercentile(abc,25.0)
        per_tile50 = stats.scoreatpercentile(abc,50.0)
        per_tile75 = stats.scoreatpercentile(abc,75.0)
        mean       = stats.mean(abc)
        std        = stats.std(abc)
        min        = abc.min()
        max        = abc.max()
        #print "per_tile: ", per_tile10, per_tile50, per_tile90
        #print "mean, std, min, max: ", mean,std,min,max
        trk_stats = (min, max, mean, std, per_tile25, per_tile50, per_tile75)
        return trk_stats

    def get_pix_speed(self):
        # Need to revise for degrees/sec when get ra/dec code available
        # calculate avg pixel/sec traveled every 10 frames, and
        # in segment group if less than 10
        # seg_size is the # of points in the segment, 0-9, 9-18, ...
        global output_fileA
        seg_size = 10
        n_segments = (self.ntrk_pts - 1)/(seg_size - 1)
        print "No. of line segments: ",n_segments
        segs = []
        for j in xrange(0,n_segments):
            segs.append(0)             
        if (self.ntrk_pts - 1 - n_segments*(seg_size - 1)) != 0:
            #print "hey: ", self.ntrk_pts - 1, n_segments*(seg_size - 1)
            # extra segment with less than seg_size frames
            n_segments = n_segments + 1
            n_in_last = self.ntrk_pts - 1 - (seg_size - 1)*(n_segments-1)
            #print "lastly:", n_in_last
        x = []; y = []
        # isolate x,y pairs
        for k in xrange(0, self.ntrk_pts):
            x.append(self.xya[k][0])
            y.append(self.xya[k][1])
            #print x[k],y[k]
        #print "x len: ",len(x)
        # create speed segments
        speed = []
        for j in xrange(0,n_segments):
            speed.append(0)
        
        #print "speed len", len(speed)
        # compute distances from line segment start to end   
        lower_index = 0
        bump = seg_size - 1
        upper_index = bump
        j = 0
        while (j < n_segments):
            #print "j,low,up: ",j, lower_index, upper_index
            # calculate segment distance
            #print "j low, up are: ",j, lower_index,upper_index
            # note this does not measure sub-segments
            x0 = x[lower_index]; y0 = y[lower_index]
            xn = x[upper_index]; yn = y[upper_index]
            #print "x0,y0 and x,y", x0,y0," -- ", xn,yn
            dist = sqrt((xn-x0)**2+(yn-y0)**2)
            #print "dist = ", dist
            speed[j] = dist/((1.0/30)*(upper_index - lower_index)) # (distance / time
            print "Pix/sec [",j,"]:", speed[j], "(",(upper_index - lower_index+1), ")",\
            ' Pixel distance: ',dist
            output_fileA.write("Pix/sec ["+str(j)+"]:"+str(speed[j])+ " ("+str((upper_index - lower_index+1))+")"+
            ' Pixel distance: '+str(dist)+'\n')
            lower_index = upper_index
            upper_index = lower_index + bump
            if upper_index > self.ntrk_pts-1:
                upper_index = self.ntrk_pts - 1
            #print "quick j,n_seg: ", j, n_segments
            j = j + 1

    def get_avg_xamp(self):
        # Find average amp in segments of size seg_size frames
        seg_size = 10
        n_segments = self.ntrk_pts/seg_size
        segs = []
        sample_size = []
        for j in xrange(0,n_segments):
            segs.append(0)
            sample_size.append(seg_size)
            
        if (n_segments*n_segments - self.ntrk_pts) != 0:
            # extra segment with less than seg_size frames
            n_segments = n_segments + 1
            sample_size.append(0)
            n_samples = self.ntrk_pts - seg_size*(n_segments - 1)
            sample_size[n_segments - 1] = n_samples
            segs.append(0)
        lower_index = 0
        for j in xrange(0, n_segments):
            upper_index = sum(sample_size[0:j+1])
            amps = []
            # calculate segment stats
            for k in xrange(0, self.ntrk_pts):
                amps.append(self.xya[k][2])
            segs[j] = sum(amps[lower_index:upper_index])
            lower_index = upper_index
            avg = segs[j]/sample_size[j]
            print "avg. amp sample size: ", j, sample_size[j]

def finish():
    print; print "Bye"
    print
    raw_input('Press Enter to Quit')
    sys.exit()

def get_track_data(file_name): #NEW
    # Collecting frame data
    track_file = open(file_name)
    track_file.readline();
    track_file.readline() # get past lines 1 and 2 to data
    # Get past the next 30 lines (ignore 1 sec of data). Start at track 0
    # An event starts at -1, but it is not possible to tell in the next
    # 60 seconds where it might end, since the untrigger value is unknown
    # Here I'm going to use 5 as an untrigger value for testing
    print
    print "wtw NEW",file_name
    for j in range(0,30):  
        line=track_file.readline()
    line_data = Track_data(line)   # Create object for line
    trk_data = CollectTrkData()
    untrigger = 5
    untrigger_switch = False
    n_data_pts = 0
    while 1:
        line = track_file.readline()  # track data lines
        if len(line) == 0:  # reach eof?
            break
        else:
            amp = int(line[5:11])
            if amp <= untrigger:
                if untrigger_switch:
                    # two untriggers in a row
                    break
                else:
                    # First untrigger event 
                    untrigger_switch = True
            fields = line_data.get_data(line)
            trk_data.collection(fields[2],fields[3],fields[4])
            n_data_pts += 1
            untrigger_switch = False
    print "wtw n:",n_data_pts
    xya = trk_data.get_collection() # (x, y, amplitude)
    
    # Now do  stats from xya, percentiles, range, std. dev., avg
    amp_stats = trk_data.get_amplitude_stats()
    fmt_amp_stats = "%8.1f %8.1f %8.1f %8.1f    %8.1f   %8.1f   %8.1f" % amp_stats
    #stats = trk_data.get_point_track_stats(n_data_pts, file_name)
    track_file.close()
    if n_data_pts == 0: 
        return
    fmt_basic_stats = "%3i " % trk_data.ntrk_pts\
        + "  %7.2f" % (trk_data.ntrk_pts*(1.0/30))\
        + " %7.2f" % trk_data.distance()
    #amp_stats = trk_data.get_amplitude_stats()
    #fmt_amp_stats = "%8.1f %8.1f %8.1f %8.1f    %8.1f   %8.1f   %8.1f" % amp_stats
    return (fmt_basic_stats, fmt_amp_stats)

def XXXget_track_data(file_name):
    global output_fileA
    track_file = open(file_name)
    track_file.readline(); abc = track_file.readline() # get past lines 1 and 2 to data
    for j in range(0,90):  # Get past the next 90 lines (ignore 3 sec of data)
        # Toss it, but maybe want to start collecing at one second
        abc=track_file.readline() 
    #print "wtw after 90: ",abc
    line = track_file.readline()
    trk_data = CollectTrkData()
    if len(line) == 0: # track has 90 or less points. short, maybe noise
        fini_stats = "%3i " % 0\
            + "  %7.2f" % (0.0)  # 3 seconds
        # effecitvely ended output for track
        return fini_stats
    while 1:
        if len(line) == 0:  # reach eof?
            break
        else:
            line_data = Track_data(line)
            fields = line_data.get_data()
            trk_data.collection(fields[2],fields[3],fields[4])
        line = track_file.readline()

    xya = trk_data.get_collection() # (x, y, amplitude)
    
    fmt_basic_stats = "%3i " % trk_data.ntrk_pts\
        + "  %7.2f" % (trk_data.ntrk_pts*(1.0/30))\
        + " %7.2f" % trk_data.distance()
    # Now do magnitude (amplitude) stats from xya, percentiles, range, std. dev., avg_stats = trk_data.get_amplitude_stats()
    fmt_amp_stats = "%8.1f %8.1f %8.1f %8.1f    %8.1f   %8.1f   %8.1f" % amp_stats
    track_file.close()
    return (fmt_basic_stats, fmt_amp_stats)

# Main program code
global output_fileA
print; print "Starting"; print
FLinit = raw_input('Enter your two initials in caps: ')
FLinit = string.capitalize(FLinit[0])+string.capitalize(FLinit[1])
print
print "UTC date-times are adjusted for the years 2006-2012"
print "   Note: The first 3.0 seconds of a track are ignored."
print "         Blank data for an event means either a short"
print "         track or a noisy first 3.0 seconds."
print "   date-time stamp is set about 1 second into the track."
print
#pathname = os.path.dirname(sys.argv[0])
now = datetime.now()
dt  = now.strftime("%Y%m%d_%H%M%S")
output_fileA = open('SentEventsA-'+FLinit+dt+'.txt','w') # Titles and events
output_fileB = open('SentEventsB-'+FLinit+dt+'.txt','w') # Events only
output_fileC = open('SentEventsC-'+FLinit+dt+'.txt','w') # Titles, events, amp stats
filenames = glob.glob('v2*.txt')
n_stop = 30000 # debug, testing -- to shorten output
seqno = 0
header_ctrl = 0
events_per_page = 20
dt_local = ""
first_page = True
for file in filenames:
    if mod(header_ctrl, events_per_page) == 0:
        if first_page == True:
            first_page = False
        else:
            output_fileA.write(chr(12))  # page/form feed
            output_fileC.write(chr(12))  # page/form feed
        header1 = "                                                  <---Track Stats--->"
        header2 = "Date/Time &    Station UTC                  Seq # Frames Time Span Pix Dst"
        header3 = "------------------- -- -------------------  ----- ------ --------- -------"
        print header1; print header2; print header3
        output_fileA.write(header1+"\n")
        output_fileA.write(header2+"\n")
        output_fileA.write(header3+"\n")
        amp_banner = "                    <--- Amplitude Statistics --->"
        amp_header = "  min      max      mean     std       pct-tile25 pct_tile50 pct_tile75"
        amp_underline= "  -------- -------- -------- --------  ---------- ---------- ----------"
        output_fileC.write(header1+amp_banner+"\n")
        output_fileC.write(header2+amp_header+"\n")
        output_fileC.write(header3+amp_underline+"\n")
    header_ctrl += 1
    seqno = seqno + 1

    dtf = file[1:16]   #date-time stamp on file 
    if dt_local == "":  # Identify if event is in a group
        dt = datetime.strptime(dtf,"%Y%m%d_%H%M%S")
        dt_last = dt.strftime("%Y/%m/%d %H:%M:%S")
        dt_local = dt_last
    else:
        dt = datetime.strptime(dtf,"%Y%m%d_%H%M%S")
        dt_local = dt.strftime("%Y/%m/%d %H:%M:%S")
    year   = int(str(dt_local)[0:4])
    spring = DST_dict[year][0]
    fall   = DST_dict[year][1]
    dt_utc = adjust_DST(dt_local, spring, fall)
    tz_adjust = 60*60.0*8  # PST
    dt_us = dt_utc
    dt_underscore = dt_us[0:10] + "_" + dt_us[11:]
    dt = datetime.strptime(dt_us,"%Y/%m/%d %H:%M:%S")
    dt_underscore = dt.strftime("%Y%m%d_%H%M%S")
    dt_utc    = adjust_ftime(dt_underscore, tz_adjust) # add 8 hours of time
    
    event_dt = str(dt_local) + " "+FLinit+" "+str(dt_utc)+("   %4i" % seqno)+"   " 
    output_fileA.write(event_dt+"")  # May need space after dt at some future time
    output_fileB.write(event_dt+"")
    output_fileC.write(event_dt+"")
    print event_dt,  # Output to shell window or console

    # Basic stats and amplitude stats
    trk_stats = get_track_data(file)  # adds track data to output
    output_fileA.write(trk_stats[0]+""+"\n")  # "" might need space in future
    output_fileB.write(trk_stats[0]+""+"\n")
    output_fileC.write(trk_stats[0]+"   "+trk_stats[1]+""+"\n")  # + amplitude stats
    print trk_stats[0]  # Output to shell window or console

    dt_last = dt_local    
    if seqno > n_stop: # debug, testing
        break

output_fileA.close()
output_fileB.close()
output_fileC.close()
finish()
