using PyPlot

# Create a type for Channels
type Channel
    name :: String
    times :: Vector
    times_clean :: Vector
    offset :: Vector
    offset_clean :: Vector
    offset_from_linear :: Vector
    drift_rate :: Float64
    t_zero_offset :: Float64
    sigma :: Float64
    sigma_from_linear :: Float64
end

# Specify channel names
channel_names = [
    # "recent-cesium-drift-longer"
    "Time Code Generator in MSR",
    "Time Code Translator in EX",
    "Time Code Translator in EY"
]
N = length(channel_names)

# Make shared axis histograms
print("Making shared axis histograms...")
histograms = figure("histograms",figsize=(10,10))
subplots_adjust(hspace=0.0) # Set the vertical spacing between axes
println(" done.")

subplot(311)
title("Histograms of 1PPS Offset in Cesium Clock Time Distribution System")

# Container for the channels
channels = Array(Channel,0)

# Do stuff for each channel
for i in 1:N
    
    name = channel_names[i]
    print("Handling channel $name...")

    # Load the offset data
    data = readdlm("$name.dat")
    data = float(data)
    x = data[:,1]    # GPS times in seconds
    y = data[:,2]    # offsets in seconds
    n = length(x)    # number of samples
    
    # Remove noise from when Dave and I were coopting the BNC cable in EY
    if name == "Time Code Translator in EY"
        print(" Munging data...")
        y = y[(x .< 1119123000) | (x .> 1119135000)]
        x = x[(x .< 1119123000) | (x .> 1119135000)]
        print(" done munging...")
    end
    
    # Remove outliers
    yclean = y[abs(y - mean(y)) .< 3std(y)]
    xclean = x[abs(y - mean(y)) .< 3std(y)]

    # Shift the times by 1.116e9 seconds; the linear
    # regression algorithm needs bigger relative
    # differences in the x values
    xclean1 = xclean - 1.116e9
    
    # Find the linear regression (slope will be correct)
    a, b = linreg(xclean1,yclean)
    
    # Find the offsets from the linear trend
    z = yclean .- [a+b*i for i in xclean1]
    z = float(z)
    
    # Remove outliers
    zclean = z[abs(z - mean(z)) .< 3std(z)]
    
    # Make a subplot for each histogram (courtesy of http://nbviewer.ipython.org/github/gizmaa/Julia_Examples/blob/master/pyplot_subplot.ipynb)
    subplot(100N + 10 + i)
    eval(parse("ax$i = gca()"))
    xticks(-1e-6:2e-7:2e-7)
    xlim(-1e-6,2e-7)
    PyPlot.plt.hist(yclean, label="Raw Offset for $name", bins=[-1e-6:1e-8:-6e-7])
    PyPlot.plt.hist(zclean, label="Linear Drift Removed", bins=[-2e-7:5e-9:2e-7])
    legend(loc="upper left",fancybox="true")
    # ax1[:set_xscale]("log") # Set the x axis to a logarithmic scale
    eval(parse("setp(ax$i[:get_xticklabels](),visible=false)")) # Disable x tick labels
    # grid("on")
    # yticks(0.1:0.2:0.9)
    # ylim(0.0,1.0)
    
    ylabel("Frequency")
    xlabel("Offset from aLIGO Timing System in Seconds")
    
    println(" done.")
    
    # Create Channel objects to hold information about each channel
    push!(channels, Channel(name, x, xclean, y, yclean, z, b, a, std(yclean), std(z)))
end
eval(parse("setp(ax$N[:get_xticklabels](),visible=true)")); # Enable x tick labels for last plot

# Make shared axis histograms
print("Making shared axis timeseries...")
timeseries = figure("timeseries",figsize=(10,15))
subplots_adjust(hspace=0.0) # Set the vertical spacing between axes
println(" done.")

subplot(311)
title("Timeseries of 1PPS Offset in Cesium Clock Time Distribution System")

# Do stuff for each channel
for i in 1:N
    
    channel = channels[i]
    name = channel.name
    print("Handling channel $name...")

    # Load the offset data
    x = channel.times    # GPS times in seconds
    y = channel.offset   # offsets in seconds
    n = length(x)        # number of samples
    a = channel.t_zero_offset  # y-intercept (with shifted x values)
    b = channel.drift_rate     # drift rate of atomic clock
    
    # Use the shifted x values we used to calculate the linear regression
    x1 = x - 1.116e9
    trend = [a+b*i for i in x1]
    trend = float(trend)
    
    # Make a subplot for each timeseries (courtesy of http://nbviewer.ipython.org/github/gizmaa/Julia_Examples/blob/master/pyplot_subplot.ipynb)
    subplot(100N + 10 + i)
    eval(parse("ax$i = gca()"))
    yticks(-9e-7:5e-8:-6.49e-7)
    ylim(-9.5e-7,-6e-7)
    plot(x, y, label="Offset for $name", color="green")
    plot(x, trend, label="Trend", color="red")
    legend(loc="upper right",fancybox="true")
    # ax1[:set_xscale]("log") # Set the x axis to a logarithmic scale
    eval(parse("setp(ax$i[:get_xticklabels](),visible=false)")) # Disable x tick labels
    # grid("on")
    # yticks(0.1:0.2:0.9)
    # ylim(0.0,1.0)
    
    xlabel("GPS Time (seconds)")
    ylabel("Offset from aLIGO Timing System (seconds)")
    
    println(" done.")
end
eval(parse("setp(ax$N[:get_xticklabels](),visible=true)")); # Enable x tick labels for last plot

for channel in channels
    println(channel.drift_rate)
end
@show mean_drift_rate = mean([channel.drift_rate for channel in channels])
@show delta_f_measured = 1078e-15
@show delta_f_calibrated = mean_drift_rate + delta_f_measured * (1 + mean_drift_rate)
@show what_to_enter_into_monitor2 = delta_f_calibrated / 1e-15