%% Time interval corresponding to ER7
er7 = [1117400416, 1118329216];

%% time periods
gps = [ 1117332016, 1117371616; ... % Jun 3 2015 02:00:00 UTC - Jun 3 2015 13:00:00 UTC
        1117479616, 1117488616; ... % Jun 4 2015 19:00:00 UTC - Jun 4 2015 21:30:00 UTC
        1117492216, 1117544416; ... % Jun 4 2015 22:30:00 UTC - Jun 5 2015 13:00:00 UTC
        1117571416, 1117591216; ... % Jun 5 2015 20:30:00 UTC - Jun 6 2015 02:00:00 UTC
        1117612816, 1117663216; ... % Jun 6 2015 08:00:00 UTC - Jun 6 2015 22:00:00 UTC
        1117668616, 1117699216; ... % Jun 6 2015 23:30:00 UTC - Jun 7 2015 08:00:00 UTC
        1117710016, 1117737016; ... % Jun 7 2015 11:00:00 UTC - Jun 7 2015 18:30:00 UTC
        1117749616, 1117767616; ... % Jun 7 2015 22:00:00 UTC - Jun 8 2015 03:00:00 UTC
        1117861216, 1117877416; ... % Jun 9 2015 05:00:00 UTC - Jun 9 2015 09:30:00 UTC
        1117886416, 1117897216; ... % Jun 9 2015 12:00:00 UTC - Jun 9 2015 15:00:00 UTC
        1117931416, 1117945816; ... % Jun 10 2015 00:30:00 UTC - Jun 10 2015 04:30:00 UTC
        1117985416, 1117989616; ... % Jun 10 2015 15:30:00 UTC - Jun 10 2015 16:40:00 UTC
        1118041216, 1118086216; ... % Jun 11 2015 07:00:00 UTC - Jun 11 2015 19:30:00 UTC
        1118302216, 1118338216; ... % Jun 14 2015 07:30:00 UTC - Jun 14 2015 17:30:00 UTC
        1120807817, 1120813217; ... % Jul 13 2015 07:30:00 UTC - Jul 13 2015 09:00:00 UTC
        1120887917, 1120893617; ... % Jul 14 2015 05:45:00 UTC - Jul 14 2015 07:20:00 UTC
        1121760017, 1121765417; ... % Jul 24 2015 08:00:00 UTC - Jul 24 2015 09:30:00 UTC
        1121898617, 1121927417; ... % Jul 25 2015 22:30:00 UTC - Jul 26 2015 06:30:00 UTC
        1121931017, 1121986817; ... % Jul 26 2015 07:30:00 UTC - Jul 26 2015 23:00:00 UTC
        1121990417, 1121994017; ... % Jul 27 2015 00:00:00 UTC - Jul 27 2015 01:00:00 UTC
        1122089417, 1122093917; ... % Jul 28 2015 03:30:00 UTC - Jul 28 2015 04:45:00 UTC
        1122096617, 1122129017; ... % Jul 28 2015 05:30:00 UTC - Jul 28 2015 14:30:00 UTC
        1122184817, 1122188417; ... % Jul 29 2015 06:00:00 UTC - Jul 29 2015 07:00:00 UTC
        1122195617, 1122219017; ... % Jul 29 2015 09:00:00 UTC - Jul 29 2015 15:30:00 UTC
        1122368417, 1122372017; ... % Jul 31 2015 09:00:00 UTC - Jul 31 2015 10:00:00 UTC
        1122382817, 1122391817; ... % Jul 31 2015 13:00:00 UTC - Jul 31 2015 15:30:00 UTC
        1122460217, 1122485417; ... % Aug 01 2015 10:30:00 UTC - Aug 01 2015 17:30:00 UTC
        1122543017, 1122552017; ... % Aug 02 2015 09:30:00 UTC - Aug 02 2015 12:00:00 UTC
        1122584417, 1122589817; ... % Aug 02 2015 21:00:00 UTC - Aug 02 2015 22:30:00 UTC
        ];
        
% round to multiple of 60 seconds!
gps = 60*round(gps/60);

% minimum range drop to consider it a glitch (Mpc)
drop = 5;

% detrending time (minutes)
detrend = 60;
    
%% loop over all periods
glitch_times = [];
for i=1:size(gps,1)
    % get range from data
    conn = nds2.connection('nds.ligo-wa.caltech.edu', 31200);
    data = conn.fetch(gps(i,1), gps(i,2), {'H1:DMT-SNSW_EFFECTIVE_RANGE_MPC.mean,m-trend'});
    range = data(1).getData();
    fs = data(1).getLength() / (gps(i,2)-gps(i,1));
    t = gps(i,1) + (0:length(range)-1)/fs;

    % detrend with a 10 minutes window
    wind = hanning(detrend);
    wind = wind / sum(wind);
    ranget = conv([ones(detrend,1) * range(1) ; range; ones(detrend,1) * range(end)], wind, 'same');
    ranget = ranget(detrend+1:end-detrend);

    % select large drops
    idx = find(range>0 & (range-ranget)<-drop);

    % make a figure
    f1 = figure();
    plot(t - gps(i,1), range, t - gps(i,1), ranget, 'm', t(idx) - gps(i,1), range(idx), 'or')
    xlabel(sprintf('Time starting from GPS %d [s]', gps(i,1)))
    ylabel('Range [Mpc]')
    legend('Range', 'Trend', 'Glitches', 'location', 'southwest')
    grid
    img = getframe(f1);
    imwrite(img.cdata, sprintf('range_drop_%d.png', gps(i,1)))

    fprintf(2, 'GPS: %d-%d: found %d glitches\n', gps(i,1), gps(i,2), length(idx));
    
    % append to the list of glitches
    glitch_times = [glitch_times, t(idx)];
end
close all
fprintf(2, 'Found %d glitches\n', length(glitch_times));
    

f = fopen('glitch_times.txt', 'w');
for i=1:numel(glitch_times)
    fprintf(f, '%d\n', glitch_times(i));
end
fclose(f);


%% Day time statistics

% this GPS time correspond to 12am local time
gps_12am = 1122879617;

% find out the time of the day of each glitch in hours
glitch_day_time = mod(glitch_times - gps_12am, 86400) / 3600;

% estimate the day coverage of the lock stretches
G = [];
for i=1:size(gps,1)
    G = [G, gps(i,1):gps(i,2)];
end
G = mod(G - gps_12am, 86400) / 3600;

% compute hourly coverage
h = 0.5:1:23.5;
hG = hist(G,h);
hG = hG/3600;

% compute glitch occurence per hour
hglitch = hist(glitch_day_time, h);
glitch_hourly_rate = hglitch ./ hG;

figure()
subplot(311)
bar(h, hG)
title(sprintf('All reasonable locks between June 1st and August 6th (total %d glitches)', sum(hglitch)))
xlim([0, 24])
xlabel('Hour of the day (LT)')
ylabel('Total hours covered by lock stretches')
set(gca, 'XTick', 0:24)
subplot(312)
bar(h, hglitch)
xlim([0, 24])
xlabel('Hour of the day (LT)')
ylabel('Number of glitches in the hour')
set(gca, 'XTick', 0:24)
subplot(313)
bar(h, hglitch./hG)
xlim([0, 24])
xlabel('Hour of the day (LT)')
ylabel('Glitch hourly rate')
set(gca, 'XTick', 0:24)
hold all
plot(h, ones(size(h)) .* mean(hglitch./hG), 'r-')
plot(h, ones(size(h)) .* (mean(hglitch./hG) + std(hglitch./hG)), 'r--')
plot(h, ones(size(h)) .* (mean(hglitch./hG) - std(hglitch./hG)), 'r--')
plot(h, ones(size(h)) .* (mean(hglitch./hG) + 2*std(hglitch./hG)), 'r:')
plot(h, ones(size(h)) .* (mean(hglitch./hG) - 2*std(hglitch./hG)), 'r:')


%% restrict only to ER7
er7_glitches = glitch_times(glitch_times > er7(1) & glitch_times < er7(2));
G = [];
for i=1:size(gps,1)
    G = [G, gps(i,1):gps(i,2)];
end
er7_G = G(G > er7(1) & G < er7(2));
er7_G = mod(er7_G - gps_12am, 86400) / 3600;
er7_glitches = mod(er7_glitches - gps_12am, 86400) / 3600;

% compute hourly coverage for ER7
h = 0.5:1:23.5;
er7_hG = hist(er7_G,h);
er7_hG = er7_hG/3600;
% compute glitch occurence per hour
er7_hglitch = hist(er7_glitches, h);
er7_glitch_hourly_rate = er7_hglitch ./ er7_hG;

figure()
subplot(311)
bar(h, er7_hG)
title(sprintf('ER7 (total %d glitches)', sum(er7_hglitch)))
xlim([0, 24])
xlabel('Hour of the day (LT)')
ylabel('Total hours covered by lock stretches')
set(gca, 'XTick', 0:24)
subplot(312)
bar(h, er7_hglitch)
xlim([0, 24])
xlabel('Hour of the day (LT)')
ylabel('Number of glitches in the hour')
set(gca, 'XTick', 0:24)
subplot(313)
bar(h, er7_hglitch./er7_hG)
xlim([0, 24])
xlabel('Hour of the day (LT)')
ylabel('Glitch hourly rate')
set(gca, 'XTick', 0:24)
hold all
plot(h, ones(size(h)) .* mean(er7_hglitch./er7_hG), 'r-')
plot(h, ones(size(h)) .* (mean(er7_hglitch./er7_hG) + std(er7_hglitch./er7_hG)), 'r--')
plot(h, ones(size(h)) .* (mean(er7_hglitch./er7_hG) - std(er7_hglitch./er7_hG)), 'r--')
plot(h, ones(size(h)) .* (mean(er7_hglitch./er7_hG) + 2*std(er7_hglitch./er7_hG)), 'r:')
plot(h, ones(size(h)) .* (mean(er7_hglitch./er7_hG) - 2*std(er7_hglitch./er7_hG)), 'r:')



%% all but ER7
ner7_glitches = glitch_times(glitch_times < er7(1) | glitch_times > er7(2));
G = [];
for i=1:size(gps,1)
    G = [G, gps(i,1):gps(i,2)];
end
ner7_G = G(G < er7(1) | G > er7(2));
ner7_G = mod(ner7_G - gps_12am, 86400) / 3600;
ner7_glitches = mod(ner7_glitches - gps_12am, 86400) / 3600;

% compute hourly coverage for ER7
h = 0.5:1:23.5;
ner7_hG = hist(ner7_G,h);
ner7_hG = ner7_hG/3600;
% compute glitch occurence per hour
ner7_hglitch = hist(ner7_glitches, h);
ner7_glitch_hourly_rate = ner7_hglitch ./ ner7_hG;

figure()
subplot(311)
bar(h, ner7_hG)
title(sprintf('All but ER7 (total %d glitches)', sum(ner7_hglitch)))
xlim([0, 24])
xlabel('Hour of the day (LT)')
ylabel('Total hours covered by lock stretches')
set(gca, 'XTick', 0:24)
subplot(312)
bar(h, ner7_hglitch)
xlim([0, 24])
xlabel('Hour of the day (LT)')
ylabel('Number of glitches in the hour')
set(gca, 'XTick', 0:24)
subplot(313)
bar(h, ner7_hglitch./ner7_hG)
xlim([0, 24])
xlabel('Hour of the day (LT)')
ylabel('Glitch hourly rate')
set(gca, 'XTick', 0:24)
hold all
plot(h, ones(size(h)) .* mean(ner7_hglitch./ner7_hG), 'r-')
plot(h, ones(size(h)) .* (mean(ner7_hglitch./ner7_hG) + std(ner7_hglitch./ner7_hG)), 'r--')
plot(h, ones(size(h)) .* (mean(ner7_hglitch./ner7_hG) - std(ner7_hglitch./ner7_hG)), 'r--')
plot(h, ones(size(h)) .* (mean(ner7_hglitch./ner7_hG) + 2*std(ner7_hglitch./ner7_hG)), 'r:')
plot(h, ones(size(h)) .* (mean(ner7_hglitch./ner7_hG) - 2*std(ner7_hglitch./ner7_hG)), 'r:')

%% get one minute of data around each drop
conn = nds2.connection('nds.ligo.caltech.edu', 31200);
    
for i=1:length(glitch_times)
    fprintf(2, '%d\n', i);
    % get DARM data
    data = conn.fetch(glitch_times(i), glitch_times(i)+60, {'H1:LSC-DARM_IN1_DQ'});
    e = data(1).getData();
    fs = length(e)/60;
    tt = (0:length(e)-1)/fs;
    
    % find the exact time using maximum of BLRMS
    [B,A] = butter(4, [30, 300]/(fs/2), 'bandpass');
    ee = filtfilt(B, A, double(e));
    ee = ee.^2;
    [m, j] = max(ee(fs:end-fs));
    j = j + fs;
    
    % high pass signal to detrend
    [B,A] = butter(4, 100/(fs/2), 'high');
    eh = filtfilt(B, A, double(e));
    % high pass signal to detrend
    [B,A] = butter(6, 1000/(fs/2), 'high');
    eh2 = filtfilt(B, A, double(e));
   
    % GPS time of the glitch
    tglitch(i) = glitch_times(i) + tt(j);
     
    % plot time series (+- deltaT seconds)
    deltaT = 0.1;
    scrsz = get(0,'ScreenSize');
    f2 = figure('Position',[1 scrsz(4)/2 scrsz(3)/3 scrsz(4)/1.2]);
    subplot(311)
    plot(tt, ee, tt(j), ee(j), 'or')
    ylim([0, 1.2*m])
    xlabel('Time [s]')
    ylabel('BLRMS 30-300 Hz')
    grid
    subplot(312)
    if j-fs*deltaT > 0 
        plot(tt(j-fs*deltaT:j+fs*deltaT), e(j-fs*deltaT:j+fs*deltaT), tt(j), e(j), 'or')
    else
        plot(tt(1:fs*2*deltaT), e(1:fs*2*deltaT), tt(j), e(j), 'or')
    end
    title(sprintf('Glitch #%d GPS %d (full signal)', i, glitch_times(i)))
    xlabel('Time [s]')
    ylabel('DARM\_IN1 [cts]')
    xlim([tt(j)-deltaT, tt(j)+deltaT])
    grid
    subplot(313)
    if j-fs*deltaT > 0 
        plot(tt(j-fs*deltaT:j+fs*deltaT), eh(j-fs*deltaT:j+fs*deltaT), tt(j), eh(j), 'or')
        hold all
        plot(tt(j-fs*deltaT:j+fs*deltaT), eh2(j-fs*deltaT:j+fs*deltaT), tt(j), eh2(j), 'or')
    else
        plot(tt(1:fs*2*deltaT), eh(1:fs*2*deltaT), tt(j), eh(j), 'or')
        hold all
        plot(tt(1:fs*2*deltaT), eh2(1:fs*2*deltaT), tt(j), eh2(j), 'or')
    end
    title(sprintf('Glitch #%d GPS %d', i, glitch_times(i)))
    xlabel('Time [s]')
    ylabel('DARM\_IN1 [cts]')
    title('High pass 100 Hz')
    xlim([tt(j)-deltaT, tt(j)+deltaT])
    grid    
    
    img = getframe(f2);
    imwrite(img.cdata, sprintf('glitch_plots/glitch_%d_darm.png', i))
    
    class(i) = input(sprintf(' %d / %d class? ', i, numel(glitch_times)));
    close all
    
    save classification2.mat class
end
