compileLoadTemplates

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

 function newTemplates = plotTemplateEstimates(utcStart,utcEnd,archiveDir,nTheta)

   XXX


   I/O:

   MAS -- 24-April-2012

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

0001 function newTemplates = compileLoadTemplates(mjdVec, numArray, estimateArray, mjdList, recalcList, oldMjdList, oldTemplates)
0002 
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %
0005 % function newTemplates = plotTemplateEstimates(utcStart,utcEnd,archiveDir,nTheta)
0006 %
0007 %   XXX
0008 %
0009 %
0010 %   I/O:
0011 %
0012 %   MAS -- 24-April-2012
0013 %
0014 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0015 
0016 
0017 nTemplates = length(mjdList);
0018 nPhase = size(estimateArray,1);
0019 nPols = size(estimateArray,3);
0020 
0021 thetaVec = ((1:nPhase) - 0.5) / nPhase * 2 * pi;
0022 
0023 newTemplates = zeros(nPhase,nTemplates,nPols);
0024 
0025 % Estimates will be filtered.
0026 [home,installeddir] = where_am_i();
0027 f = open([home '/' installeddir '/reduc/load/filters/template_filter_1024.mat']);
0028 fLen = length(f.Num);
0029 
0030 
0031 % We'll ignore the data between thetas of 1 and 3 radians for flagging and
0032 % scaling purposes.
0033 phaseFlag = zeros(nPhase,1);
0034 phaseFlag(round(nPhase/2/pi):round(3*nPhase/2/pi)) = 1;
0035 
0036 
0037 
0038 % Loop over the templates.
0039 
0040 for k=1:nTemplates
0041 
0042     if recalcList(k) == 0
0043         % Not changed.  Keep the old template for this MJD.
0044 
0045         [poo, oldIndex] = min(abs(oldMjdList - mjdList(k)));
0046 
0047         newTemplates(:,k,:) = oldTemplates(:,oldIndex,:);
0048         
0049         
0050     else
0051         % Changed!  Calculate the new one.
0052 
0053         % Identify the start and stop times for this template.
0054         startMJD = mjdList(k);
0055         if k < length(mjdList)
0056             endMJD = mjdList(k+1) - 1 / 12;
0057         else
0058             endMJD = max(mjdVec);
0059         end
0060 
0061         fprintf(['\nNow calculating template beginning \n' ...
0062             mjd2string(startMJD) ...
0063             '\nand ending ' mjd2string(endMJD) '\n']);
0064 
0065         % Identify the template estimates that will be used.
0066         [poo, startIndex] = nanmin(abs(mjdVec - startMJD));
0067         [poo, endIndex] = nanmin(abs(mjdVec - endMJD));
0068 
0069         % Grab the estimates of interest.
0070         estimateSubArray = estimateArray(:,startIndex:endIndex,:);
0071         numSubArray = numArray(:,startIndex:endIndex);
0072         
0073         
0074         % Filter the estimates.
0075         estimateFilter = filter(f.Num,1,repmat(estimateSubArray, [3 1 1]));
0076         estimateSubArray = estimateFilter(ceil(fLen/2)+nPhase:ceil(fLen/2)+2*nPhase-1,:,:);
0077         numVec = sum(numSubArray);
0078         
0079         
0080         % Do the median/stddev rejection.
0081         estimateSubSub = estimateSubArray(phaseFlag == 0,:,[1 6 7 8]);
0082         flagVec = template_reject(estimateSubSub);
0083 %         disp(flagVec);
0084         
0085 
0086         % Calculate the mean templates (ignoring flagged estimates).
0087         templateMean = template_mean(estimateSubArray,flagVec,numVec);
0088 
0089 
0090         % Do the mean/stddev rejection.
0091         estimateSubSub = estimateSubArray(phaseFlag == 0,:,[1 6 7 8]);
0092         templateMeanSub = templateMean(phaseFlag == 0,:,[1 6 7 8]);
0093         flagVec = template_reject(estimateSubSub, templateMeanSub);
0094 %         disp(sum(isnan(estimateSubSub(:))));
0095 %         disp(sum(isnan(templateMeanSub(:))));
0096 %         disp(sum(isnan(flagVec(:))));
0097 %         disp(flagVec);
0098 
0099 
0100         % Calculate the mean templates (again!).
0101         templateMean = template_mean(estimateSubArray,flagVec,numVec);
0102 %         disp(sum(isnan(templateMean(:))));
0103 
0104 
0105         % Loop back through, calculating scaling factor and DC offset for each
0106         % estimate.
0107         estimateSubSub = estimateSubArray(phaseFlag == 0,:,:);
0108         templateMeanSub = templateMean(phaseFlag == 0,:,:);
0109         [mArray bArray] = template_fit(estimateSubSub, templateMeanSub);
0110 %         disp(sum(isnan(estimateSubSub(:))));
0111 %         disp(sum(isnan(templateMeanSub(:))));
0112 
0113 
0114         % Scale!
0115         estimateSubScale = estimateSubArray ./ ...
0116             repmat(mArray, [nPhase 1 1]) - repmat(bArray, [nPhase 1 1]);
0117 %         disp(sum(isnan(estimateSubScale(:))));
0118         
0119         
0120         % Do the mean/stddev rejection once more, but taking scaling into account.
0121         estimateSubSub = estimateSubScale(phaseFlag == 0,:,[1 6 7 8]);
0122         templateMeanSub = templateMean(phaseFlag == 0,:,[1 6 7 8]);
0123         flagVec = template_reject(estimateSubSub, templateMeanSub);
0124 %         disp(sum(isnan(estimateSubSub(:))));
0125 %         disp(sum(isnan(flagVec(:))));
0126 %         disp(sum(isnan(numVec(:))));
0127 %         disp(flagVec);
0128 
0129 
0130         % Recalculate the mean templates.
0131         if sum(flagVec == 0) >= 2
0132             % Calculate the templates...
0133             newTemplates(:,k,:) = template_mean(estimateSubScale,flagVec,numVec);        
0134         else
0135             % If everything is flagged, then... get rid of the flagging.
0136             newTemplates(:,k,:) = template_mean(estimateSubScale,flagVec*0,numVec);        
0137         end
0138 
0139         
0140         
0141         % Plot the latest templates, for the user's interest.
0142         gcf;
0143 
0144         i1 = plot(thetaVec,newTemplates(:,k,1),'color','b');
0145         hold on;
0146         q3 = plot(thetaVec,newTemplates(:,k,6),'color','g');
0147         u3 = plot(thetaVec,newTemplates(:,k,7),'color','m');
0148         i2 = plot(thetaVec,newTemplates(:,k,8),'color','r');
0149 
0150         xlabel('Phase (radians)');
0151         ylabel('Amplitude');
0152 
0153         legend([i1 q3 u3 i2],'I1','Q3','U3','I2');
0154         hold off;
0155         input('Sample templates for filtered data.  Press enter.');
0156         
0157     end
0158 end
0159 
0160 end
0161 
0162 
0163 % Used for rejecting outliers from the template estimates.  Can optionally
0164 % provide a templateModel.  If not, then the median of the estimates is
0165 % used.
0166 function flagVec = template_reject(estimateArray, templateModel)
0167 
0168 nPhase = size(estimateArray,1);
0169 
0170 
0171 if nargin < 2
0172     % MEDIAN-BASED REJECTION
0173     templateModel = nanmedian(estimateArray,2);
0174 end
0175 
0176 
0177 % Offer a quick normalization.
0178 templateModel = ...
0179     (templateModel - repmat(nanmean(templateModel),[nPhase 1 1])) ./ ...
0180     repmat(nanstd(templateModel),[nPhase 1 1]);
0181 estimateArray = ...
0182     (estimateArray - repmat(nanmean(estimateArray),[nPhase 1 1])) ./ ...
0183     repmat(nanstd(estimateArray),[nPhase 1 1]);
0184 
0185 
0186 flagVec = zeros(1,size(estimateArray,2));
0187 
0188 
0189 for k=1:5
0190     templateStd = nanstd(estimateArray(:,~flagVec,:),[],2);
0191 
0192 %     figure;
0193 %     plot(templateModel(:,1));
0194 %     hold on;
0195 %     plot(templateStd(:,1),'r');
0196 %     input('asdfasfd');
0197     
0198     
0199     flagVec = ...
0200         nanmax(nanmax( ...
0201         abs(estimateArray - ...
0202         repmat(templateModel,[1 size(estimateArray,2) 1])) > ...
0203         5*repmat(templateStd,[1 size(estimateArray,2) 1]),[],3),[],1) ...
0204         | ...
0205         max(max(isnan(estimateArray),[],3),[],1);
0206 %     figure;
0207 %     plot(templateModel(:,1));
0208 %     hold on;
0209 %     plot(estimateArray(:,flagVec,1));
0210 %     figure;
0211 %     plot(flagVec,'r');
0212 %     input('adfasd');
0213 end
0214 
0215 end
0216 
0217 
0218 % Calculate the weighted mean of the template estimates.  Just one line,
0219 % but it's a long line.
0220 function templateMean = template_mean(estimateArray,flagVec,numVec)
0221 
0222 nPhase = size(estimateArray,1);
0223 
0224 % Weighted by numArray.  It's as simple as this, yo.
0225 templateMean = ...
0226     nansum(estimateArray(:,~flagVec,:) .* ...
0227     repmat(numVec(~flagVec),[nPhase 1 20]),2) ./ ...
0228     repmat(nansum(numVec(~flagVec),2),[nPhase 1 20]);
0229 % figure;
0230 % plot(squeeze(templateMean));
0231 
0232 end
0233 
0234 
0235 % Fit the multiplicative and additive offsets to a given template estimate.
0236 function [mArray bArray] = template_fit(estimateArray, templateMean)
0237 
0238 % figure;
0239 % plot(estimateArray(:,:,1));
0240 % disp(estimateArray(:,20,1));
0241 % input('asdfsa');
0242 
0243 % disp(size(templateMean));
0244 %
0245 % figure;
0246 % plot(squeeze(templateMean));
0247 % disp(templateMean(:,1));
0248 % input('385743958');
0249 
0250 
0251 nPhase = size(estimateArray,1);
0252 
0253 
0254 mArray = zeros(1,size(estimateArray,2),size(estimateArray,3));
0255 bArray = zeros(1,size(estimateArray,2),size(estimateArray,3));
0256 
0257 for k=1:size(estimateArray,3)
0258 
0259     A = [ones(nPhase,1) templateMean(:,k)];
0260 
0261     for l=1:size(estimateArray,2)
0262 
0263 %         figure;
0264 %         plot(A(:,2));
0265 %         hold on;
0266 %         plot(squeeze(estimateArray(:,l,k)),'r');
0267 %         input('adfasf');
0268         
0269         x = A \ estimateArray(:,l,k);
0270 
0271 %         disp(x);
0272 %         input('asdfas');
0273         mArray(1,l,k) = x(2);
0274         bArray(1,l,k) = x(1);
0275     end
0276 end
0277 
0278 end
0279