--- a
+++ b/Analysis/jPCA_ForDistribution/phaseSpace.m
@@ -0,0 +1,373 @@
+% For making publication quality rosette plots
+% useage:
+%       phaseSpace(Projection, Summary)
+%       phaseSpace(Projection, Summary, params)
+%
+% You can limit the conditions you plot either by
+%   1) including only some entries of Projection (e.g., Projection(1:2))
+%   or 2) by using params.conds2plot to restrict (e.g., params.conds2plot = 1:2);
+%   In the former case, scaling etc is based just on the passed points.
+%   In the latter scaling is based on all the points.
+%
+% outputs are [colorStruct, hf, haxP, vaxP] = phaseSpace(Projection, Summary)
+%   
+%   'colorStruct' is a structure (one per dataset) of cells of linecolors (one per condition)
+%   that you might wish to pass to another function (e.g., one that plots the rosette PSTH's 
+%   or the hand trajectories).
+%
+%    'hf' is the fig#, 'haxP' and 'vaxP' are the axis parameters.
+%
+%   *** ALL of the outputs, except the last, pertain only the LAST graph plotted. ***
+%
+% rosetteData comes from multiRosetteScript
+%
+% params can have the following fields:
+%       .times    These override the default times (those corresponding to scores; e.g., the orginal
+%               times that were used to build the space.  If empty, the defaults are used.
+%                 Note that zero is movement onset.  Thus .times will probably start strongly
+%               negatively.  
+%                 Thus, you might pass it -1550:10:150 if you wanted to start way back at the
+%               beginning.
+%                 A nice feature is that only those times that match times in 'scoresExtraTime' are
+%               used.  Thus, if you pass -100000:1000000 things will still work fine.
+%               Scalings are based on all times in Projection.projAllTimes.  This is nice for movies
+%               as the scaling won't change as a function of the times you plot.
+%
+%       .planes2plot   list of the jPC planes you want plotted.  Default is [1].  [1,2] would also be reasonable.
+%
+%       .arrowSize        The default is 5
+%       .arrowGain        FOR MOVIES: sets velocity dependence of arrow size (0 to not grow when faster).
+%       .plotPlanEllipse  Controls whether the ellipse is plotted.  The default is 'true'
+%       .useAxes          whether axes should be plotted.  Default is 'true'
+%       .useLabel         whether to label with dataset. Default is 'true'
+%       .planMarkerSize   size of the plan dot.  Default is 6.
+%       .lineWidth        width of the trajectories.  Default is 0.85.
+%       .arrowMinVel      minimum velocity for plotting an arrow.
+%       .rankType         default is 'eig', but you can override with 'varCapt'.  The first plane will
+%                         then be the jPC plane that captured the most variance (often associated with the 
+%                         largest eigenvalue but not always
+%       .conds2plot       which conditions to plot (scalings will still be based on all the conds in 'Projection')
+%       .substRawPCs      use PC projections rather than jPC projections
+%       .crossCondMean    if present and == 1, plot the cross condition mean in cyan.
+%       .reusePlot        if present and == 1, do cla then reuse the plot
+%       .dataRanges       normally this is set automatically, but you can decide yourself what the
+%                         range should be.  You should supply one entry per plane to be plotted.  You can also supply
+%                         just the first, and then the defaults will be used after that.
+%
+function [colorStruct, haxP, vaxP] = phaseSpace(Projection, Summary, params)
+
+
+%% some basic parameters
+axLimScale = 1.35; 
+axisSeparation = 0.20;  % separated by 20% of the maximum excursion (may need more if plotting future times, which aren't used to compute farthestLeft or farthestDown)
+
+numPlanes = length(Summary.varCaptEachPlane);  % total number of planes provided (may only plot a subset)
+
+%% set defaults and override if 'params' is included as an argument
+
+% allows for the use of times other than the original ones that correspond to 'scores' (those that
+% were used to create the projection and do the analysis)
+overrideTimes = [];
+if exist('params', 'var') && isfield(params,'times')
+    overrideTimes = params.times;
+end
+
+arrowSize = 5;
+if exist('params', 'var') && isfield(params,'arrowSize')
+    arrowSize = params.arrowSize;
+end
+
+arrowGain = 0;
+if exist('params', 'var') && isfield(params,'arrowGain')
+    arrowGain = params.arrowGain;
+end
+
+% Default is we plot the ellipse if we have 6 or more conditions
+% NOTE: we still plot it even if asked to only plot one cond, so long as we HAVE more than 6 to
+% build the ellipse off.  It matters whether length(Projection) >= 6, not whether conds2plot >= 6
+if length(Projection) >= 6
+    plotPlanEllipse = true;
+else
+    plotPlanEllipse = false; 
+    axLimScale = 1.3*axLimScale;
+end
+if exist('params', 'var') && isfield(params,'plotPlanEllipse')
+    plotPlanEllipse = params.plotPlanEllipse;
+end
+
+useAxes = true;
+if exist('params', 'var') && isfield(params,'useAxes')
+    useAxes = params.useAxes;
+end
+
+useLabel = true;
+if exist('params', 'var') && isfield(params,'useLabel')
+    useLabel = params.useLabel;
+end
+
+planMarkerSize = 6;
+if exist('params', 'var') && isfield(params,'planMarkerSize')
+    planMarkerSize = params.planMarkerSize;
+end
+
+lineWidth = 0.85;
+if exist('params', 'var') && isfield(params,'lineWidth')
+    lineWidth = params.lineWidth;
+end
+
+arrowMinVel = [];
+if exist('params', 'var') && isfield(params,'arrowMinVel')
+    arrowMinVel = params.arrowMinVel;
+end
+
+planes2plot = [1];  % this is a list of which planes to plot 
+if exist('params', 'var') && isfield(params,'planes2plot')
+    planes2plot = params.planes2plot;
+end
+
+rankType = 'eig';
+if exist('params', 'var') && isfield(params,'rankType')
+    rankType = params.rankType;
+end
+
+reusePlot = 0;
+if exist('params', 'var') && isfield(params,'reusePlot')
+    reusePlot = params.reusePlot;
+end
+if length(planes2plot) > 1, reusePlot = 0; end  % cant reuse if we are plotting more than one thing.
+
+numConds = length(Projection);
+conds2plot = 1:numConds;
+if exist('params', 'var') && isfield(params,'conds2plot')
+    if ~strcmp(params.conds2plot,'all')
+        conds2plot = params.conds2plot;
+    end
+end
+
+% if asked, plot the cross-condition mean on the same plot.
+crossCondMean = false;
+if exist('params', 'var') && isfield(params,'crossCondMean')
+    crossCondMean = params.crossCondMean;
+end
+
+
+% If asked, substitue the raw PC projections.
+substRawPCs = 0;
+if exist('params', 'var') && isfield(params,'substRawPCs') && params.substRawPCs
+    substRawPCs = 1;
+    % just overwrite
+    for c = 1:numConds
+        Projection(c).proj = Projection(c).tradPCAproj;
+        Projection(c).projAllTimes = Projection(c).tradPCAprojAllTimes;
+    end
+    Summary.varCaptEachPlane = sum(reshape(Summary.varCaptEachPC,2,numPlanes));
+end
+
+if strcmp(rankType, 'varCapt')  && substRawPCs == 0  % we WONT reorder if they were PCs
+    [~, sortIndices] = sort(Summary.varCaptEachPlane,'descend');
+    planes2plot_Orig = planes2plot;  % keep this so we can label the plane appropriately
+    planes2plot = sortIndices(planes2plot);  % get the asked for planes, but by var accounted for rather than eigenvalue
+end
+
+% the range of the data will set the size of the plot unless you manually override
+dataRanges = max(abs(vertcat(Projection.proj)));
+dataRanges = max(reshape(dataRanges,2,numPlanes));  % one range per plane
+
+if exist('params', 'var') && isfield(params,'dataRanges')
+    for i = 1:length(params.dataRanges)
+        dataRanges(i) = params.dataRanges(i);  % only override those values that are specified
+    end
+end
+    
+
+
+arrowEdgeColor = 'k';
+
+for pindex = 1:length(planes2plot)
+    
+    % get some useful indices
+    plane = planes2plot(pindex);  % which plane to plot
+    d2 = 2*plane;  % indices into the dimensions
+    d1 = d2-1;
+
+    % set the limits of the figure
+    axLim = axLimScale * dataRanges(plane) * [-1 1 -1 1];
+    axisLength = 0.5;
+
+    
+    for c = 1:numConds
+        % Always taken from the 1st element (NOT the first that will be plotted)
+        % This way the ellipse doesn't depend on which times you choose to plot 
+        planData(c,:) = Projection(c).proj(1,[d1,d2]);  
+    end
+    % need this for ellipse plotting
+    if numConds > 1   % used to be only if plotPlanEllipse==1.  Now we always use the same scaling regardless of whether we use the plan ellipse.
+        ellipseRadii = 2*var(planData).^0.5;  % we may plot an ellipse for the plan activity
+
+        % these will be altered further below based on how far the data extends left and down
+        farthestLeft = -ellipseRadii(1);  % used figure out how far the axes need to be offset
+        farthestDown = -ellipseRadii(2);  % used figure out how far the axes need to be offset
+    else 
+        temp = vertcat(Projection.proj);
+        temp = temp(:,[d1,d2]);
+        farthestLeft = -1.2*max(temp(:));  % used figure out how far the axes need to be offset
+        farthestDown = -1.2*max(temp(:));  % used figure out how far the axes need to be offset
+    end
+
+    %% deal with the color scheme
+
+    % ** colors graded based on PLAN STATE
+    % These do NOT depend on which times you choose to plot (only on which time is first in Projection.proj).
+    htmp = redgreencmap(numConds, 'interpolation', 'linear');
+    [~,newColorIndices] = sort(planData(:,1));
+
+    htmp(newColorIndices,:) = htmp;
+
+    for c = 1:numConds  % cycle through conditions, and assign that condition's color
+        lineColor{c} = htmp(c,:);
+        arrowFaceColor{c} = htmp(c,:);
+        planMarkerColor{c} = htmp(c,:);
+    end
+
+    % override colors if asked
+    if exist('params', 'var') && isfield(params,'colors')
+        lineColor = params.colors;
+        arrowFaceColor = params.colors;
+        planMarkerColor = params.colors;
+        disp('hi');
+    end
+    
+    colorStruct(pindex).colors = lineColor;
+
+    %% Plot the rosette itself
+    if reusePlot == 0, blankFigure(axLim); else cla; end
+
+    % first deal with the ellipse for the plan variance (we want this under the rest of the data)
+    if plotPlanEllipse, circle([0 0], ellipseRadii, 0.6*[1 1 1], 1); end
+
+    % cycle through conditions
+    for c = 1:numConds
+
+        if isempty(overrideTimes)  % if we are going with the original times (those that were used to create the projection and do the analysis)        
+            P1 = Projection(c).proj(:,d1);
+            P2 = Projection(c).proj(:,d2);       
+        else
+            useTimes = ismember(Projection(c).allTimes, overrideTimes);
+            P1 = Projection(c).projAllTimes(useTimes,d1);
+            P2 = Projection(c).projAllTimes(useTimes,d2);
+        end
+
+        if ismember(c,conds2plot)
+            plot(P1, P2, 'color', lineColor{c}, 'lineWidth', lineWidth);
+
+            if planMarkerSize>0
+                plot(P1(1), P2(1), 'ko', 'markerSize', planMarkerSize, 'markerFaceColor', planMarkerColor{c});
+            end
+
+            % for arrow, figure out last two points, and (if asked) supress the arrow if velocity is
+            % below a threshold.
+            penultimatePoint = [P1(end-1), P2(end-1)];
+            lastPoint = [P1(end), P2(end)];
+            vel = norm(lastPoint - penultimatePoint);
+            if isempty(arrowMinVel) || vel > arrowMinVel
+                aSize = arrowSize + arrowGain * vel;  % if asked (e.g. for movies) arrow size may grow with vel
+                arrowMMC(penultimatePoint, lastPoint, [], aSize, axLim, arrowFaceColor{c}, arrowEdgeColor);
+            else
+                plot(lastPoint(1), lastPoint(2), 'ko', 'markerSi', arrowSize, 'markerFac', arrowFaceColor{c}, 'markerEdge', arrowEdgeColor);
+            end
+        end
+
+        % axis locations will be based on the original set of times used to make the scores
+        % and not on the actual times used.  Here we get the leftmost and bottommost point
+        if isfield(Projection, 'projAllTimes')
+            farthestLeft = min(farthestLeft, min(Projection(c).projAllTimes(:,d1)));
+            farthestDown = min(farthestDown, min(Projection(c).projAllTimes(:,d2)));
+        else
+            farthestLeft = min(farthestLeft, min(Projection(c).proj(:,d1)));
+            farthestDown = min(farthestDown, min(Projection(c).proj(:,d2)));
+        end
+    end
+
+    plot(0,0,'b+', 'markerSi', 7.5);  % plot a central cross
+    
+    
+    %% if asked we will also plot the cross condition mean
+    if crossCondMean && length(Summary.crossCondMean) > 1
+        meanColor = [0 1 1];
+        
+        if isempty(overrideTimes)  % if we are going with the original times (those that were used to create the projection and do the analysis)        
+            P1 = Summary.crossCondMean(:,d1);
+            P2 = Summary.crossCondMean(:,d2);       
+        else
+            useTimes = ismember(Projection(c).allTimes, overrideTimes);
+            P1 = Summary.crossCondMeanAllTimes(useTimes,d1);
+            P2 = Summary.crossCondMeanAllTimes(useTimes,d2);
+        end
+        
+        plot(P1, P2, 'color', meanColor, 'lineWidth', 1.2*lineWidth);  % make slightly thicker than for rest of data.
+        if planMarkerSize>0
+            plot(P1(1), P2(1), 'ko', 'markerSize', planMarkerSize, 'markerFaceColor', meanColor);
+        end
+        
+        % for arrow, figure out last two points, and (if asked) supress the arrow if velocity is
+        % below a threshold.
+        penultimatePoint = [P1(end-1), P2(end-1)];
+        lastPoint = [P1(end), P2(end)];
+        vel = norm(lastPoint - penultimatePoint);
+
+        aSize = arrowSize + arrowGain * vel;  % if asked (e.g. for movies) arrow size may grow with vel
+        arrowMMC(penultimatePoint, lastPoint, [], aSize, axLim, meanColor, arrowEdgeColor);
+  
+    end
+
+    %% make axes
+    if useAxes
+        clear axisParams;
+
+        extraSeparation = axisSeparation*(min(farthestDown,farthestLeft));
+
+        % general axis parameters
+        axisParams.tickLocations = [-axisLength, 0, axisLength];
+        axisParams.longTicks = 0;
+        axisParams.fontSize = 10.5;
+
+        % horizontal axis
+        axisParams.axisOffset = farthestDown + extraSeparation;
+        axisParams.axisLabel = 'projection onto jPC_1 (a.u.)';
+        axisParams.axisOrientation = 'h';
+        haxP = AxisMMC(-axisLength, axisLength, axisParams);
+
+        % vertical axis
+        axisParams.axisOffset = farthestLeft + extraSeparation;
+        axisParams.axisLabel = 'projection onto jPC_2 (a.u.)';
+        axisParams.axisOrientation = 'v';
+        axisParams.axisLabelOffset = 1.9*haxP.axisLabelOffset;
+        vaxP = AxisMMC(-axisLength, axisLength, axisParams);
+    end
+
+    
+    % plot a label at the top
+    if substRawPCs == 0, planeType = 'jPCA'; else planeType = 'PCA'; end
+    if useLabel
+        if substRawPCs == 1
+            titleText = sprintf('raw PCA plane %d', plane);
+        elseif strcmp(rankType, 'varCapt')
+            letters = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ';
+            titleText = sprintf('jPCA plane %s (var capt ranked)', letters(planes2plot_Orig(pindex)));
+        else
+            titleText = sprintf('jPCA plane %d (eigval ranked)', plane);
+        end
+        titleText2 = sprintf('%d%% of var captured', round(100*Summary.varCaptEachPlane(plane)));
+        text(0,0.99*axLim(4),titleText, 'horizo', 'center');
+        text(0,0.88*axLim(4),titleText2, 'horizo', 'center', 'fontSize', 8.5);
+    end
+
+
+
+end  % done looping through planes
+
+end  % end of the main function
+
+
+
+