fold_eval

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  function r = fold_eval()

   24-May-2012 (MAS): Requires FILTERED or POLONLY data, or 
       CLASSIC data that has had stokes run on it.  But, honestly, just
       use FILTERED data.

   03-May-2012 (MAS): Now assumes that pipeline (alpha and load)
   has been run on the data.

   25-Jan-2012 (MAS): Modified to work within the automated tuning.

   01-Nov-2011 (MAS): Created.

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

0001 function foldResult = fold_eval(d)
0002 
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %  function r = fold_eval()
0005 %
0006 %   24-May-2012 (MAS): Requires FILTERED or POLONLY data, or
0007 %       CLASSIC data that has had stokes run on it.  But, honestly, just
0008 %       use FILTERED data.
0009 %
0010 %   03-May-2012 (MAS): Now assumes that pipeline (alpha and load)
0011 %   has been run on the data.
0012 %
0013 %   25-Jan-2012 (MAS): Modified to work within the automated tuning.
0014 %
0015 %   01-Nov-2011 (MAS): Created.
0016 %
0017 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0018 
0019 
0020 azRateErr = 2e-3;
0021 
0022 
0023 % Get the feature vector.
0024 featVec = reshape(repmat(d.array.frame.features',100,1),1,100*length(d.array.frame.features));
0025 flagVec = d.antenna0.receiver.flags;
0026 
0027 % We also care about the scanning rate...
0028 azRate = diff(d.antenna0.servo.az);
0029 azRate = cat(1,azRate,azRate(end));
0030 
0031 % And where exactly in azimuth are we?
0032 azVec = round(4*d.antenna0.servo.az)/4;
0033 azVec(azVec >= 360) = azVec(azVec >= 360) - 360;
0034 azVec(azVec < 0) = azVec(azVec < 0) + 360;
0035 
0036 
0037 % The data!
0038 dataVec = d.antenna0.receiver.data(:,[1 6 7 8]);
0039 
0040 
0041 %clear d;
0042 
0043 
0044 
0045 % Prepare these:
0046 dataSig = zeros(length(dataVec),4);
0047 dataNum = zeros(length(dataVec),1);
0048 
0049 foldResult = zeros(length(dataVec),1);
0050 
0051 checkVec = ones(1,length(dataVec));
0052 
0053 
0054 
0055 % Do the all-sky scans:
0056 %-----------------------------
0057 
0058 indList = featVec == 2^0 & checkVec;
0059 
0060 % Do we have any?
0061 if sum(indList) > 0
0062     
0063     % Alright.  Split scans into contiguous sections via the feature vector
0064     % and the scan speed.
0065     
0066     % Grab the indices where we have the right scan rate:
0067     medAzRate = median(abs(azRate(indList)));
0068     rateList = (abs(azRate) > (medAzRate-azRateErr)) & ...
0069         (abs(azRate) < (medAzRate+azRateErr));
0070 
0071     rateList = rateList & indList';
0072     
0073 
0074     % Deglitch this list:
0075     rateList(2:end) = rateList(2:end) & rateList(1:end-1);
0076     rateList(1:end-1) = rateList(1:end-1) & rateList(2:end);
0077 
0078     rateList(21:end) = rateList(21:end) & rateList(1:end-20);
0079     rateList(1:end-20) = rateList(1:end-20) & rateList(21:end);
0080         
0081     rateList(2:end) = rateList(2:end) | rateList(1:end-1);
0082     rateList(1:end-1) = rateList(1:end-1) | rateList(2:end);
0083     
0084     for k=1:2
0085         rateList(21:end) = rateList(21:end) | rateList(1:end-20);
0086         rateList(1:end-20) = rateList(1:end-20) | rateList(21:end);
0087     end
0088     
0089     rateList(2:end) = rateList(2:end) | rateList(1:end-1);
0090     rateList(1:end-1) = rateList(1:end-1) | rateList(2:end);
0091     
0092 
0093     checkVec = checkVec & ~rateList';
0094 
0095 
0096     % High pass filter the data...
0097     dataHP = zeros(size(dataVec));
0098     dataHP(:,1) = dataVec(:,1) - smooth(dataVec(:,1),15);
0099     dataHP(:,2) = dataVec(:,2) - smooth(dataVec(:,2),15);
0100     dataHP(:,3) = dataVec(:,3) - smooth(dataVec(:,3),15);
0101     dataHP(:,4) = dataVec(:,4) - smooth(dataVec(:,4),15);
0102 
0103     dataVec(rateList,:) = dataHP(rateList,:);
0104 
0105     
0106     
0107     % Subtract the medians from the relevant regions.
0108     dataVec = median_subtract(dataVec, rateList);
0109 
0110             
0111     % Fold on azimuth.
0112     [pSig, pNum, ~, ~] = outlier_sort(dataVec, azVec, 0.25, rateList, 3);
0113 
0114 
0115     % And bring the results back into the overall vectors:
0116     dataSig = dataSig + pSig;
0117     dataNum = dataNum + pNum;
0118 
0119 end
0120 
0121 
0122 
0123 
0124 
0125 
0126 
0127 % Do the NCP scans:
0128 %-----------------------------
0129 
0130 indList = featVec == 2^16 & checkVec;
0131 %checkVec = checkVec & ~indList;
0132 %indList = featVec == 2^16;
0133 
0134 
0135 
0136 % Do we have any?
0137 if sum(indList) > 0
0138     
0139     % Alright.  Split scans into contiguous sections via the feature vector
0140     % and the scan speed.
0141     
0142     % Grab the indices where we have the right scan rate:
0143     medAzRate = median(abs(azRate(indList)));
0144     rateList = (abs(azRate) > (medAzRate-azRateErr)) & ...
0145         (abs(azRate) < (medAzRate+azRateErr));
0146 
0147     rateList = rateList & indList';
0148     
0149 
0150     % Deglitch this list:
0151     rateList(2:end) = rateList(2:end) & rateList(1:end-1);
0152     rateList(1:end-1) = rateList(1:end-1) & rateList(2:end);
0153 
0154     rateList(21:end) = rateList(21:end) & rateList(1:end-20);
0155     rateList(1:end-20) = rateList(1:end-20) & rateList(21:end);
0156         
0157     rateList(2:end) = rateList(2:end) | rateList(1:end-1);
0158     rateList(1:end-1) = rateList(1:end-1) | rateList(2:end);
0159     
0160     for k=1:2
0161         rateList(21:end) = rateList(21:end) | rateList(1:end-20);
0162         rateList(1:end-20) = rateList(1:end-20) | rateList(21:end);
0163     end
0164     
0165     rateList(2:end) = rateList(2:end) | rateList(1:end-1);
0166     rateList(1:end-1) = rateList(1:end-1) | rateList(2:end);
0167     
0168 
0169     checkVec = checkVec & ~rateList';
0170 
0171     
0172     % High pass filter the data...
0173     dataHP = zeros(size(dataVec));
0174     dataHP(:,1) = dataVec(:,1) - smooth(dataVec(:,1),60);
0175     dataHP(:,2) = dataVec(:,2) - smooth(dataVec(:,2),60);
0176     dataHP(:,3) = dataVec(:,3) - smooth(dataVec(:,3),60);
0177     dataHP(:,4) = dataVec(:,4) - smooth(dataVec(:,4),60);
0178 
0179     dataVec(rateList,:) = dataHP(rateList,:);
0180     
0181 
0182     % Subtract the medians from the relevant regions.
0183     dataVec = median_subtract(dataVec, rateList);
0184 
0185             
0186     % Fold on azimuth.
0187     [pSig, pNum, ~, ~] = outlier_sort(dataVec, azVec, 0.25, rateList, 3);
0188 
0189 
0190     % And bring the results back into the overall vectors:
0191     dataSig = dataSig + pSig;
0192     dataNum = dataNum + pNum;
0193     
0194 end
0195 
0196 
0197 
0198 % Noise events:
0199 %-----------------------------
0200 
0201 % On and off times...
0202 onList = bitand(flagVec,4)' > 0;
0203 offList = (featVec == 2^11) & ~onList;
0204 
0205 % Get rid of transition regions:
0206 transLen = 13;
0207 onList(1+transLen:end) = onList(1:end-transLen) & onList(1+transLen:end);
0208 onList(1:end-transLen) = onList(1:end-transLen) & onList(1+transLen:end);
0209 offList(1+transLen:end) = offList(1:end-transLen) & offList(1+transLen:end);
0210 offList(1:end-transLen) = offList(1:end-transLen) & offList(1+transLen:end);
0211 % offList(2:end) = offList(1:end-1) & offList(2:end);
0212 % offList(1:end-1) = offList(1:end-1) & offList(2:end);
0213 
0214 
0215 % figure;
0216 % plot(onList);
0217 % hold on;
0218 % plot(offList,'r');
0219 % plot(checkVec,'g');
0220 % hold off;
0221 
0222 % Treat them equally, as we're going to median subtract, anyhow.
0223 indList = (onList | offList) & checkVec;
0224 % figure;
0225 % plot(indList);
0226 % input('adfadsfa');
0227 %checkVec = checkVec & ~indList;
0228 
0229 
0230 
0231 % Do we have any?
0232 if sum(indList) > 0
0233 
0234     % Alright, let's try to avoid transients near the beginning and end:
0235     rateList = smooth(indList',25) > 0.97;
0236 
0237 
0238     checkVec = checkVec & ~rateList';
0239 
0240     
0241     % Subtract the medians from the relevant regions.
0242     dataVec = median_subtract(dataVec, rateList);
0243 
0244 %     figure;
0245 %     plot(dataVec(rateList,:));
0246 %     input('adfasf');
0247 
0248     % Fold on azimuth.
0249     [pSig, pNum, ~, ~] = ...
0250         outlier_sort(dataVec, ones(size(dataVec,1),1), 1, rateList, 3);
0251 
0252 
0253     % And bring the results back into the overall vectors:
0254     dataSig = dataSig + pSig;
0255     dataNum = dataNum + pNum;
0256     
0257     
0258 %     figure;
0259 %     plot(dataSig(rateList,:));
0260 %     input('afsdfas');
0261     
0262 end
0263 
0264 
0265 
0266 
0267 % Now convert to Final result:
0268 %----------------------------------
0269 dataSigMax = max(dataSig,[],2);
0270 %foldResult = zeros(size(dataSigMax));
0271 foldResult(dataSigMax == 0) = -1; 
0272 foldResult(dataSigMax < 3) = 1; 
0273 foldResult(dataSigMax >= 3 & dataSigMax < 4) = 2; 
0274 foldResult(dataSigMax >= 4 & dataSigMax < 5) = 3; 
0275 foldResult(dataSigMax >= 5) = 4; 
0276 
0277 % figure;
0278 % plot(dataSigMax(rateList));
0279 % hold on;
0280 % plot(foldResult(rateList),'r');
0281 %
0282 % input('asdfasdfas');
0283 
0284 
0285 % Make sure that this is not an isolated spike or a short stretch of
0286 % clean data between RFI.
0287 foldResultSmooth = smooth(foldResult > 1,11);
0288 foldResult(foldResultSmooth > 0 & foldResultSmooth < 5/11) = -1;
0289 
0290 % figure;
0291 % plot(foldResult(rateList));
0292 % input('asdfasdfasdfasdfa');
0293 % figure;
0294 % plot(foldResult);
0295 % hold on;
0296 % plot(rateList-2,'r');
0297 % input('asdfasdfasdfasdfa');
0298 
0299     
0300 % Get rid of data without enough repeats...
0301 foldResult(dataNum < 10) = -1;
0302     
0303         
0304 
0305 end
0306 
0307 
0308 
0309 
0310 % Sort over a given parameter and look for outliers.
0311 function [pSig, pNum, pMean, pStd] = ...
0312     outlier_sort(dataArr, sortVec, sortStep, pList, nClip)
0313 
0314 
0315 % Prep these outputs...
0316 pNum = zeros(length(pList),1);
0317 pMean = zeros(length(pList),size(dataArr,2));
0318 pStd = zeros(length(pList),size(dataArr,2));
0319 pSig = zeros(length(pList),size(dataArr,2));
0320 
0321 
0322 % Loop over the sorting parameter.
0323 for k=min(sortVec(pList)):sortStep:max(sortVec(pList))
0324     sortInd = pList & sortVec == k;
0325 
0326 
0327     % Do we have any, here?  If not, move on.
0328     if sum(sortInd) < 3
0329         continue
0330     end
0331 
0332 
0333     % Take the mean and stddev at this sort value:
0334     sortMean = median(dataArr(sortInd,:));
0335     sortStd = std(dataArr(sortInd,:));
0336 
0337 %    nClip = 1;
0338     
0339     % Refine the mean and stddev by clipping outliers.
0340     clipInd = sortInd;
0341     for l=1:nClip
0342          clipTry = sum( ...
0343              abs(dataArr(sortInd,:) - ...
0344              repmat(sortMean,sum(sortInd),1)) ./ ...
0345              repmat(sortStd,sum(sortInd),1) > 3, 2) > 0;
0346 %        clipTry = sum( ...
0347 %            abs(dataArr(sortInd,:) - ...
0348 %            repmat(sortMean,sum(sortInd),1)) ./ ...
0349 %            repmat(sortStd,sum(sortInd),1) > 0.9, 2) > 0;
0350 
0351 
0352         clipInd(sortInd) = clipInd(sortInd) & ~clipTry;
0353 
0354         sortMean = median(dataArr(clipInd,:));
0355         sortStd = std(dataArr(clipInd,:));
0356 
0357     end
0358 
0359 
0360 
0361     % Populate the output vectors:
0362     pNum(sortInd) = sum(sortInd);
0363 
0364 
0365     % And the output vector of real interest:
0366     pSig(sortInd,:) = ...
0367         abs(dataArr(sortInd,:) - ...
0368         repmat(sortMean,sum(sortInd),1)) ./ ...
0369         repmat(sortStd,sum(sortInd),1);
0370 
0371 
0372 end
0373 
0374 
0375 
0376 end
0377 
0378 
0379 
0380 function dataArr = median_subtract(dataArr, sortList)
0381 
0382 
0383 % OK, so when do the scans start and stop?
0384 startVec = find(diff(sortList) > 0);
0385 stopVec = find(diff(sortList) < 0);
0386 
0387 
0388 % Check for anomalies.
0389 if (length(stopVec) < 2) || (length(startVec) < 2)
0390     disp('PROBLEM!?!?  Not enough start/stop times.');
0391     return
0392 end
0393 
0394 if stopVec(1) < startVec(1)
0395     startVec = cat(1,1,startVec);
0396 end
0397 if startVec(end) > stopVec(end)
0398     stopVec = cat(1,stopVec,length(sortList));
0399 end
0400 
0401 if length(stopVec) ~= length(startVec)
0402     disp('PROBLEM!?!?  Start and stop vectors are not of same length.');
0403     return
0404 end
0405 
0406 
0407 % Now we do the median subtractions...
0408 for k=1:length(startVec)
0409 
0410     startInd = startVec(k);
0411     endInd = stopVec(k);
0412 
0413     dataArr(startInd:endInd,:) = dataArr(startInd:endInd,:) - ...
0414         repmat(median(dataArr(startInd:endInd,:)),endInd-startInd+1,1);
0415 end
0416 
0417 
0418 
0419 end