plotTemplateEstimates

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

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

   XXX


   I/O:

   MAS -- 20-April-2012

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

0001 function [mjdList changeList] = ...
0002     plotTemplateEstimates(mjdVec, estimateArray, mjdList)
0003 
0004 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0005 %
0006 % function templateArray = plotTemplateEstimates(utcStart,utcEnd,archiveDir,nTheta)
0007 %
0008 %   XXX
0009 %
0010 %
0011 %   I/O:
0012 %
0013 %   MAS -- 20-April-2012
0014 %
0015 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0016 
0017 %disp(mjdList);
0018 
0019 
0020 % figure;
0021 % plot(squeeze(estimateArray(:,1,:)));
0022 % disp(estimateArray(:,1,1));
0023 % input('asdfasdfa');
0024 
0025 
0026 % These are hardcoded.  Deal with it.
0027 polNames = {'I1'; 'Q3'; 'U3'; 'I2'};
0028 
0029 
0030 
0031 % Need to zero-pad for plotting purposes.
0032 n_phase = size(estimateArray,1);
0033 n_pol = length(polNames);
0034 n_pad = round(12*(mjdVec(end) - mjdVec(1)))+1;
0035 
0036 if n_pad > length(mjdVec)
0037     % Easy enough.
0038     mjdVec_pad = mjdVec(1) + (0:n_pad-1)/12;
0039 
0040     % Initialize this guy.
0041     estimateArray_pad = zeros(n_phase,n_pad,n_pol);
0042 
0043     j = 0;
0044     dt = round(diff(12*mjdVec));
0045     for k=1:length(mjdVec)
0046 %        disp(dt(k));
0047         if k > 1 && dt(k-1) > 1
0048             % Find the skippers.
0049             j = j + dt(k-1) - 1;
0050         end
0051         estimateArray_pad(:,k+j,:) = estimateArray(:,k,:);
0052     end
0053 else
0054     % No padding needed.
0055     mjdVec_pad = mjdVec;
0056     estimateArray_pad = estimateArray;
0057 end
0058 
0059 
0060 
0061 % Welcome the user to this most wonderful of experiences.
0062 displayintro();
0063 
0064 
0065 % Y tick labels, yo.
0066 thetaVec = ((1:n_phase) - 0.5) / n_phase * 2 * pi;
0067 
0068 
0069 % Quickly rescale the colour scale.
0070 % estimateMax = mean(estimateArray(1:round(n_phase/2/pi),:,:));
0071 % estimateMin = mean(estimateArray(round(3*n_phase/2/pi):round(4*n_phase/2/pi),:,:));
0072 % scaleEstimate = abs(estimateMax - estimateMin);
0073 
0074 scaleEstimate = 2.0 * median(abs(estimateArray_pad));
0075 
0076 estimateArray_pad = estimateArray_pad ./ repmat(scaleEstimate, [n_phase 1 1]);
0077 
0078 % Preparing the outputs.
0079 %mjdList = mjdList;
0080 changeList = zeros(size(mjdList));
0081 
0082 
0083 % These variables are used by the while loop to determine behaviour.
0084 % Setting them up with default values so that the first loop iteration will
0085 % perform as needed.
0086 whichpol = 1;
0087 whichpol_last = 0;
0088 zoomrange = [min(mjdVec_pad) max(mjdVec_pad)];
0089 deletebound = 0;
0090 createbound = 0;
0091 nbound = length(mjdList);
0092 
0093 
0094 % Here's the loop.  WHILE style, baby.
0095 loopval = 1;
0096 while loopval
0097     
0098     % Grabs the current figure.  Makes a new one if needed.
0099     gcf;
0100     axis ij;
0101     hold on;
0102     ylabel('Phase (radians)');
0103     xlabel('MJD');
0104 
0105     
0106     % Depending upon the values in the control variables, a given iteration
0107     % of the loop will do different things.
0108     
0109     % First, if needed, we re-plot a different polarization.
0110     if whichpol ~= whichpol_last
0111 
0112         % Let the user know which polarization s/he's looking at.
0113         title(polNames{whichpol});
0114         
0115 % %         % Calculate the color scaling limit.
0116 %         estimatePol = estimateArray(:,:,whichpol);
0117 % %         climit = std(estimatePol(:));
0118 % %         figure;
0119 % %         figure;
0120 %         plot(squeeze(estimateArray(:,1,2)));
0121 %         figure;
0122 %         plot(estimatePol(:));
0123 %         input('asdfasdfa');
0124         
0125         
0126         % Plot a colourmap of the templates.
0127         imagesc(mjdVec_pad,thetaVec,estimateArray_pad(:,:,whichpol), [-1.0 1.0]);
0128 
0129         % The template boundaries are shown as vertical, magenta, dashed
0130         % lines.
0131         handleVec = zeros(nbound,1);
0132         for k=1:nbound
0133             handleVec(k) = ...
0134                 plot([mjdList(k) mjdList(k)], thetaVec([1 end]), 'm--', 'LineWidth', 4);
0135         end
0136         
0137         % Keeping the control variables in order.
0138         whichpol_last = whichpol;
0139         
0140         
0141     % In this case, the user has decided to create a new load boundary.
0142     elseif createbound ~= 0
0143 
0144         % Round off the boundary to nearest 2-hour.
0145         createbound = round(createbound * 12) / 12;
0146         
0147         % Keep the user apprised of the situation.
0148         disp(['Creating new boundary at ' mjd2string(createbound)]);
0149         
0150         % Add the new boundary to the list.
0151         mjdList(end+1) = createbound;
0152         changeList(end+1) = 1;
0153         % Plot the new boundary.
0154         handleVec(end+1) = ...
0155             plot([createbound createbound], thetaVec([1 end]),'m--','LineWidth',4);
0156 
0157         % Keep track of the number of boundaries.
0158         nbound = nbound + 1;
0159         
0160         % Keeping the control variables in order.
0161         createbound = 0;
0162 
0163 
0164     % In this case, the user has decided to delete an existing load
0165     % boundary.
0166     elseif deletebound ~= 0 && nbound > 0
0167         
0168         % Get the closest boundary to the mouse click.
0169         [poo, delInd] = min(abs(mjdList - deletebound));
0170 
0171         if min(abs(mjdList(delInd)-mjdVec(1)),abs(mjdList(delInd)-mjdVec(end))) <= 1/24 
0172             disp('Sorry.  Cannot delete boundary so close to plot edge.')
0173 
0174         else
0175             % Verify that the user actually wants to delete this boundary.
0176             delQuery = ...
0177                 input(['Delete boundary at ' mjd2string(mjdList(delInd)) ...
0178                 '? (Y/N)  '],'s');        
0179 
0180             % Delete it.
0181             if strcmpi(delQuery,'y')
0182                 delete(handleVec(delInd));
0183                 handleVec(delInd) = [];
0184                 mjdList(delInd) = [];
0185 
0186                 % This is part of keeping track of changes in template length.
0187                 if changeList(delInd) == 0 && delInd > 1
0188                     changeList(delInd-1) = 2;
0189                 end
0190                 changeList(delInd) = [];
0191 
0192                 nbound = nbound - 1;
0193             end
0194         end
0195         % Keeping the control variables in order.
0196         deletebound = 0;
0197     end
0198 
0199     
0200     % Set the axis limits.
0201     ylim([0 2*pi]);
0202     xlim(zoomrange);
0203 
0204     
0205     hold off;
0206     
0207     
0208     % Now wait for input.
0209     [x, poo, button] = ginput(1);
0210     
0211     % If we don't check for this, then we can crash during the switch
0212     % statement.
0213     if isempty(button)
0214         continue
0215     end
0216     
0217     
0218     % What to do with the input, eh?
0219     switch button
0220 
0221         % q: quit
0222         case 113
0223             fprintf('\nYou have chosen to quit the interactive stage.\n');
0224             quitQuery = input( ...
0225                 'Are you sure? (Y/N)  ','s');
0226             if strcmpi(quitQuery,'y')
0227                 loopval = 0;
0228             end
0229             
0230         % z: zoom
0231         case 122
0232             disp('Zooming...');
0233             % Get the second click.
0234             [x2, poo, button] = ginput(1);
0235             
0236             if abs(x2-x) < 1/20. 
0237                 disp('  ... zero zoom range.  Try again.');
0238             % Set the control variable.  The zooming will occur during the
0239             % next loop iteration.
0240             elseif button == 122 || button == 1
0241                 zoomrange = sort([x x2]);
0242             end
0243 
0244         % u: unzoom
0245         case 117
0246             % Set the control variable.  The unzooming will occur during
0247             % the next loop iteration.
0248             zoomrange = [min(mjdVec_pad) max(mjdVec_pad)];
0249 
0250         % d: delete boundary
0251         case 100
0252             % Record the mouse position.  Boundary will be deleted during
0253             % the next loop iteration.
0254             deletebound = x;
0255             
0256         % c: create boundary
0257         case 99
0258             % Record the mouse position.  Boundary will be created during
0259             % the next loop iteration.
0260             createbound = x;
0261             
0262         % ]: next polarization
0263         case 93
0264             % Set the control variable.  Polarization will be changed
0265             % during the next loop iteration.
0266             if whichpol < 4
0267                 whichpol = whichpol + 1;
0268             else
0269                 whichpol = 1;
0270             end
0271             
0272         % [: previous polarization
0273         case 91
0274             % Set the control variable.  Polarization will be changed
0275             % during the next loop iteration.
0276             if whichpol > 1
0277                 whichpol = whichpol - 1;
0278             else
0279                 whichpol = 4;
0280             end
0281             
0282         % p: select polarization
0283         case 112
0284             % Query the user for his/her preferred polarization.
0285             whichpolS = ...
0286                 input('Which polarization? (1=I1, 2=Q3, 3=U3, 4=I2)  ','s');
0287 
0288             % Set the control variable.  Polarization will be changed
0289             % during the next loop iteration.
0290             switch whichpolS
0291                 case '1'
0292                     whichpol = 1;
0293                 case '2'
0294                     whichpol = 2;
0295                 case '3'
0296                     whichpol = 3;
0297                 case '4'
0298                     whichpol = 4;
0299                 otherwise
0300                     disp('Invalid choice.');
0301             end
0302             
0303         % t: give utc time of mouse position
0304         case 116
0305             % Tell the user where the mouse is pointing.
0306             disp(['Mouse position: ' mjd2string(x)]);
0307             
0308         % h: help
0309         case 104
0310             % Display the help message.
0311             displayhelp();
0312     end
0313 end
0314 
0315 
0316 % Sort the boundaries.
0317 % disp(size(mjdList));
0318 % disp(size(changeList));
0319 [mjdList,sortI,sortJ] = unique(mjdList,'last');
0320 %[mjdList sortI] = sort(mjdList);
0321 changeList = changeList(sortI);
0322 
0323 % Note those load templates with changing length.
0324 lenChange = find(changeList(1:end-1) == 0 & changeList(2:end) == 1);
0325 if ~isempty(lenChange)
0326     changeList(lenChange) = 2;
0327 end
0328 
0329 % if changeList(end) == 0
0330 %     changeList(end) = 2;
0331 % end
0332 changeList(end) = 2 - changeList(end);
0333 
0334 close(gcf);
0335 
0336 end
0337 
0338 
0339 
0340 function r = displayintro()
0341 
0342 disp('====================================');
0343 disp('= Welcome to the Interactive       =');
0344 disp('= Load Boundary Identification     =');
0345 disp('====================================');
0346 disp(' ');
0347 disp('You are cordially requested to view the');
0348 disp('2-hour load templates for the filtered');
0349 disp('I1, Q3, U3, and I2 channels.  Previously');
0350 disp('identified boundaries are denoted by');
0351 disp('magenta lines.  Please use your judgement ');
0352 disp('to add or remove cold load boundaries.');
0353 disp(' ');
0354 disp('Fewer boundaries means higher SNR for');
0355 disp('the templates.  Boundaries are necessary, ');
0356 disp('however, when changes in the ');
0357 disp('template structure occur.');
0358 disp(' ');
0359 disp('Press h for help.');
0360 disp(' ');
0361 disp('Take care.');
0362 disp(' ');
0363 disp('====================================');
0364 disp(' ');
0365 
0366 r = 1;
0367 end
0368 
0369 
0370 
0371 function r = displayhelp()
0372 
0373 disp('====================================');
0374 disp('= Interactive Load Boundary Help   =');
0375 disp('====================================');
0376 disp(' ');
0377 disp('Position your mouse pointer over the ');
0378 disp('plot window and use keyboard to input');
0379 disp('commands');
0380 disp(' ');
0381 disp('c -- Create a new load boundary at the');
0382 disp('       chosen location.');
0383 disp('d -- Delete the load boundary closest ');
0384 disp('       to the cursor.');
0385 disp('z -- Zoom in mjd.  Press at the ');
0386 disp('       desired extrema.');
0387 disp('u -- Unzoom to full length.');
0388 disp('[ -- View the previous polarization.');
0389 disp('] -- View the next polarization.');
0390 disp('p -- Choose a polariation to view.');
0391 disp('t -- UTC at the current mouse position.');
0392 disp('h -- This help message.');
0393 disp('q -- Quit the interactive step and move on.');
0394 disp(' ');
0395 disp('====================================');
0396 disp(' ');
0397 
0398 r = 1;
0399 end