remove hardcode and update input params

README.txt

@@ -9,7 +9,7 @@ Run the compiled executable file on execute nodes.
 -------------------
 The following files need to be put in the same directory for compiling.
 Project scripts:
-    - search_O1.m (main script)
+    - run_search.m (main script)
     - viterbi_colFLT.m
     - argmax.m
     - binaryfilter.m
@@ -47,17 +47,3 @@ atlas-matlab-init <executable file name> <parameter1, parameter2, ...>
 EXAMPLE: (Run as Condor jobs)
 Refer to sample submission file viterbi_search.sub
 
-4. Notes
---------
-The main script search_O1.m is an example of the search. The parameters can 
-be adjusted according to the requirement of different searches. For the purpose
-of parallel search, the frequency band is separated into 1-Hz sub-bands and this script takes two parameters.
-
-Parameter 1 - start frequency for all the parallel jobs
-Parameter 2 - Job offset index
-
-EXAMPLE:
-./search_O1 100 0 - Search a 1-Hz sub-band from 100 to 101 Hz.
-./search_O1 100 1 - Search a 1-Hz sub-band from 101 to 102 Hz. 
-
-

binaryfilter.m

@@ -29,7 +29,7 @@ M = length(f);
 % Identify the indices to place the matched filter in the sub-band.
 % k indicates the interpolated indices at query points n/P, 
 % using nearest interpolation. Set extrapval to 1 for out of range values.
-k = interp1(f-mean(f),1:M,n/P,'nearest',1);
+k = interp1(f-f0,1:M,n/P,'nearest',1);
 % Initialise the matched filter for the whole sub-band
 b = zeros(M,1);
 % Place the filter

job_scripts/viterbi_search.sub

 Universe   = vanilla
 Executable = /home/ling.sun/O1/Viterbi/matlab/job_scripts/atlas-matlab-init
-Arguments  = /home/ling.sun/O1/Viterbi/matlab/job_scripts/search_O1 81 $(Process)
+Arguments  = /home/ling.sun/O1/Viterbi/matlab/job_scripts/run_search 81 $(Process) /home/ling.sun/O1/Viterbi/fstat/ /home/ling.sun/O1/Viterbi/output/outViterbi-$(Process).dat 1.805 0.8 0.05 68023.70496 13
 output     = out.$(Cluster).$(Process)
 error      = err.$(Cluster).$(Process)
 log        = log.$(Cluster).$(Process)
 notification    = Never
 accounting_group = ligo.prod.o1.cw.directedbinary.sideband
-Queue   20
+Queue   5

readfstats.m

-function [f, X, N, T] = readfstats(dirName, startFreq, steps)
+function [f, X, N, T] = readfstats(dirName, bandIdx, steps)
 % READFSTATS Read Fstat files for all steps
 %
 % Note
 % ----
 % Fstat file name should follow the format:
-%     Fstat-<start freq>-<step index>.dat
+%     Fstat-<freq band index>-<step index>.dat
 %     <step index> start from 0.
 % 
 % I/O Spec
 % ========
 % Input:
 %     dirName   - Directory of Fstat files
-%     startFreq - Start frequency marking the Fstat file name
+%     bandIdx   - Frequency band index marking the Fstat file name
 %     steps     - Total number of steps 
 % 
 % Output: 
@@ -22,7 +22,7 @@ function [f, X, N, T] = readfstats(dirName, startFreq, steps)
 
 X = [];
 for n = 1:steps
-    fileName = ['Fstat-' num2str(startFreq) '-' num2str(n-1) '.dat']
+    fileName = ['Fstat-' num2str(bandIdx) '-' num2str(n-1) '.dat']
 
 	try
 		fstatData = load([dirName fileName]);

search_O1.m → run_search.m

-function search_O1(freq,idx)
-% SEARCH_O1 This script is for O1 Sco X-1 search. 
-% The search is separated into parallel 1-Hz sub-jobs.
-% The start frequency for each sub-job is the overal start frequency plus job index.
+function run_search(bandIdx,bandOffset,fstatdir,outfile,a0median,a0uncert,a0delta,orbitPeriod,viterbiSteps)
+% RUN_SEARCH Run Viterbi tracking using Fstat (Sco X-1) 
+% 
+% Note
+% ----
+% The script is written for the convenience of running parallel sub-jobs.
+% Frequency band index for each sub-job is the input band index plus job offset.
+% (i.e. Fstat band index = bandIdx + bandOffset)
+% The input band index is only used to index Fstat file. 
+% The real frequencies in each band are read from Fstat.
 %
 % I/O Spec
 % ========
 % Input:
-%     freq - Start frequency for all the parallel jobs
-%     idx  - Condor job offset index
-%     Note: Start frequency for each sub-job: freq + idx
+%     bandIdx      - Start Fstat band index for all the parallel jobs
+%     freqOffset   - Sub-job offset 
+%     fstatdir     - Fstat file directory
+%     outfile      - Output file directory and name
+%     a0median     - Median value of asini searched
+%     a0uncert     - Uncertainty of asini (e.g. 0.8 stands for +/-80% uncertainty)
+%     a0delta      - asini uncertainty step size 
+%     orbitPeriod  - Binary orbit period
+%     viterbiSteps - Total Viterbi tracking steps
 %
 % Output:
-%     Save output files to the output directory.
+%     Save results to the specified output file.
 
 format long g
 
-startFreq = str2num(freq) + str2num(idx)
-
-% Set directories and file name
-fstatdir = '/home/ling.sun/O1/Viterbi/fstat/';
-outdir = '/home/ling.sun/O1/Viterbi/output/';
-outfile = strcat(outdir,'Viterbi-', num2str(startFreq),'.dat')
-
 % Set parameters
-a0median = 1.805     % median asini being searched
-g = [-0.8:0.05:0.8]; % asini grid (+/-80% uncertainty)
-P = 68023.70496      % orbit period
-steps = 13           % total steps
+fstatBandIdx = str2double(bandIdx) + str2double(bandOffset)
+P = str2double(orbitPeriod)      
+steps = str2double(viterbiSteps)
+
+% Set asini search grid
+a0median = str2double(a0median)
+grid = -1*str2double(a0uncert):str2double(a0delta):str2double(a0uncert);
+a0 = (1+grid)*a0median;
 
 % Read Fstats
 tic
-[f,X] = readfstats(fstatdir, startFreq, steps);
+[f,X] = readfstats(fstatdir, fstatBandIdx, steps);
 toc
 
-fmean = f(1)+0.5; % Take the mean value of each 1-Hz sub-band
-a0 = (1+g)*a0median; % asini values being searched
+% Take the mean value of each sub-band
+fmean = mean(f) 
 
-for m=1:length(g)
+for m=1:length(grid)
     disp(sprintf('a0=%g',a0(m)))
 
     tic
@@ -56,12 +65,12 @@ end
 % [maxsc, a0ind] = max(sc)
 % path = path0(a0ind,:)
 
-% Save output file for each 1-Hz sub-band
+% Save output file for each sub-band
 fid = fopen(outfile,'wt');
 fprintf(fid,'%%StartFreq\t%g\n',f(1));
 fprintf(fid,'%%Score\ta0\tPath\n');
 
-for m=1:length(g)
+for m=1:length(grid)
     fprintf(fid,'%g\t%g\t',sc(m),a0(m));
     for n=1:steps
         fprintf(fid,'%d\t',path0(m,n));