% TESTFALSERATE Test variation of false rate with threshold
%
% TESTFALSERATE tests the dependence of measured white noise false
% event rate on significance threshold, signal space, and maximum
% mismatch parameter.
%
% TESTFALSERATE repeatedly applies the Q transform to ideal stationary
% white noise and records all tiles exceeding a significance
% corresponding to a roughly 10 Hz white noise false event rate.
%
% As it runs, TESTFALSERATE writes the normalized energies of events
% to a log file testfalserate/testfalserate.log.  Upon completion,
% TESTFALSERATE also writes these normalized energies results to a MAT
% file testfalserate/testfalserate.mat and produces the following
% plot.
%
%   testfalserate/testfalserate_1.png
%     Measured and predicted false event rates vs. snr threshold
%
% See also QTILE, QTRANSFORM, QTHRESHOLD. and QSELECT.

% Shourov K. Chatterji
% shourov@ligo.caltech.edu

% $Id: testfalserate.m,v 1.1 2007/04/12 05:42:34 shourov Exp $

% path to qtransform code base
path('..', path);

% create output directory
unix('mkdir -p testfalserate');

% open log file
logFileName = 'testfalserate/testfalserate.log';
logFileID = fopen(logFileName, 'w');

% number of trials
numberOfTrials = 100;

% signal space parameters
timeRange = 64;
qRange = [sqrt(11) 100];
frequencyRange = [0 Inf];
sampleFrequency = 4096;

% desired tiling accuracy
maximumMismatch = 0.2;

% tile the signal space
tiling = qtile(timeRange, qRange, frequencyRange, sampleFrequency, ...
               maximumMismatch);

% time vector for signal construction
time = 0 : 1 / sampleFrequency : timeRange - 1 / sampleFrequency;

% initialize random number generator
rand('state', sum(100 * clock));
randn('state', sum(100 * clock));

% initialize results vectors
normalizedEnergies = [];

% begin loop over trials
for trialNumber = 1 : numberOfTrials,

  % construct noise
  noise = normrnd(0, sqrt(tiling.sampleFrequency / 2), 1, length(time));

  % transform to one-sided frequency domain signal
  noise = fft(noise);
  noise = noise(1 : length(noise) / 2 + 1);

  % apply q transform
  calibration = [];
  outlierFactor = 2.0;
  transform = qtransform(noise, tiling, calibration, outlierFactor);

  % identify significant events
  falseEventRate = 10;
  events = qthreshold(transform, tiling, 0, falseEventRate);

  % down select overlapping events
  events = qselect(events);

  % report results
  fprintf(logFileID, '%9.3e\n', events{1}.normalizedEnergy);

  % append result to results vectors
  normalizedEnergies = [normalizedEnergies; ...
                        events{1}.normalizedEnergy(:)];

% end loop over trials
end

% save data
save testfalserate/testfalserate.mat tiling normalizedEnergies

% determine measured and predicted event rates
sortedEnergies = sort(normalizedEnergies);
livetime = (tiling.duration - 2 * tiling.transientDuration) * ...
           numberOfTrials;
measuredRate = (length(normalizedEnergies) : -1 : 1) / livetime;
predictedRate = exp(-sortedEnergies) * 1.5 * ...
                tiling.numberOfIndependents / tiling.duration;
snrs = sqrt(2 * sortedEnergies);

% display measured and predicted event rates
clf;
set(gca, 'FontSize', 18);
handle = semilogy(snrs, predictedRate, 'b.-', ...
                  snrs, measuredRate, 'r.-');
axis([floor(min(snrs)) ceil(max(snrs)) ...
      10^floor(log10(1/livetime)) 10]);
set(gca, 'YTick', 10.^(floor(log10(1/livetime)) : 1));
axis([floor(min(snrs)) ceil(max(snrs)) ...
      10^ceil(log10(1/livetime)) 10]);
set(gca, 'YTick', 10.^(ceil(log10(1/livetime)) : 1));
grid on;
set(handle, 'MarkerSize', 16);
set(handle, 'LineWidth', 2);
xlabel('Signal to noise ratio threshold');
ylabel('False event rate');
title('Measured vs. predicted false event rates');
legend('Predicted', 'Measured', 3);
drawnow;
print -dpng testfalserate/testfalserate_1.png
unix(['convert -resize 500x500 testfalserate/testfalserate_1.png ' ...
      'testfalserate/testfalserate_1_thumbnail.png']);