results_rfi

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  function  [...] = results_rfi(xxx)

   Feb 1, 2012 (MAS) - Created.

   Generates the output plots.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

0001 function r = results_rfi(startMjd,endMjd,basedir)    
0002 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0003 %  function  [...] = results_rfi(xxx)
0004 %
0005 %   Feb 1, 2012 (MAS) - Created.
0006 %
0007 %   Generates the output plots.
0008 %
0009 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0010 
0011 
0012 % Read in the statistical flagging file.
0013 statFile = open([basedir 'stat_flag_results.mat']);
0014 
0015 disp([basedir 'stat_flag_results.mat']);
0016 %disp('hello');
0017 
0018 statUtc = statFile.statUtcVec;
0019 statFlag = statFile.statFlagArray_f;
0020 
0021 %disp('hello');
0022 
0023 
0024 % Read in the candidates file.
0025 candidateFile = open([basedir 'rfi_candidates.mat']);
0026 
0027 %disp('hello');
0028 
0029 candidateUtc = candidateFile.candidateUtc;
0030 candidateLen = candidateFile.candidateLen;
0031 candidateSet_f = candidateFile.candidateSet_f;
0032 candidateSetA_f = candidateFile.candidateSetA_f;
0033 candidateSetZ_f = candidateFile.candidateSetZ_f;
0034 candidateSet_p = candidateFile.candidateSet_p;
0035 candidateSetA_p = candidateFile.candidateSetA_p;
0036 candidateSetZ_p = candidateFile.candidateSetZ_p;
0037 candidateVer = candidateFile.candidateVer;
0038 
0039 %disp('hello');
0040 
0041 % Read in the recommended parameters.
0042 paramFile = open([basedir 'params.mat']);
0043 pArray_f = paramFile.pArray_f;
0044 pArray_p = paramFile.pArray_p;
0045 
0046 
0047 % How many parameter sets to we have?
0048 nSets = size(pArray_f,2);
0049 
0050 
0051 % Create the RFI plot directories.
0052 rfiplotdir = cell(nSets,2);
0053 for k=1:nSets
0054     
0055     rfiplotdir(k,1) = {[basedir 'rfiplots_filtered_' int2str(k) '/']};
0056     
0057     if ~exist(char(rfiplotdir(k,1)),'dir')
0058         mkdir(char(rfiplotdir(k,1)));
0059     end
0060     
0061     
0062     rfiplotdir(k,2) = {[basedir 'rfiplots_polonly_' int2str(k) '/']};
0063     
0064     if ~exist(char(rfiplotdir(k,2)),'dir')
0065         mkdir(char(rfiplotdir(k,2)));
0066     end
0067 end
0068 
0069 disp('hello');
0070 
0071 % Loop over and produce plots of flagged events...
0072 %================================
0073 
0074 canPlot =  (candidateVer ~= 0) | ...
0075     (max(candidateSet_f,[],2) > 0) | (max(candidateSet_p,[],2) > 0);
0076 
0077 canUtcPlot = candidateUtc(canPlot);
0078 canLenPlot = candidateLen(canPlot);
0079 canSetFPlot = candidateSet_f(canPlot,:);
0080 canSetPPlot = candidateSet_p(canPlot,:);
0081 canVerPlot = candidateVer(canPlot);
0082 
0083 
0084 
0085 % How many time sections will we use?
0086 nTimes = round(10 * (endMjd - startMjd));
0087 
0088 
0089 % Initialize the inner iterator:
0090 nPlot = sum(canPlot);
0091 
0092 l = 0;
0093 while (1)% l=1:nPlot
0094     l = l + 1;
0095     
0096     % Make sure we haven't already plotted these fellas...
0097     filePre = [char(rfiplotdir(1,2)) mjd2string(canUtcPlot(l))];
0098     if ~(exist([filePre '_R.png'],'file') || ...
0099             exist([filePre '_F.png'],'file') || ...
0100             exist([filePre '_M.png'],'file') || ...
0101             exist([filePre '_O.png'],'file'))
0102         break;
0103     end
0104     
0105     if (l == nPlot)
0106         l = l + 1;
0107         break;
0108     end
0109 end
0110 
0111 
0112 
0113 % Loop over the time segments.
0114 disp('Looping over the time segments...');
0115 for k=1:nTimes
0116 
0117     if l > nPlot
0118         disp('All plots generated...');
0119         break;
0120     end
0121     
0122 
0123     startSec = mjd2string(startMjd + 0.1 * (k-1) - 5/24/60);
0124     endSec = mjd2string(startMjd + 0.1 * k + 5/24/60);
0125     
0126     disp(['Data from ' startSec ' to ' endSec]);
0127 
0128     
0129     % Make sure we have got something here, before reading in...
0130     if canUtcPlot(l) >= tstr2mjd(endSec)
0131         disp('No plots for this data range...');
0132         continue;
0133     end
0134         
0135     
0136     % Read in the data, but just the interesting parts...
0137     [d, pointFlag, galacticFlag, sunFlag, moonFlag] = ...
0138         read_posflag(startSec,endSec);
0139     
0140     utcVec = d.antenna0.receiver.utc;
0141     I1Vec = d.antenna0.receiver.data(:,1);
0142     QVec = 0.5*(d.antenna0.receiver.data(:,2) + ...
0143         d.antenna0.receiver.data(:,4));
0144     UVec = 0.5*(d.antenna0.receiver.data(:,3) - ...
0145         d.antenna0.receiver.data(:,5));
0146     I2Vec = d.antenna0.receiver.data(:,6);
0147     
0148 %    dataLen = length(utcVec);
0149 
0150 
0151     % Loop over candidates.
0152     while (l <= nPlot && (canUtcPlot(l)+max(10.5*canLenPlot(l),3000)/8640000) < tstr2mjd(endSec))
0153 
0154         disp(['Candidate ' int2str(l) ' out of ' int2str(nPlot)]);
0155 
0156         canBuffer = max(10*canLenPlot(l),3000);
0157         
0158         
0159         % Candidate Starting Index:
0160         canStart = max(1,ceil((canUtcPlot(l) - utcVec(1)) * 8640000 - canBuffer));
0161         % Candidate Ending Index:
0162         canEnd = min(length(utcVec),floor((canUtcPlot(l) - utcVec(1)) * 8640000 + canBuffer));
0163 %        % Candidate Centre Index:
0164 %        canCen = round((candidateUtc(l) - utcVec(1)) * 8640000) + 1;
0165         
0166 %        disp(length(I1Vec));
0167 %        disp(canEnd);
0168         I1Med = median(I1Vec(canStart:canEnd));
0169         I1Range = range(I1Vec(canStart:canEnd));
0170         QMed = median(QVec(canStart:canEnd));
0171         QRange = range(QVec(canStart:canEnd));
0172         UMed = median(UVec(canStart:canEnd));
0173         URange = range(UVec(canStart:canEnd));
0174         I2Med = median(I2Vec(canStart:canEnd));
0175         I2Range = range(I2Vec(canStart:canEnd));
0176         
0177         maxRange = max(max([I1Range QRange URange I2Range]),5);
0178         [~, fStart] = min(abs(statUtc - utcVec(canStart)));
0179         [~, fEnd] = min(abs(statUtc - utcVec(canEnd)));
0180 
0181         % The Time vector.
0182         tVec = (utcVec(canStart:canEnd) - canUtcPlot(l)) * 24 * 3600;
0183         
0184 
0185         % Now produce a plot for each of the parameter sets...
0186         for m=1:nSets
0187 
0188             if canVerPlot(l) == 0 && ...
0189                     canSetFPlot(l,m) == 0 && canSetPPlot(l,m) == 0
0190                 continue
0191             end
0192 
0193             figure('visible','off');
0194 
0195             % Generate the axes...
0196             axis([min(tVec), max(tVec), -1.5*maxRange, 1.5*maxRange]);
0197             hold on;
0198 
0199 
0200             % Plot the position flags.
0201             flagplot(tVec,sunFlag(canStart:canEnd),[0.8 0.7 0.7],maxRange);
0202             flagplot(tVec,moonFlag(canStart:canEnd),[0.7 0.7 0.8],maxRange);
0203             flagplot(tVec,galacticFlag(canStart:canEnd),[0.7 0.8 0.7],maxRange);
0204             flagplot(tVec,pointFlag(canStart:canEnd),[0.7 0.7 0.7],maxRange);
0205 
0206 
0207             % Plot the candidate.
0208             plot(tVec,I1Vec(canStart:canEnd) - I1Med + maxRange,'b');
0209             plot(tVec,QVec(canStart:canEnd) - QMed + maxRange/3,'r');
0210             plot(tVec,UVec(canStart:canEnd) - UMed - maxRange/3,'g');
0211             plot(tVec,I2Vec(canStart:canEnd) - I2Med - maxRange,'m');
0212 
0213             % Show the statistical flagging.
0214             plot(86400*(statUtc(fStart:fEnd)-canUtcPlot(l)), ...
0215                 maxRange*(statFlag(fStart:fEnd,m)-0.5), 'k');
0216 
0217 
0218             legend('Sun', 'Moon', 'Galaxy', 'Point Sources', ...
0219                 'I1', 'Q', 'U', 'I2', 'Flag');
0220             xlabel('Time (s)');
0221             ylabel('Power (s)');
0222             title(mjd2string(canUtcPlot(l)));
0223 
0224             hold off;
0225 
0226             
0227             % First for the FILTERED flagging.
0228             if canVerPlot(l) == 1 && canSetFPlot(l,m) == 1
0229                 outfilename = [char(rfiplotdir(m,1)) mjd2string(canUtcPlot(l)) '_R.png']; % Caught
0230             elseif canVerPlot(l) == 1 && canSetFPlot(l,m) == 0
0231                 outfilename = [char(rfiplotdir(m,1)) mjd2string(canUtcPlot(l)) '_M.png']; % Missed
0232             elseif canVerPlot(l) == 0 && canSetFPlot(l,m) == 1
0233                 outfilename = [char(rfiplotdir(m,1)) mjd2string(canUtcPlot(l)) '_F.png']; % False Positive
0234             else
0235                 outfilename = [char(rfiplotdir(m,1)) mjd2string(canUtcPlot(l)) '_O.png']; % Other
0236             end
0237 
0238             print('-dpng', outfilename);
0239 
0240             % Then for the POLONLY flagging.
0241             if canVerPlot(l) == 1 && canSetPPlot(l,m) == 1
0242                 outfilename = [char(rfiplotdir(m,2)) mjd2string(canUtcPlot(l)) '_R.png']; % Caught
0243             elseif canVerPlot(l) == 1 && canSetPPlot(l,m) == 0
0244                 outfilename = [char(rfiplotdir(m,2)) mjd2string(canUtcPlot(l)) '_M.png']; % Missed
0245             elseif canVerPlot(l) == 0 && canSetPPlot(l,m) == 1
0246                 outfilename = [char(rfiplotdir(m,2)) mjd2string(canUtcPlot(l)) '_F.png']; % False Positive
0247             else
0248                 outfilename = [char(rfiplotdir(m,2)) mjd2string(canUtcPlot(l)) '_O.png']; % Other
0249             end
0250 
0251             print('-dpng', outfilename);
0252 
0253             close(gcf);
0254         end
0255 
0256         
0257         l = l + 1;
0258 
0259     end
0260     
0261 end
0262 
0263 
0264 % XXX
0265 
0266 %
0267 % % Loop over the time segments.
0268 % disp('Looping over the time segments...');
0269 % for k=1:nTimes
0270 %
0271 %
0272 % %for k=1:nPlot
0273 %
0274 %     canBuffer = max(10*canLenPlot(k),3000);
0275 %
0276 %     canMjdStart = canUtcPlot(k) - canBuffer/8640000.;
0277 %     canMjdEnd = canUtcPlot(k) + canBuffer/8640000.;
0278 %
0279 %     canStrStart = mjd2string(canMjdStart);
0280 %     canStrEnd = mjd2string(canMjdEnd);
0281 %
0282 %     disp(['RFI candidate from ' canStrStart ' to ' canStrEnd]);
0283 %
0284 %     filePre = [char(rfiplotdir(1)) mjd2string(canUtcPlot(k))];
0285 %     if exist([filePre '_R.png'],'file') || ...
0286 %             exist([filePre '_F.png'],'file') || ...
0287 %             exist([filePre '_M.png'],'file') || ...
0288 %             exist([filePre '_O.png'],'file')
0289 %         disp('This event already plotted');
0290 %         continue;
0291 %     end
0292 %
0293 %     % Read in the data:
0294 %     [d, pointFlag, galacticFlag, sunFlag, moonFlag] = ...
0295 %         read_posflag(canStrStart,canStrEnd);
0296 %
0297 %     % The time vector.
0298 %     tVec = (d.antenna0.receiver.utc - canUtcPlot(k)) * 86400;
0299 %
0300 %     I1Vec = d.antenna0.receiver.data(:,1);
0301 %     QVec = 0.5*(d.antenna0.receiver.data(:,2) + ...
0302 %         d.antenna0.receiver.data(:,4));
0303 %     UVec = 0.5*(d.antenna0.receiver.data(:,3) - ...
0304 %         d.antenna0.receiver.data(:,5));
0305 %     I2Vec = d.antenna0.receiver.data(:,6);
0306 %
0307 %     I1Med = median(I1Vec);
0308 %     I1Range = range(I1Vec);
0309 %     QMed = median(QVec);
0310 %     QRange = range(QVec);
0311 %     UMed = median(UVec);
0312 %     URange = range(UVec);
0313 %     I2Med = median(I2Vec);
0314 %     I2Range = range(I2Vec);
0315 %
0316 %     maxRange = max(max([I1Range QRange URange I2Range]),5);
0317 %     [~, fStart] = min(abs(statUtc - canMjdStart));
0318 %     [~, fEnd] = min(abs(statUtc - canMjdEnd));
0319 %
0320 %
0321 %     % Now produce a plot for each of the parameter sets...
0322 %     for l=1:nSets
0323 %
0324 %         if canVerPlot(k) == 0 && canSetPlot(k,l) == 0
0325 %             continue
0326 %         end
0327 %
0328 %         figure('visible','off');
0329 % %        figure();
0330 %
0331 %         % Generate the axes...
0332 %         axis([min(tVec), max(tVec), -1.5*maxRange, 1.5*maxRange]);
0333 %         hold on;
0334 %
0335 %
0336 %         % Plot the position flags.
0337 %         flagplot(tVec,sunFlag,[0.8 0.7 0.7],maxRange);
0338 %         flagplot(tVec,moonFlag,[0.7 0.7 0.8],maxRange);
0339 %         flagplot(tVec,galacticFlag,[0.7 0.8 0.7],maxRange);
0340 %         flagplot(tVec,pointFlag,[0.7 0.7 0.7],maxRange);
0341 %
0342 %
0343 %         % Plot the candidate.
0344 %         plot(tVec,I1Vec - I1Med + maxRange,'b');
0345 %         plot(tVec,QVec - QMed + maxRange/3,'r');
0346 %         plot(tVec,UVec - UMed - maxRange/3,'g');
0347 %         plot(tVec,I2Vec - I2Med - maxRange,'m');
0348 %
0349 %         % Show the statistical flagging.
0350 %         plot(86400*(statUtc(fStart:fEnd)-canUtcPlot(k)), ...
0351 %             maxRange*(statFlag(fStart:fEnd,l)-0.5), 'k');
0352 %
0353 %
0354 %         legend('Sun', 'Moon', 'Galaxy', 'Point Sources', ...
0355 %             'I1', 'Q', 'U', 'I2', 'Flag');
0356 % %        legend('Galaxy', 'Point Sources', ...
0357 % %            'I1', 'Q', 'U', 'I2', 'Flag');
0358 %         xlabel('Time (s)');
0359 %         ylabel('Power (s)');
0360 %         title(mjd2string(canUtcPlot(k)));
0361 %
0362 %         hold off;
0363 %
0364 %         if canVerPlot(k) == 1 && canSetPlot(k,l) == 1
0365 %             outfilename = [char(rfiplotdir(l)) mjd2string(canUtcPlot(k)) '_R.png']; % Caught
0366 %         elseif canVerPlot(k) == 1 && canSetPlot(k,l) == 0
0367 %             outfilename = [char(rfiplotdir(l)) mjd2string(canUtcPlot(k)) '_M.png']; % Missed
0368 %         elseif canVerPlot(k) == 0 && canSetPlot(k,l) == 0
0369 %             outfilename = [char(rfiplotdir(l)) mjd2string(canUtcPlot(k)) '_F.png']; % False Positive
0370 %         else
0371 %             outfilename = [char(rfiplotdir(l)) mjd2string(canUtcPlot(k)) '_O.png']; % Other
0372 %         end
0373 %
0374 %         print('-dpng', outfilename);
0375 %
0376 %     end
0377 % end
0378 
0379 
0380 
0381 % Loop over the parameter sets...
0382 %================================
0383 
0384 
0385 
0386 % Prep the output arrays.
0387 falseRate = zeros(nSets,2);
0388 falseErr = zeros(nSets,2);
0389 catchRate = zeros(nSets,2);
0390 catchErr = zeros(nSets,2);
0391 
0392 
0393 for k=1:nSets
0394 
0395     disp(['Parameter set: ' int2str(k)]);
0396 
0397        
0398     for l=1:2
0399     
0400         switch l
0401             case 1
0402                 candidateSet = candidateSet_f;
0403                 candidateSetA = candidateSetA_f;
0404                 candidateSetZ = candidateSetZ_f;
0405             case 2
0406                 candidateSet = candidateSet_p;
0407                 candidateSetA = candidateSetA_p;
0408                 candidateSetZ = candidateSetZ_p;
0409         end
0410         
0411         % False positive rate.
0412         falseRate(k,l) = sum((candidateSet(:,k) & (candidateVer == 0)) .* ...
0413             candidateLen) / (24*3600*100*(endMjd-startMjd));
0414 
0415         falseA = sum((candidateSetA(:,k) & (candidateVer == 0)) .* ...
0416             candidateLen) / (24*3600*100*(endMjd-startMjd));
0417         falseZ = sum((candidateSetZ(:,k) & (candidateVer == 0)) .* ...
0418             candidateLen) / (24*3600*100*(endMjd-startMjd));
0419 
0420         falseErr(k,l) = abs((falseZ - falseA) / 2);
0421 
0422 
0423         % RFI catch rate.
0424         catchRate(k,l) = sum(candidateSet(:,k) & (candidateVer == 1)) ./ ...
0425             sum(candidateVer == 1);
0426 
0427         catchA = sum(candidateSetA(:,k) & (candidateVer == 1)) ./ ...
0428             sum(candidateVer == 1);
0429         catchZ = sum(candidateSetZ(:,k) & (candidateVer == 1)) ./ ...
0430             sum(candidateVer == 1);
0431 
0432         catchErr(k,l) = abs((catchZ - catchA) / 2);
0433     
0434     end
0435 end
0436 
0437 
0438 
0439 % Write out the statistics for the recommended parameters..
0440 save([basedir 'params.mat'], 'pArray_f', 'pArray_p', ...
0441     'falseRate', 'falseErr', 'catchRate', 'catchErr');
0442 
0443 
0444 
0445 
0446 
0447 
0448 %     % Plot the histogram vs. feature flags.
0449 %     fracFeat = zeros(length(features),1);
0450 %     for l=1:length(features)
0451 %
0452 %         nFeat = sum(candidateFeat == 2^features(l));
0453 %         nRFI = sum(statflagVec(featVec(:,k) == 1));
0454 %
0455 %         if nFeat > 0
0456 %             fracFeat(k) = nRFI / nFeat;
0457 %         else
0458 %             fracFeat(k) = 0;
0459 %         end
0460 %     end
0461 %
0462 % figure('Units', 'pixels', ...
0463 %     'Position', [100 100 500 375]);
0464 %
0465 % yMat = count(1:6,:);
0466 % figure; bar(yMat);
0467 %
0468 % bar(1:13,fracFeat);
0469 %
0470 % axis([0 14 0 0.02])
0471 %
0472 %
0473 % hold on;
0474 %
0475 %
0476 % title([startTime ' - ' endTime]);
0477 %
0478 % %hXLabel = xlabel('Feature');
0479 % ylabel('Fraction Flagged');
0480 % hold off;
0481 %
0482 % %set(hXLabel, ...
0483 % %    'Rotation'   , 45);
0484 %
0485 %
0486 % featureNames = {'source', 'elscan', 'calibrator', 'abscal', 'blank', ...
0487 %     'skydip', 'radio point cross', 'radio point scan', 'opt point', ...
0488 %     'beammap', 'noise', 'ncp', 'noise event'};
0489 %
0490 %
0491 % % Place the text labels
0492 % t = text(1:13,zeros(1,13),featureNames);
0493 % set(t,'HorizontalAlignment','left','VerticalAlignment','top', ...
0494 % 'Rotation',-30);
0495 %
0496 % set(gca, ...
0497 %   'XTick'       , 1:13, ...
0498 %   'XTickLabel'       , {});
0499 %
0500 % set(gcf, 'PaperPositionMode', 'auto');
0501 % print('-dpng',[histdir startTime '-' endTime '.png']);
0502 % close;
0503 %
0504 %
0505 %
0506 % % Plots of flagged data.
0507 %
0508 % % Plots of missed data.
0509 %
0510 % % Plots of histogram of flags vs. feature flags.
0511 
0512 
0513 
0514 r = 1;
0515 end
0516 
0517 
0518 
0519 %%%%%%%%%%%%%%%%%%%%%%%%%
0520 function flagplot(tVec,flagVec,flagColour,maxRange)
0521 
0522 
0523 % Do we even have any such flags...?
0524 if max(flagVec) > 0
0525 
0526 
0527     % Identify the flagged regions...
0528     dFlag = diff(flagVec);
0529 
0530     flagStart = find(dFlag > 0);
0531     flagEnd = find(dFlag < 0);
0532 
0533     if isempty(flagStart)
0534         flagStart = 1;
0535     end
0536     if isempty(flagEnd)
0537         flagEnd = length(flagVec);
0538     end
0539     if flagEnd(1) < flagStart(1)
0540         flagStart = cat(1,1,flagStart);
0541     end
0542     if flagEnd(end) < flagStart(end)
0543         flagEnd = cat(1,flagEnd,length(flagVec));
0544     end
0545 
0546 
0547     nFlag = length(flagStart);
0548 
0549     % Prepare the vertices for the filled region...
0550     xVec = zeros(4*nFlag,1);
0551     yVec = zeros(4*nFlag,1);
0552     for k=1:nFlag
0553         xVec(4*k-3) = tVec(flagStart(k));
0554         xVec(4*k-2) = tVec(flagStart(k));
0555         xVec(4*k-1) = tVec(flagEnd(k));
0556         xVec(4*k) = tVec(flagEnd(k));
0557 
0558         yVec(4*k-3) = -1.5*maxRange;
0559         yVec(4*k-2) = 1.5*maxRange;
0560         yVec(4*k-1) = 1.5*maxRange;
0561         yVec(4*k) = -1.5*maxRange;
0562     end
0563 
0564 else
0565 
0566 
0567     xVec = [0 0];
0568     yVec = [-2*maxRange -2*maxRange];
0569 
0570 end
0571 
0572 
0573 fill(xVec,yVec,flagColour,'EdgeColor',flagColour);
0574 
0575 end