text:
%
% Tools for using nmds multidimensional scaling plots and data structure
% the data structure variable in this code is called c (arbitrarily)
% See cluster_nmdsfig for a function that uses many of these.
%
% By Tor Wager, Aug-Sept 2006
% V 1.2, April 2007
%
% ---------------------------------------------------------------------
% --
% Get data matrix, rank [optional], adjust for nuisance [optional]
% -----------------------------------------------------------------------
% c = nmdsfig_tools('get_data_matrix',c,cl,fldname,colm,contrasts,doranks)
% c = nmdsfig_tools('get_data_matrix',c,cl,'CONTRAST.avg',1,[],1);
%
% define input c.covs_nointerest with covariates to remove them.
% define input c.outcome to also remove covs from outcome.
%
% -----------------------------------------------------------------------
% get correlations, significance, and number of dimensions
% If groups are entered and ANCOVA is specified, handle that
% -----------------------------------------------------------------------
% c = nmdsfig_tools('get_correlations',c);
%
% if c.ancova_codes is input, does ancova analysis and returns stats
%
% -----------------------------------------------------------------------
% get cluster solution for a set of data (objectsByDimsData or
% GroupSpace)
% -----------------------------------------------------------------------
% c = cluster_solution(c, objectsByDimsData, n_cluster_range, niter, names)
%
% names is cell of one name per object
% to set up input, use get_correlations, and then this:
% [c.GroupSpace,c.obs,c.implied_dissim] = shepardplot(c.D,[]);
% c = nmdsfig_tools('cluster_solution',c, c.GroupSpace, 2:5, 1000, []);
%
% -----------------------------------------------------------------------
% Make a connectivity figure; with optional add btwn vars and side plots
% Uses stats from correlations
% -----------------------------------------------------------------------
% nmdsfig_tools('nmdsfig_plot',c, addbtwn, dosideplots, fillstr)
% addbtwn, dosideplots are 1 or 0, fillstr is 'fill' or something else
% uses colors in c.colors if it exists.
% -----------------------------------------------------------------------
% Make a connectivity figure; WITH outcome var and covariates on figure
% -----------------------------------------------------------------------
% DO THIS USING c STRUCTURE AFTER CLUSTER SOLUTION IS DONE
% nmdsfig_tools('nmdsfig_plot_withcovs',c, cl, dorank)
% needs fields: c.outcome c.covs_nointerest
%
% -----------------------------------------------------------------------
% Apply a classification to a set of data and r matrix, get class avgs
% -----------------------------------------------------------------------
% c = nmdsfig_tools('apply_clusters',c)
% c.class_avg_dat contains class averages
% needs c.ClusterSolution.classes, c.r
%
% Can turn off interactive naming
% c = nmdsfig_tools('apply_clusters',c, 0);
%
% -----------------------------------------------------------------------
% Get class avgs from some other data matrix using an existing c struct
% -----------------------------------------------------------------------
% class_avg_dat = nmdsfig_tools('get_class_avgs',c,dat)
% needs c.ClusterSolution.classes, dat (data matrix)
%
% -----------------------------------------------------------------------
% Predict behavior from either class average data or regions
% -----------------------------------------------------------------------
% c = nmdsfig_tools('predict_behavior',c,'classes');
% c = nmdsfig_tools('predict_behavior',c,'regions');
%
% -----------------------------------------------------------------------
% Add lines to a graph
% -----------------------------------------------------------------------
% handles = nmdsfig_tools('drawlines',coords,significance_mtx, bendval);
% specify pos/neg colors and styles:
% nmdsfig_tools('drawlines',c.GroupSpace,c_compare.sig,[1 0 0;0 0 1],{':',':'}, .05);
%
% -----------------------------------------------------------------------
% % Remove lines from a graph
% -----------------------------------------------------------------------
% nmdsfig_tools('removelines');
%
% -----------------------------------------------------------------------
% Add stars where stepwise params are significant
% -----------------------------------------------------------------------
% nmdsfig_tools('stepwise_stars',c,'classes')
% nmdsfig_tools('stepwise_stars',c,'regions')
%
% or where robust regression with behavior is significant
% nmdsfig_tools('robust_correl_stars',c,wh_reg)
%
% -----------------------------------------------------------------------
% Get unique colors for each object, but similar colors within a class
% -----------------------------------------------------------------------
% colors = nmdsfig_tools('class_colors',classes,basecolors)
%
%
% -----------------------------------------------------------------------
% Connect two points x and y in 3-D space with a straight or curved line
% -----------------------------------------------------------------------
% out = nmdsfig_tools('connect3d',x,y, varargin)
% x and y are [x,y,z] coordinate triples, or clusters/region objects with mm_center field.
% out.h is handles, out.coords are also returned
%
% Optional inputs: Keyword followed by value:
% color,thickness,bendpercent, nstreamlines, streamlineshift, nsamples
% bendpercent can be scalar or [x y z] triple for bend in each direction
% streamlineshift is multiplier for variation in x y z loc across streamlines
%
% Note: for 3-D graphs, use cluster_nmdsfig_glassbrain.m
%
% -----------------------------------------------------------------------
% Connect two points x and y in 2-D space with a straight or curved line
% -----------------------------------------------------------------------
% out = nmdsfig_tools('connect2d',x,y,color,thickness,bendpercent)
% x and y are [x,y] coordinate doubles, or clusters with mm_center field.
% out.h is handles, out.coords are also returned
% bendpercent can be scalar or [x y] pair for bend in each direction
% e.g.,
% out = nmdsfig_tools('connect2d', [-.3 0], [-.1 .1], 'b', 3, .1);
%
% SEQUENCES FOR ANALYSIS
% and examples
% -----------------------------------------------------------------------
% See also: cluster_nmdsfig, for a batch sequence
% -----------------------------------------------------------------------
% Use average data across conditions to get basic structure:
% c = nmdsfig_tools('get_data_matrix',c,cl,'CONTRAST.avg',1,[],1);
% c = nmdsfig_tools('get_correlations',c);
% nmdsfig_tools('nmdsfig_plot',c, 0, 0);
% c = nmdsfig_tools('apply_clusters',c);
% c = nmdsfig_tools('predict_behavior',c,'regions');
% c = nmdsfig_tools('predict_behavior',c,'classes');
%
% Get new data and impose the class structure on it:
% c = nmdsfig_tools('get_data_matrix',c,cl,'CONTRAST.data',1:2,[1 -1],0);
% c = nmdsfig_tools('predict_behavior',c,'regions');
% c = nmdsfig_tools('get_correlations',c);
% c = nmdsfig_tools('apply_clusters',c);
% c = nmdsfig_tools('predict_behavior',c,'classes')
%
% Create a basic figure based on the averages,
% and then use new data to get lines
% add stars where there are sig. correls with behavior
% c = nmdsfig_tools('get_data_matrix',c,cl,'CONTRAST.avg',1,[],1);  % ranks, for robustness
% c = nmdsfig_tools('get_correlations',c);
%f1 = nmdsfig(c.GroupSpace,'classes',c.ClusterSolution.classes,'sig',c.STATS.sigmat,'names',c.names,'sizescale',[4 16]);
% nmdsfig_tools('removelines');
% c = nmdsfig_tools('get_data_matrix',c,cl,'CONTRAST.data',1:2,[1 -1],0);
% c = nmdsfig_tools('get_correlations',c);
% nmdsfig_tools('drawlines',c.GroupSpace,c.STATS.sigmat);
% nmdsfig_tools('robust_correl_stars',c,1);
%
% Get correl diffs among class averages
% get_data_matrix adjusts for covs of no interest, ranks, and centers
%c = nmdsfig_tools('get_data_matrix',c,cl,'CONTRAST.data',2,[],0); yON = c.dat; % get rank data
%c = nmdsfig_tools('get_data_matrix',c,cl,'CONTRAST.data',1,[],0); yOFF = c.dat; % get rank data
%yONavgs = nmdsfig_tools('get_class_avgs',c,yON);
%yOFFavgs = nmdsfig_tools('get_class_avgs',c,yOFF);
% ycomp_struct = correl_compare_dep(yONavgs,yOFFavgs,.05,0); % data already
% %ranked, so no need to redo
%
%
% Complete example for pre-appraisal mediation data:
% For set-up, see mediation_brain, mediation_brain_results and
% mediation_brain_results_detail
%
%     c = []; c.covs_nointerest = [SETUP.X SETUP.covariates]; c.outcome = SETUP.Y; doranks = 1;
%     for i = 1:length(cl), c.names{i} = cl(i).shorttitle; end
%     c = nmdsfig_tools('get_data_matrix',c,cl,'timeseries',1,[],doranks);
%     c = nmdsfig_tools('get_correlations',c);
%     [c.GroupSpace,c.obs,c.implied_dissim] = shepardplot(c.D,[]);
%     c = nmdsfig_tools('cluster_solution',c, c.GroupSpace, 2:5, 1000, []);
%     c = nmdsfig_tools('apply_clusters',c);
%     c = nmdsfig_tools('predict_behavior',c,'classes');
%     cluster_orthviews_classes(cl,c.ClusterSolution.classes, EXPT.overlay, 'sagittal', 0);
%     nmdsfig_tools('nmdsfig_plot',c, 0, 0, 'fill');
% c_complete = nmdsfig_tools('nmdsfig_plot_withcovs',c, cl, doranks)

function [c, varargout] = nmdsfig_tools(meth,varargin)

c = [];

% First define required inputs
% -----------------------------------------------------------------------
switch meth

    case 'get_data_matrix'
        % This should execute the function:
        % c = get_data_matrix(c,cl,fldname,colm,contrasts,doranks)
        innames = {'c' 'cl' 'fldname' 'colm' 'contrasts' 'doranks'};

    case 'nmdsfig_plot'
        % nmdsfig_plot(c, addbtwn, dosideplots)
        innames = {'c' 'addbtwn' 'dosideplots' 'fillstr'};

    case 'nmdsfig_plot_withcovs'
        innames = {'c' 'cl' 'doranks'};

    case {'apply_clusters','get_correlations'}
        % c = apply_clusters(c)
        innames = {'c', 'dointeractive', 'corrtype'};
        if length(varargin) < 2, varargin{2} = []; end % optional; do interactive 1/0
        if length(varargin) < 3, varargin{3} = 'r'; end % optional; corrtype, 'r' = Pearson't (default)

    case {'cluster_solution'}
        innames = {'c' 'objectsByDimsData' 'n_cluster_range' 'niter' 'names'};

        % c = cluster_solution(c, objectsByDimsData, n_cluster_range, niter, names);

    case 'get_class_avgs'
        % class_avg_dat = nmdsfig_tools('get_class_avgs',c,dat)
        innames = {'c' 'dat'};

    case 'predict_behavior'
        %c = predict_behavior(c,meth)
        innames = {'c' 'predmeth'};

    case 'drawlines'
        % [hhp,hhn] = drawlines(pc,sigmat,sigcol,legmat,[color],[style])

        innames = {'pc','sigmat','color','style','bendval'};
        % optional: color/style format strings
        if length(varargin) < 3, varargin{3} = []; end % optional
        if length(varargin) < 4, varargin{4} = []; end % optional
        if length(varargin) < 5, varargin{5} = []; end % optional
        if length(varargin) > 5, error('Too many input arguments.'); end

    case 'removelines'
        innames = {};

    case 'stepwise_stars'
        innames = {'c' 'meth'};

    case 'robust_correl_stars'
        innames = {'c' 'wh_reg'};

    case 'class_colors'
        innames = {'classes' 'basecolors'};
        if length(varargin) < 2, varargin{2} = []; end % optional

    case 'connect3d'
        % h = nmdsfig_tools('connect3d',x,y,color,thickness,bendpercent,[nsamples])
        % innames = {'x' 'y' 'color','thickness','bendpercent','nsamples'};
        % This has been updated to use inputParser object (2019)
        innames = varargin;
        %if length(varargin) < 6, varargin{6} = []; end % optional
        % Skip all the stuff below
        c = connect3d(varargin{:});
        return

    case 'connect2d'
        % out = connect2d(x,y,color,thickness,bendpercent,varargin)
        innames = {'x' 'y' 'color','thickness','bendpercent','nsamples'};
        if length(varargin) < 6, varargin{6} = []; end % optional

    otherwise
        error('Unknown method string.')
end

% this generic code block reads varargin into variables for
% do it this way partly to keep track of vars by name
% -----------------------------------------------------------------------
if length(varargin) < length(innames)
    error('Wrong number of input arguments.  Check help on method you are using.')
end
for i = 1:length(innames)
    eval([innames{i} ' = varargin{' num2str(i) '};']);
end

% Now execute subfunction
% -----------------------------------------------------------------------
str = ['c = ' meth '('];
for i = 1:length(varargin)
    if i ~= 1, str = [str ',']; end
    str = [str innames{i}];
end
str = [str ');'];

%disp(str)
eval(str)

return

% -----------------------------------------------------------------------
% -----------------------------------------------------------------------
%
%
% Sub-functions
%
%
% -----------------------------------------------------------------------
% -----------------------------------------------------------------------

% -----------------------------------------------------------------------
% -----------------------------------------------------------------------
% get data from selected field
% -----------------------------------------------------------------------
% -----------------------------------------------------------------------

function c = get_data_matrix(c,cl,fldname,colm,contrasts,doranks)

if ~isfield(c, 'outcome')
    c.outcome = [];
end

if ~isfield(c, 'covs_nointerest')
    c.covs_nointerest = [];
end

dat = [];
for i = 1:length(cl)
    str = ['dat1 = cl(i).' fldname '(:,[' num2str(colm) ']);'];
    eval(str)

    if ~isempty(contrasts), dat1 = dat1 * contrasts'; end

    if doranks, dat1 = rankdata(dat1); end
    dat=[dat dat1];
end

if doranks
    disp('Returning rank data in c.dat');
    c.dat_descrip = 'Dat and outcome data: Ranked';

    if ~isempty(c.outcome)
        disp('Ranking outcome data.')
        c.outcome = rankdata(c.outcome);
    end

    if ~isempty(c.covs_nointerest)
        disp('Ranking covariate data.')
        for j = 1:size(c.covs_nointerest, 2)
            c.covs_nointerest(:,j) = rankdata(c.covs_nointerest(:,j));
        end
    end

else
    disp('Returning data in c.dat');
    c.dat_descrip = 'Dat and outcome data: Original (not ranked) ';
end

c.doranks = doranks;
%if length(varargin) > 0,dat = [c.outcome dat];end

% remove covariates of no interest, if any
% if data is not rank data, use robust regression
% otherwise, use OLS
% -----------------------------------------------------------------------

if isempty(c.covs_nointerest)
    % do nothing
    disp('nmdsfig_tools get_data_matrix : no covariates entered.')
    c.dat_descrip = [c.dat_descrip ' raw (not adjusted)'];
else
    disp('nmdsfig_tools get_data_matrix : Removing covariate(s) of no interest (using IRLS).')
    c.dat_descrip = [c.dat_descrip ' Adjusted for covariates with IRLS'];
    if ~isfield(c, 'outcome'), c.outcome = []; end

    if ~isempty(c.outcome)
        % partialcorr gets adjusted x, y, and correls; here just use it to
        % get adjusted x
        %c.outcome  = partialcor([c.outcome c.covs_nointerest],ones(size(c.outcome,1),1),1);
        c.outcome =  adjust_y(c.covs_nointerest,c.outcome,doranks);
    end

    for i = 1:size(dat,2)
        %dat(:,i) = partialcor([dat(:,i) c.covs_nointerest],ones(size(dat,1),1),1);
        [dat(:,i), b(:, i)] =  adjust_y(c.covs_nointerest,dat(:,i),doranks);
    end

    c.nuisance_betas = b;
end

c.dat = dat;

return

% -----------------------------------------------------------------------
% -----------------------------------------------------------------------

% % Adjust x and y for regressors of no interest

% -----------------------------------------------------------------------
% -----------------------------------------------------------------------
function [r, b] = adjust_y(x,y,doranks)

if doranks
    X = [x ones(size(x,1),1)];
    b = pinv(X) * y;
    r = y - X * b;
else
    [b,stats]=robustfit(x,y);
    r = stats.resid;
end

return

% -----------------------------------------------------------------------
% -----------------------------------------------------------------------

% Make an NMDSfig figure from a c structure

% -----------------------------------------------------------------------
% -----------------------------------------------------------------------

function figname = nmdsfig_plot(c, addbtwn, dosideplots, fillstr)

% Required fields: (See cluster_nmdsfig, for example)
% requires c. ...
% STATS.sigmat/sigmat2/siginteract
% outcome
% ClusterSolution.classes
% names
% GroupSpace (stim. coords)
% X
% r
%
% btwn_coords (for addbtwn option)
% factor_scores (for sideplots option)

f1 = create_figure;

colors = {'yo'  'bo' 'go' 'ro' 'co' 'mo' 'ko' 'r^' 'g^' 'b^' 'y^' 'c^' 'm^' 'k^'};

if isfield(c,'colors'), colors = c.colors; end

% for side plots
%subplot('position',[0.3 0.3 0.5 0.6]);
if ~isfield(c,'outcome'), c.outcome = []; end

if ~isempty(c.outcome)  % plot interactions between beh and regional covariance
    if isfield(c.STATS,'siginteract')
        % we have a group x cov interaction in ANCOVA
        nmdsfig(GroupSpace,'classes',c.ClusterSolution.classes,'names',c.names,'sig',c.STATS.sigmat, ...
            'legend',{'Pos' 'Neg' 'Pos x Beh' 'Neg x Beh'},'sig2',c.STATS.siginteract, ...
            'colors', colors, fillstr);
        figname = 'cluster_nmdsfig_ancova_interact';
    else
        nmdsfig(c.GroupSpace,'classes',c.ClusterSolution.classes,'names',c.names,'sig',c.STATS.sigmat,'sig2', ...
            c.STATS.sigmat2,'legend',{'Pos' 'Neg'},'sizescale',[4 16], ...
            'colors', colors, fillstr);
        figname = 'cluster_nmdsfig_mdsfig';

    end
else
    nmdsfig(c.GroupSpace,'classes',c.ClusterSolution.classes,'names',c.names,'sig',c.STATS.sigmat,'sig2', ...
        c.STATS.sigmat2,'legend',{'Pos' 'Neg'},'sizescale',[4 16], ...
        'colors', colors, fillstr);
    figname = 'cluster_nmdsfig_mdsfig';
end

try
    nmdsfig_legend(c.GroupSpace,c.r);
catch
end

save_figure('cluster_nmdsfig_legend',550);

% behavior
if addbtwn
    if ~isfield(c, 'btwn_coords'), error('You must enter c.btwn_coords to use btwn-subjects breakdown plot.  See cluster_nmdsfig.'); end

    hold on; plot(c.btwn_coords(1),c.btwn_coords(2),'ko','MarkerSize',12,'LineWidth',3);
end

% MDSfig sideplots
% -----------------------------------------------------------------------

dosideplots = 0;

if dosideplots

    % X-axis subplot: Component 1
    subplot('position',[0.25 0.1 0.6 0.1]);     hold on;
    set(gca,'FontSize',16);
    fact = c.factor_scores(:,1);
    if ~isempty(c.outcome),
        fact=[fact c.outcome]; fact=sortrows(fact,2);
        xlab='Activation score for each participant: Component 1';

        lo = find(fact(:,2)<median(fact(:,2)));
        hi = find(fact(:,2)>median(fact(:,2)));

        plot(lo,fact(lo,1),'ko','LineWidth',2);
        plot(hi,fact(hi,1),'k^','LineWidth',2);
        set(gca,'XTick',[1 max(lo)+.5 size(fact,1)],'XTickLabel',{'Lowest' 'Median' 'Highest'})
        title(xlab)
        xlabel('Behavior');
        ylabel('Score')
    else
        xlab='Scores on Component 1';
        plot(fact(:,1),'ko','LineWidth',2);
        title(xlab); xlabel('Index'); ylabel('Score');
    end

    % Y-axis subplot: Component 2
    subplot('position',[0.1 0.3 0.1 0.6]);   hold on;
    set(gca,'FontSize',16);
    fact = c.factor_scores(:,2);
    if ~isempty(c.outcome)
        fact=[fact c.outcome]; fact=sortrows(fact,2);
        xlab='Component 2';

        lo = find(fact(:,2)<median(fact(:,2)));
        hi = find(fact(:,2)>median(fact(:,2)));

        plot(fact(lo,1),lo,'ko','LineWidth',2);
        plot(fact(hi,1),hi,'k^','LineWidth',2);
        set(gca,'YTick',[1 max(lo)+.5 size(fact,1)],'YTickLabel',{'Lowest' 'Median' 'Highest'})
        title(xlab)
        ylabel('Behavior');
        xlabel('Score')
    else
        xlab='Scores on Component 1';
        plot(fact(:,1),'ko','LineWidth',2);
        title(xlab); ylabel('Index'); xlabel('Scores: Comp. 2')
    end

end % if dosideplots

figure(f1)
axis off
drawnow

save_figure('cluster_nmdsfig',500);

return

% -----------------------------------------------------------------------
% -----------------------------------------------------------------------

% Re-plot an nmds figure with outcome and covs added (but not part of
% clustering solution); Use existing completed c structure

% -----------------------------------------------------------------------
% -----------------------------------------------------------------------
function c_complete = nmdsfig_plot_withcovs(c, cl, doranks)

extravars = [c.outcome c.covs_nointerest];
nextra = size(extravars, 2);

% get data again, and do not remove covariates
c_complete = c;
c_complete = rmfield(c_complete, 'covs_nointerest');
c_complete = nmdsfig_tools('get_data_matrix',c_complete,cl,'timeseries',1,[],doranks);

disp('Creating c_complete.dat with [Outcome Covs Vars]');
c_complete.dat = [extravars c_complete.dat];

c_complete = nmdsfig_tools('get_correlations',c_complete);
[c_complete.GroupSpace,c_complete.obs,c_complete.implied_dissim,c_complete.stress] = shepardplot(c_complete.D,[]);

c_complete.ClusterSolution = c.ClusterSolution;

cl_id = max(c_complete.ClusterSolution.classes) + 1;  % class ID for extra points

c_complete.ClusterSolution.classes = [cl_id .* ones(nextra, 1); c_complete.ClusterSolution.classes];

extranames = {};
if ~isempty(c.outcome), extranames{1} = 'Outcome'; end
if ~isempty(c.covs_nointerest)
    for i = 1:size(c.covs_nointerest, 2)
        extranames{end + 1} = ['Cov' num2str(i)];
    end
end

c_complete.names = [extranames c.names];

nmdsfig_tools('nmdsfig_plot',c_complete, 0, 0, 'nofill');

saveas(gcf,'NMDS_complete_flat_plot_orig','png');

% direct effects, but should exclude OUTCOME?
[c_complete.direct_mtx, c_complete.mediated_mtx, c_complete.mediators] = matrix_direct_effects(c_complete.STATS.sigmat2, c_complete.dat);

nmdsfig_tools('removelines');
nmdsfig_tools('drawlines',c_complete.GroupSpace, c_complete.direct_mtx);
saveas(gcf,'NMDS_complete_flat_plot_direct','png');

return

% -----------------------------------------------------------------------
% -----------------------------------------------------------------------

% Get correlation matrix, depending on ancova status

% -----------------------------------------------------------------------
% -----------------------------------------------------------------------

function c = get_correlations(c, varargin)

disp('Saving correlation matrix in c.r')

if ~isfield(c, 'ancova_codes'), c.ancova_codes = []; end

corrtype = 'r';
if length(varargin) > 1, corrtype = varargin{2}; end
if ~strcmp(corrtype, 'r') && ~isempty(c.ancova_codes)
    warning('nmdsfig_tools:incompat_inputs', 'Corr type cannot be entered with ANCOVA covariates.');
end

disp(['Correlation/association metric: ' corrtype])

if isempty(c.ancova_codes)
    doancova = 0;
    % no between covariates, just use correlations
    if strcmp(corrtype, 'r')
        [r,p]=corrcoef(c.dat); r(r>.9999) = 1;
    elseif strcmp(lower(corrtype), 'spearman')
        [r,p] = corr(c.dat, 'type', 'Spearman');
    else
        [r,t,p] = correlation(corrtype, c.dat);
    end

    c.STATS.corrtype = corrtype;

    [pthr,sig] = fdr_correct_pvals(p,r);

    c.STATS.p_thr = pthr;

else
    %doancova = input('Use covariate scores as covariate in ANCOVA? (inter-region correls then done w/i group): 1/0: ');
    doancova = 1;
end

if doancova

    disp('Using ANCOVA to model between-subjects effects')
    disp('Inter-region correlations within groups; interaction tests group x brain covariance interaction')
    % main effect of grp in this case is hi vs lo behavior on region Y (brain act.) after
    % controlling for region X -- a bit hard to interpret

    if length(unique(c.ancova_codes) > 2),
        disp('Binarizing ANCOVA code regressor into two groups.');
        c.ancova_codes = mediansplit(c.ancova_codes);
    end

    [b,t,p] = ancova(c.ancova_codes,[],c.dat); % first cell is grp (diff), 2nd is r, 3rd is r diff x group
    pall = p;   % save all slopes

    r = b{2};   % overall slope
    r = r + eye(size(r));   % add 1 to diagonal
    p = pall{2};
    [c.STATS.pthr(1),sig] = fdr_correct_pvals(p,r);

    ri = b{3};  % behavior x covariance interaction
    %ri = ri + eye(size(ri));   % add 1 to diagonal
    pi = pall{3};
    [c.STATS.pthr(2),sigi] = fdr_correct_pvals(pi,ri);

    c.STATS.siginteract = sigi;
    c.STATS.b = b; c.STATS.t = t; c.STATS.p = pall;
    c.STATS.btp_descrip = 'first cell is grp (diff), 2nd is r, 3rd is r diff x group.  betas, t-vals, p-vals';

else
    c.STATS.p = p;
end

c.doancova = doancova;
c.r = r;
c.STATS.sigmat = sig;
c.STATS.sigmat2 = (p < .05) .* sign(r);
c.STATS.p_descrip = 'sig. p-values are FDR corrected.';
c.STATS.p_descrip2 = 'sigmat2 is p < .05 uncorrected.';

fprintf('Positive connections (q < .05 FDR): %3.0f\n', sum(c.r(find(c.STATS.sigmat(:))) > 0))
fprintf('Negative connections (q < .05 FDR): %3.0f\n', sum(c.r(find(c.STATS.sigmat(:))) < 0))

disp('Getting dissimilarity matrix from correlation matrix, saving in c.D')
% scale correlations (or beta weights) so that they are a similarity matrix between 0 and 1.
c.D = (1 - c.r) ./ 2;
%c.D=distosim(c.D);
% Make sure it's a dissim. matrix
c.D=simtodis(c.D);

return

function [pthr,sig] = fdr_correct_pvals(p,r)

psq = p; psq(find(eye(size(p,1)))) = 0;
psq = squareform(psq);
pthr = FDR(p,.05);
if isempty(pthr), pthr = 0; end

sig = sign(r) .* (p < pthr);

return

% -----------------------------------------------------------------------
% -----------------------------------------------------------------------

% Find cluster solution

% -----------------------------------------------------------------------
% -----------------------------------------------------------------------

function c = cluster_solution(c, objectsByDimsData, n_cluster_range, niter, names)
% -----------------------------------------------------------------------
% testcluster:
%
% Permute objects in space, get null hypothesis cluster quality, and test
% observed cluster solution against this.
% -----------------------------------------------------------------------

% pval is p-values for each number of clusters tested
% classes = class ("network") assignments for best clustering solution

c.ndims = size(objectsByDimsData, 2);
c.GroupSpace = objectsByDimsData;

[c.ClusterSolution.pvals, ...
    c.ClusterSolution.classes, ...
    c.ClusterSolution.classnames, ...
    c.ClusterSolution.X, ...
    c.ClusterSolution.include, ...
    c.ClusterSolution.names, ...
    c.ClusterSolution.stats] = ...
    testclustnew(objectsByDimsData, n_cluster_range, c.ndims, niter, names, 'keep', 'average');

c.ClusterSolution.niter = niter;

return

% -----------------------------------------------------------------------
% -----------------------------------------------------------------------

% Apply clusuter solution and get class average data

% -----------------------------------------------------------------------
% -----------------------------------------------------------------------

function c = apply_clusters(c, varargin)

dointeractive = 1;
if length(varargin) > 0, dointeractive = varargin{1}; end
%
% requires c. ...
% ClusterSolution.classes
%
% uses c. ...
% outcome
%
% creates
% c.ClusterSolution.classnames
% c.class_avg_data

% add names of regions if missing
if ~isfield(c, 'names') || isempty(c.names) || length(c.names) < size(c.dat, 2)
    disp('Missing or invalid region names; adding.');
    for ii = 1:size(c.dat, 2)
        c.names{ii} = ['V' num2str(ii)];
    end
end

% name 'networks'
% -----------------------------------------------------------------------

if ~isfield(c,'ClusterSolution'), error('You need a ClusterSolution field first! ... try testclustnew'), end

if ~isfield(c,'dat'), error('Enter data for each variable in .dat field first'), end

if ~isfield(c,'colors') || isempty(c.colors)
    %c.colors = {'yo' 'bo' 'go' 'ro' 'co' 'mo' 'ko' 'r^' 'g^' 'b^' 'y^' 'c^' 'm^' 'k^'};
    c.colors = {[1 1 0] [0 0 1] [0 1 0] [1 0 0] [1 .5 0] [0 1 1] [1 0 1] [.5 1 0] [0 .5 1] [0 1 .5] [1 0 .5]};
end

if ~isfield(c,'GroupSpace')
    disp(c)
    fieldnm = input('Type name of GroupSpace field with MCA/MDS/PCA coordinates: ', 's');
    myGroupSpace = c.(fieldnm);
else
    myGroupSpace = c.GroupSpace;
end

if ~isfield(c,'doancova')
    c.doancova = 0;
end

classnames = [];
if isfield(c,'APPLY_CLUSTER') && isfield(c.APPLY_CLUSTER,'names')
    classnames = c.APPLY_CLUSTER.names;
end

if (isempty(classnames) || length(classnames) < max(c.ClusterSolution.classes)) && (isempty(dointeractive) || dointeractive)
    donames = input('Name Clusters (Classes)? (1/0) ');

    if donames
        for i = 1:max(c.ClusterSolution.classes)
            disp(['Network ' num2str(i)])
            tmp = find(c.ClusterSolution.classes == i);
            for j = 1:length(tmp)
                fprintf(1,'%s . ',c.names{tmp(j)});
                if j > 4, fprintf(1,'\n'),end
            end
            fprintf(1,'\n')
            classnames{i} = input('Name this network: ','s');
        end

    else
        for i = 1:max(c.ClusterSolution.classes)
            classnames{i} = ['Net' num2str(i)];
        end
    end
end

% if we have a between-subjects covariate, enter it, otherwise empty
if ~isfield(c,'outcome'), c.outcome = []; end

% Apply this clustering solution to the clusters
% Get individual 'network' scores and use them as predictors of behavior
% in multiple regression
% -----------------------------------------------------------------------

% now run for the average across task conditions
% average correlations within and between networks, and test
if isfield(c, 'r')
    mycorrmtx = c.r; disp('Using .r field for correlation matrix.');
elseif isfield(c, 'corr')
    mycorrmtx = c.corr; disp('Using .corr field for correlation matrix.');
else
    error('Need .r or .corr field in input structure.')
end

c.APPLY_CLUSTER = apply_cluster_solution(c.ClusterSolution.classes,...
    mycorrmtx,...
    'names',classnames,'bcov',c.outcome, 'n', size(c.dat, 1), 'dointeractive', dointeractive);

c.APPLY_CLUSTER.names = classnames;

% matrix plot is not meaningful here; close
close;
%save_figure('cluster_nmdsfig_rmatrix',400);

% standardize variables?  No, let's let low variance vars contribute less
% to avgs

c.class_avg_dat = get_class_avgs(c,c.dat);
[c.APPLY_CLUSTER.avgr,c.APPLY_CLUSTER.avgp] = corrcoef(c.class_avg_dat);
c.APPLY_CLUSTER.avgsig = c.APPLY_CLUSTER.avgp < .055 .* sign(c.APPLY_CLUSTER.avgr);

% get class (network) average stimulus locations (group space) and plot
% -----------------------------------------------------------------------

Gs = myGroupSpace(:,1:c.ndims);
cla = c.ClusterSolution.classes;
for i = 1:max(cla)
    wh = find(cla == i);
    if length(wh)>1;
        classGs(i,:) = mean(Gs(wh,:));
    else
        classGs(i,:) = Gs(wh,:);
    end
end
c.APPLY_CLUSTER.Gs = classGs;
c.APPLY_CLUSTER.classes = 1:max(cla);

% significance of r among class averages (done above)
% [r,p]=corrcoef(c.class_avg_dat); r(r>.9999) = 1;
% sig = (p < .05) .* sign(r);
%[pthr,sig] = fdr_correct_pvals(p,r);

% Figures
if strcmp(c.doancova,'Yes'), create_wide_figure; subplot(1,2,1); else create_figure; end
set(gca,'FontSize',14)
nmdsfig(classGs(:,1:2),'names',c.APPLY_CLUSTER.names,'classes',c.APPLY_CLUSTER.classes, ...
    'sig', c.APPLY_CLUSTER.avgsig, ... %c.APPLY_CLUSTER.group_stats.sigu,'sizescale',[4 16]);
    'sizescale',[4 16],'colors',c.colors);

title(['Correlations among class averages (p <= .05 two-tailed)'])

if strcmp(c.doancova,'Yes'), subplot(1,2,2);
    set(gca,'FontSize',14)
    mdsfig(classGs(:,1:2),classnames,cla,c.STATS.siginteract);
    title(['Group by cov interaction: ANCOVA'])
end

save_figure('cluster_nmdsfig_class_lines_uncor',550);

return

% -----------------------------------------------------------------------
% -----------------------------------------------------------------------

% Get class average data

% -----------------------------------------------------------------------
% -----------------------------------------------------------------------

function class_avg_dat = get_class_avgs(c,dat)

class_avg_dat = [];

[nobs,nvars] = size(dat);
nclasses = length(c.ClusterSolution.classes);
if nvars ~= nclasses, error('Number of classes does not match num. columns in dat'); end

for i = 1:max(c.ClusterSolution.classes)
    wh = find(c.ClusterSolution.classes == i);
    if length(wh) > 1
        class_avg_dat = [class_avg_dat nanmean(dat(:,wh)')'];
    else
        class_avg_dat = [class_avg_dat dat(:,wh)];
    end
end

return

% -----------------------------------------------------------------------
% -----------------------------------------------------------------------

% Stepwise regressions on behavior

% -----------------------------------------------------------------------
% -----------------------------------------------------------------------

function c = predict_behavior(c,meth)
% needs:
% c.
% outcome
% meth is 'classes' or 'regions'

if isempty(c.outcome)
    disp('Between-subjects c.outcome is empty.  No correlations with behavior run.');
    return
end

switch meth
    case 'classes'

        % Add logistic regression for categorical 1/0 DVs!
        % see also stepwise_tor.m which does this stuff.
        disp(' ')
        disp('------------------------------------------')
        disp('Using networks to predict behavior')
        disp('------------------------------------------')
        disp(' ')
        str = correlation_to_text([c.outcome c.class_avg_dat],[],[{'Behav'} c.APPLY_CLUSTER.class_names]);

        %[c.ClusterSolution.STEPWISE.b,c.ClusterSolution.STEPWISE.se,c.ClusterSolution.STEPWISE.pval, ...
        %c.ClusterSolution.STEPWISE.inmodel,stats] = ...
        %stepwisefit(c.class_avg_dat,c.outcome);
        disp(' ')
        c.class_STEPWISE = stepwise_tor(c.class_avg_dat,c.outcome,c.APPLY_CLUSTER.class_names);
        STEP = c.class_STEPWISE;
        names = c.APPLY_CLUSTER.class_names;

    case 'regions'

        %N = stats.dfe+stats.df0;
        %r2 = (stats.SStotal-stats.SSresid) ./ stats.SStotal;
        %adjr2 = 1 - (1-r2)*((N - 1) ./ (N - stats.df0 - 1));
        %fprintf(1,'Omnibus F(%3.0f,%3.0f) = %3.2f, RMSE = %3.2f, p = %3.4f, Adj R^2 = %3.2f\n', ...
        %stats.df0,stats.dfe,stats.fstat, ...
        %stats.rmse,stats.pval,adjr2);
        %stats.adjr2 = adjr2;
        %c.ClusterSolution.STEPWISE.stats = stats;

        disp(' ')
        disp('------------------------------------------')
        disp('Using individual regions to predict behavior'),
        disp('------------------------------------------')
        disp(' ')
        str = correlation_to_text([c.outcome c.dat],[],[{'Behav'} c.names],1);
        disp(' ')
        c.indiv_STEPWISE = stepwise_tor(c.dat,c.outcome,c.names);
        STEP = c.indiv_STEPWISE;
        names = c.names;

    otherwise error('Unknown method.  Valid = classes or regions')
end

% ------------------------------------------
% make plots
% ------------------------------------------
doplot = 1;

if ~isfield(c,'colors')
    c.colors = {'yo' 'bo' 'go' 'ro' 'co' 'mo' 'ko' 'r^' 'g^' 'b^' 'y^' 'c^' 'm^' 'k^'};
end
while length(c.colors) < size(c.dat, 2), c.colors = [c.colors c.colors]; end

wh = logical(STEP.inmodel);

switch meth
    case 'classes'
        X = c.class_avg_dat(:, wh);
    case 'regions'
        X = c.dat(:, wh);
end

names = names(:, wh);

colors_to_use = c.colors(wh);

nsig = size(X,2);

if doplot && nsig > 0

    disp('Plotting partial correlation scatterplots for significant variables.');

    % Check for existing figure, and create if necessary
    create_figure('class_scatterplots', 1, nsig);

    scn_export_papersetup(300);

    for i = 1:nsig
        subplot(1, nsig, i);

        prplot(c.outcome, X, i, 0, colors_to_use(i));
        title(names{i})
        xlabel(' ');
        ylabel(' ');

    end

end

return

% -----------------------------------------------------------------------
% -----------------------------------------------------------------------

% Draw lines on nmdsfig graph

% -----------------------------------------------------------------------
% -----------------------------------------------------------------------

function line_handles = drawlines(pc,sigmat,varargin)
% function handles_out = drawlines(pc,sigmat,[color],[style],[bendval])
%
% pc is stim coords
% sigmat is matrix of which lines to draw
% sigcol is multiplier to scale line width (0 is fixed width)
%   sigcol is a scalar between zero and one, lower is 'more salient'
% legmat is not used now (for legend stuff)
%
% Example:
% [hh3,hh4] = drawlines(pc,wh,sigcol(i),legmat,[.5 .5 .5; 0 1 1],[{':'} {':'}]);

% sigcol is a scalar between zero and one, lower is 'more salient'
sigcol = 0;
hold on

% linew
lw = 2;

color(1,:) = [0 0 0];   % first line color, 'positive'
color(2,:) = [0 1 1];   % second line color, 'negative'
style{1} = '-';         % first line style
style{2} = '-';         % second line style

if length(varargin) > 0 && ~isempty(varargin{1}), color = varargin{1}; end
if length(varargin) > 1 && ~isempty(varargin{2}), style = varargin{2}; sigcol = 0; end

bendval = .1;
if length(varargin) > 2 && ~isempty(varargin{3}), bendval = varargin{3};  end

hhp=[];hhn=[];
for i = 1 : size(pc,1)
    for j = i+1 : size(pc,1)

        if sigmat(i,j) > 0 || sigmat(i,j) < 0
            out = nmdsfig_tools('connect2d', [pc(i,1) pc(i,2)], [pc(j,1) pc(j,2)], 'b', lw - (lw*sigcol)-1, bendval);

            if sigmat(i,j) > 0
                hhp(end+1) = out.h;
                %hhp(end+1) = line([pc(i,1) pc(j,1)],[pc(i,2) pc(j,2)]);
                set(hhp,'Color',color(1,:) + [1 1 1] * abs(sigcol),'LineStyle',style{1},'LineWidth',lw - (lw*sigcol)-1)
            elseif sigmat(i,j) < 0
                hhn(end+1) = out.h;
                %hhn(end+1) = line([pc(i,1) pc(j,1)],[pc(i,2) pc(j,2)]);
                set(hhn,'Color',color(2,:) + abs([sigcol sigcol 0]),'LineStyle',style{2},'LineWidth',lw - (lw*sigcol)-1)
            end

        end
    end
end
%legend ([hhp hhn],legmat);

% store in gui data for later deletion, etc.
figh = gcf;
linehandles = [hhp hhn];
data = guidata(figh);
if ~isfield(data,'linehandles')
    data.linehandles = linehandles;
else
    data.linehandles = [data.linehandles linehandles];
end
guidata(figh,data);

line_handles.hhp = hhp;
line_handles.hhn = hhn;

return

% -----------------------------------------------------------------------
% -----------------------------------------------------------------------

% Remove lines from nmdsfig plot

% -----------------------------------------------------------------------
% -----------------------------------------------------------------------
function doneok = removelines

doneok = 0;
figh = findobj('Tag', 'nmdsfig');
data = guidata(figh);
if isfield(data,'linehandles') && ~isempty(data.linehandles)
    delete(data.linehandles(ishandle(data.linehandles)))
    data.linehandles = [];
    guidata(figh,data);
    doneok = 1;
else
    disp('Cannot find line handles to remove in guidata.  Doing nothing.')
end

return

% -----------------------------------------------------------------------
% -----------------------------------------------------------------------

% Put stars on figure where there are sig. relations in stepwise regression
% also do stars for robust correlations in glm (sep. function)
% -----------------------------------------------------------------------
% -----------------------------------------------------------------------
function sig = stepwise_stars(c,meth)

switch meth
    case 'regions'
        myfield = 'indiv_STEPWISE';
        myf2 = 'GroupSpace';

    case 'classes'
        myfield = 'class_STEPWISE';
        myf2 = 'APPLY_CLUSTER.Gs';

    otherwise
        error('Unknown method.  regions or classes.')
end

sig = c.(myfield).inmodel' .* c.(myfield).t;
plot_stars(c,sig,myf2);

return

function sig = robust_correl_stars(c,wh_reg)

[b,t,p,sig,F,fp,fsig,stat] = robust_reg_matrix([c.X ones(size(c.X,1),1)],c.dat,1);
sig = sig(wh_reg,:)' .* t(wh_reg,:)';
plot_stars(c,sig,'GroupSpace');

return

function plot_stars(c,sig,myfield)

wh = find(sig);
nsig = length(wh);

if nsig == 0
    disp('No significant results to plot.')
    return
end

% coordinates
pc = c.(myfield);
pc = pc(wh,1:2);
signval = (sig(wh) > 0) + 1;    % 1 for negative, 2 for positive
colors = {'bp' 'rp'};

for i = 1:nsig
    plot(pc(i,1),pc(i,2),colors{signval(i)},'MarkerSize',18,'LineWidth',3);
end

return

% -----------------------------------------------------------------------
% -----------------------------------------------------------------------

% Get unique colors for each object, but similar colors within a class

% -----------------------------------------------------------------------
% -----------------------------------------------------------------------

function colors = class_colors(classes,basecolors)

nclasses = length(unique(classes));
nobj = length(classes);

if nargin < 2 || isempty(basecolors), basecolors = {[1 0 0] [0 1 0] [0 0 1] [1 1 0] [0 1 1] [1 0 1]}; end

while nclasses > length(basecolors), basecolors = [basecolors basecolors]; end

colors = cell(1,nobj);
mmax = .4;       % max for off-color value
%omin = .8;       % min for on-color value
rands = .25;

for i = 1:max(classes)

    wh = find(classes == i);

    n = length(wh);

    x = sort(linspace(0,mmax,n),'descend')';

    c = repmat(x,1,3);

    % random variation
    c = c + rands .* (rand(size(c)) .* sign((rand(size(c)) > .5) - .5));
    c(c<0) = 0;
    c(c>1) = 1;

    % add basecolor back in
    whone = find(basecolors{i} > 0);
    c(:,whone) = 1;

    for j = 1:n
        colors(wh(j)) = {c(j,:)};
    end

end

%     for i = 1:nobj
% fprintf(1,'%3.0f\t%3.2f %3.2f %3.2f\n',MCA.ClusterSolution.classes(i),colors{i})
% end
return

% -----------------------------------------------------------------------
% -----------------------------------------------------------------------

% Utility functions

% -----------------------------------------------------------------------
% -----------------------------------------------------------------------

% This is now its own separate function
% % function f1 = create_figure(tagname)
% %     % checks for old figure, and creates new one
% %
% %     if nargin < 1 || isempty(tagname)
% %         tagname = 'nmdsfig';
% %     end
% %
% %     old = findobj('Tag', tagname);
% %
% %     if ~isempty(old)
% %         disp(['Clearing old image: ' tagname]);
% %         clf(old);
% %     end
% %
% %     scnsize = get(0,'ScreenSize');
% %
% %     xdim = min(scnsize(3)./2, 700);
% %     ydim = min(scnsize(4)./2, 700);
% %
% %     f1 = figure('position',round([50 50 xdim ydim]),'color','white');
% %     set(f1, 'Tag', tagname);
% %
% %     set(gca, 'FontSize', 16);
% %
% %
% %     return

function f1 = create_wide_figure
scnsize = get(0,'ScreenSize');

myposition = [50 50 scnsize(3)-100 scnsize(4)/2];
myposition(3) = min(myposition(3), 1200);
myposition(4) = min(myposition(4), 500);

f1 = figure('position',myposition,'color','white');
return

function save_figure(myname,npts)
if nargin < 2, npts = 400; end
disp(['Saving figure: ' myname '.png'])
f1 = gcf;
scn_export_papersetup(npts);
saveas(f1,myname,'png');
return

% -----------------------------------------------------------------------
% -----------------------------------------------------------------------

% Connect two points in 3-D space with a straight or curved line

% -----------------------------------------------------------------------
% -----------------------------------------------------------------------

function out = connect3d(varargin)

% ----------------------------------------------------------------------
% Parse inputs
% ----------------------------------------------------------------------

p = inputParser;

% Validation functions - customized for each type of input
% ----------------------------------------------------------------------

valfcn_scalar = @(x) validateattributes(x, {'numeric'}, {'nonempty', 'scalar'});

valfcn_number = @(x) validateattributes(x, {'numeric'}, {'nonempty'}); % scalar or vector

% Validation: Region object, atlas, structure, or [x1 x2 x3] triplet
valfcn_custom = @(x) isstruct(x) || isa(x, 'region') || isa(x, 'atlas') || (~isempty(x) && all(size(x) - [1 3] == 0) && all(isnumeric(x)));

% Validation: [x1 x2 x3] triplet
valfcn_xyz = @(x) validateattributes(x, {'numeric'}, {'nonempty', 'size', [1 3]});

% Required inputs
% ----------------------------------------------------------------------
p.addRequired('x', valfcn_custom);
p.addRequired('y', valfcn_custom);

% Optional inputs
% ----------------------------------------------------------------------
% Pattern: keyword, value, validation function handle

p.addParameter('color', [.9 .2 0], valfcn_xyz);
p.addParameter('bendpercent', .15, valfcn_number); % can be scalar or vector
p.addParameter('thickness', .1, valfcn_scalar);
p.addParameter('nstreamlines', 30, valfcn_scalar);
p.addParameter('streamlineshift', 10, valfcn_scalar);
p.addParameter('nsamples', [], valfcn_scalar);

% Parse inputs and distribute out to variable names in workspace
% ----------------------------------------------------------------------
% e.g., p.parse([30 1 0], [-40 0 10], 'bendpercent', .1);
p.parse(varargin{:});

IN = p.Results;
fn = fieldnames(IN);

for i = 1:length(fn)
    str = sprintf('%s = IN.(''%s'');', fn{i}, fn{i});
    eval(str)
end

% Adjust input vars as needed
% ----------------------------------------------------------------------

if isa(x, 'atlas'), x = atlas2region(x); end
if isa(y, 'atlas'), y = atlas2region(y); end

% x, y inputs can be [x,y,z] coordinate triplets or clusters.
if isstruct(x) || isa(x, 'region'), x = x.mm_center; end
if isstruct(y) || isa(y, 'region'), y = y.mm_center; end

% make x, y, z bend percents
if length(bendpercent) == 1
    bendpercent = repmat(bendpercent,1,3);
end

% make 3 coords, so we can bend if we want to
xcoords = [x(1); x(1)+(y(1)-x(1))./2 + bendpercent(1)*x(1); y(1)];
ycoords = [x(2); x(2)+(y(2)-x(2))./2 + bendpercent(2)*x(2); y(2)];
zcoords = [x(3); x(3)+(y(3)-x(3))./2 + bendpercent(3)*x(3); y(3)];

n = length(xcoords);

if isempty(nsamples)
    nsamples = 10 * n;
end

% Plot streamlines
% ----------------------------------------------------------------------

for i = 1:nstreamlines % number of streamlines

    randval = unifrnd(-.15, .15, 4, 1);

    mycolor = color + randval(1);
    mycolor(mycolor > 1) = 1;
    mycolor(mycolor < 0) = 0;

    if any(bendpercent)
        % bow out: curved line

        myx = xcoords;
        myy = ycoords;
        myz = zcoords;

        % random bend
        myx(2) = myx(2) + streamlineshift*randval(2);
        myy(2) = myy(2) + streamlineshift*randval(2);
        myz(2) = myz(2) + streamlineshift*randval(2);

        t = 1:n;
        ts = 1:((n-1)/(nsamples-1)):n;          % spline grid

        myx = spline(t, myx, ts);
        myy = spline(t ,myy, ts);
        myz = spline(t, myz, ts);
    end

    % random shift
    myx = myx + streamlineshift * randval(2);
    myy = myy + streamlineshift * randval(3);
    myz = myz + streamlineshift * randval(4);

    h(i) = plot3(myx, myy, myz, 'Color', mycolor,'LineWidth', thickness);

end

% h = plot3(xcoords,ycoords,zcoords,'Color',color,'LineWidth',thickness);

out.h = h;
out.xcoords = xcoords;
out.ycoords = ycoords;
out.zcoords = zcoords;

return

% -----------------------------------------------------------------------
% -----------------------------------------------------------------------

% Connect two points in 2-D space with a straight or curved line

% -----------------------------------------------------------------------
% -----------------------------------------------------------------------
function out = connect2d(x,y,color,thickness,bendpercent,varargin)

% x, y inputs can be [x,y,z] coordinate triplets or clusters.
if isstruct(x), x = x.mm_center(1:2); end
if isstruct(y), y = y.mm_center(1:2); end

% make x, y bend percents
if length(bendpercent) == 1
    bendpercent = repmat(bendpercent,1,2);
end

% make 2 coords, so we can bend if we want to
xcoords = [x(1); x(1)+(y(1)-x(1))./2 + bendpercent(1)*x(1); y(1)];
ycoords = [x(2); x(2)+(y(2)-x(2))./2 + bendpercent(2)*x(2); y(2)];

if any(bendpercent)
    % bow out: curved line
    n = length(xcoords);

    nsamples = [];
    if length(varargin) > 0
        nsamples = varargin{1};
    end

    if isempty(nsamples)
        nsamples = 10 * n;
    end

    t = 1:n;
    ts = 1:((n-1)/(nsamples-1)):n;          % spline grid

    xcoords = spline(t,xcoords,ts);
    ycoords = spline(t,ycoords,ts);

end

h = plot(xcoords,ycoords,'Color',color,'LineWidth',thickness);

out.h = h;
out.xcoords = xcoords;
out.ycoords = ycoords;

return


text:
classdef SDUpdater < handle
    properties (Hidden)
        handles;            %Graphics objects
        BoardRDTable;       % Table of boards / reference designs
        VideoPSPDir;        % Root PSP directory
        SDList;             % List of SD objects
        BoardVariantTable = containers.Map();  % Table of boards / variants / zip files
    end

    properties (Hidden, Constant)
        FWMODE_ZIP = 'FAT32 ZIP';
        FWMODE_IMG = 'Disk Image';
        DEVICETYPE_ZYNQ = 'Zynq';
        DEVICETYPE_ALTERASOC = 'Altera SoC';
    end

    properties
        AppName = 'SD Updater'
        FWMode@hdlbsp.util.SDUpdater.FWModes = hdlbsp.util.SDUpdater.FWModes.FAT32_ZIP;
        DeviceType@hdlbsp.util.SDUpdater.DeviceTypes = hdlbsp.util.SDUpdater.DeviceTypes.ZYNQ;
        GUIEnabled@logical = true;
    end

    properties (Access = protected)
        FWUpdater
        privActiveVariant
        privActiveBoard
        privActiveSDDrive = 'N/A'
    end

    properties (Dependent)
        ActiveSDDrive
        ActiveZip
        ActiveVariant
        ActiveBoard
    end

    properties (Hidden, Dependent)
        WipeSD         % Preserve SD card contents on udpate
        ActiveBoardObj
        ActiveVariantList
        DriveList
    end

    methods (Hidden)
        function app = SDUpdater(varargin)
            % Configure the app
            app.configureApp;

            % Load the variants
            app.loadVariants;

            % Launch the app
            app.launchApp;
        end

        function delete(app)
            try %#ok
                delete(app.handles.Fig);
            end
        end

        function createGUI(app)
            LeftTextAlgin = 10;
            LeftMenuAlign = 150;
            MenuHeight = 27;
            LineSize = 40;
            CurrentVertPos = LineSize*7+10;
            CheckboxWidth = 50;

            % Create the GUI
            app.handles.Fig = figure(...
                'Color',[0.94 0.94 0.94],...
                'Position',[532 624 500 CurrentVertPos],...
                'MenuBar','none',...
                'Name',app.AppName,...
                'NumberTitle','off',...
                'Resize','off',...
                'Visible','on',...
                'CloseRequestFcn', @app.callback_closeApp);

            % Create the SD card list
            CurrentVertPos = CurrentVertPos - LineSize;

            app.handles.menuSDCard = app.uimenu(...
                'N/A', [LeftMenuAlign CurrentVertPos 61 MenuHeight],...
                @app.callback_menuSDCard);

            app.handles.txtSDCard = app.uitext(...
                'Drive', [LeftTextAlgin CurrentVertPos]);

            app.handles.btnRefresh = uicontrol(...
                'Parent',app.handles.Fig,...
                'Callback',@app.callback_updateSDList,...
                'Position',[69 CurrentVertPos 60 MenuHeight],...
                'String','Refresh');

            % Create the write button
            app.handles.btnWrite = uicontrol(...
                'Parent',app.handles.Fig,...
                'Callback',@app.callback_WriteSD,...
                'Position',[390 CurrentVertPos 101 MenuHeight],...
                'String','Write SD');

            % Create the board list
            CurrentVertPos = CurrentVertPos - LineSize;

            app.handles.menuBoard = app.uimenu(...
                'N/A', [LeftMenuAlign CurrentVertPos 337 MenuHeight],...
                @app.callback_menuBoard);

            app.handles.txtBrd = app.uitext(...
                'Board', [LeftTextAlgin CurrentVertPos]);

            % Create the reference design list
            CurrentVertPos = CurrentVertPos - LineSize;

            app.handles.menuVariant = app.uimenu(...
                'N/A', [LeftMenuAlign CurrentVertPos 337 MenuHeight],...
                 @app.callback_menuVariant);

            app.handles.txtVariant = app.uitext(...
                'Variant', [LeftTextAlgin CurrentVertPos]);

            if isequal(app.FWMode, hdlbsp.util.SDUpdater.FWModes.FAT32_ZIP)
                % Create the preserve SD checkbox
                CurrentVertPos = CurrentVertPos - LineSize;

                app.handles.txtWipeSD = app.uitext(...
                    'Wipe SD Card', [LeftTextAlgin CurrentVertPos]);

                app.handles.checkWipeSD = uicontrol(...
                    'Parent',app.handles.Fig,...
                    'Position',[LeftMenuAlign CurrentVertPos CheckboxWidth MenuHeight],...
                    'Callback',@app.callback_WipeSD,...
                    'Style','checkbox',...
                    'Value',app.WipeSD);
            end

        end

        function handle = uitext(app, Text, Position)
            % Set the width/height to 1
            Position(3:4) = 1;

            % Create the text control
            handle = uicontrol(...
                'Parent',app.handles.Fig,...
                'Position',Position,...
                'HorizontalAlignment', 'Left',...
                'String',Text,...
                'Style','text');

            % Resize based on the string
            ext = get(handle, 'Extent');
            Position(3:4) = (ext(3:4))+3;
            set(handle, 'Position', Position);

        end

        function handle = uimenu(app, Text, Position, Callback)
            if nargin < 4
                Callback = {};
            end
            handle = uicontrol(...
                'Parent',app.handles.Fig,...
                'BackgroundColor', [1 1 1 ],...
                'Position',Position,...
                'String',Text,...
                'BusyAction', 'cancel',...
                'Callback',Callback,...
                'Style','popupmenu',...
                'Value',1);
        end

        function handle = uiedit(app, Text, Position, Callback)
            if nargin < 4
                Callback = {};
            end
            handle = uicontrol(...
                'Parent',app.handles.Fig,...
                'BackgroundColor', [1 1 1 ],...
                'Position',Position,...
                'String',Text,...
                'BusyAction', 'cancel',...
                'Callback',Callback,...
                'Style','edit',...
                'Value',1);
        end

    end
    %% GUI Callbacks
    % Syntax: method(app,hObject,eventdata)
    methods (Hidden)

        function callback_closeApp(app,varargin)
            delete(app.handles.Fig);
            delete(app);
        end

        function callback_WriteSD(app,varargin)
            app.enableUI(false);
            cleanup = onCleanup(@()app.enableUI(true));

            if isequal(app.getSelectedItem(app.handles.menuSDCard), 'N/A')
                error('You must select a valid drive to program the SD card');
            end

            if app.WipeSD
                warnVal = 'wipe all';
            else
                warnVal = 'overwrite';
            end
            WarnStr = sprintf(['WARNING: Updating the SD Card will %s existing contents\n'...
                    'Are you sure you want to proceed?'], warnVal);

            ButtonName = questdlg(WarnStr, ...
                                     'Update SD Card', ...
                                     'No');
            if isequal(ButtonName, 'Yes')
                app.writeSDCard();
            end
        end

        function callback_updateSDList(app,varargin)
            app.enableUI(false);
            cleanup = onCleanup(@()app.enableUI(true));

            driveList = app.updateSDList;
            set(app.handles.menuSDCard, 'String', driveList)
        end

        function callback_menuSDCard(app,varargin)
            app.privActiveSDDrive = app.getSelectedItem(app.handles.menuSDCard);
        end

        function callback_menuBoard(app,varargin)
            app.ActiveBoard = app.getSelectedItem(app.handles.menuBoard);
            set(app.handles.menuVariant, 'String', app.ActiveVariantList);
        end

        function callback_menuVariant(app,varargin)
            app.ActiveVariant = app.getSelectedItem(app.handles.menuVariant);
        end

        function callback_WipeSD(app,varargin)
            app.WipeSD = get(app.handles.checkWipeSD, 'Value');
        end
    end
    %% Set / Get Properties
    methods
        function val = get.WipeSD(app)
            if isequal(app.FWMode, hdlbsp.util.SDUpdater.FWModes.FAT32_ZIP)
                val = getpref('HDLBSP_COMMON', 'SDUpdater_WipeSD', 0);
            else
                val = 1;
            end
        end

        function set.WipeSD(app, val) %#ok
            if isequal(val, 0) || isequal(val, 1)
                setpref('HDLBSP_COMMON', 'SDUpdater_WipeSD', val);
            else
                error('Must specify 0 or 1 for WipeSD');
            end
        end

        function variantList = get.ActiveVariantList(app)
            boardObj = app.ActiveBoardObj();
            variantList = keys(boardObj.Variants);
        end

        function brdObj = get.ActiveBoardObj(app)
            brdObj = app.BoardVariantTable(app.ActiveBoard);
        end

        function set.ActiveVariant(app, varName)
            if ~ismember(varName, app.ActiveVariantList)
                error('%s is not an available variant.', varName);
            end
            app.privActiveVariant = varName;
        end

        function varName = get.ActiveVariant(app)
            varName = app.privActiveVariant;
        end

        function set.ActiveBoard(app, boardName)
            if ~ismember(boardName, keys(app.BoardVariantTable))
                error('%s is not an available board.', boardName);
            end
            app.privActiveBoard = boardName;
        end

        function boardName = get.ActiveBoard(app)
            boardName = app.privActiveBoard;
        end

        function zipFile = get.ActiveZip(app)
            brdObj = app.ActiveBoardObj;
            varName = app.ActiveVariant;
            zipFile = brdObj.Variants(varName);
        end

        function driveList = get.DriveList(app)
            driveList = app.FWUpdater.getDriveList;
            if isempty(driveList)
                driveList = {'N/A'};
            end
        end

        function set.ActiveSDDrive(app, sdDrive)
            if ~ismember(sdDrive, app.DriveList)
                error('%s is not an available SD card', sdDrive);
            end
            app.privActiveSDDrive = sdDrive;
        end

        function sdDrive = get.ActiveSDDrive(app)
            sdDrive = app.privActiveSDDrive;
        end

    end

    %% API Functions
    methods (Abstract, Access = protected)
        % Configure the App settings
        configureApp(app)
        % Load the SD card variants
        loadVariants(app)
    end

    methods
        % Update the SD card contents
        function writeSDCard(app)
            app.FWUpdater.Drive = app.ActiveSDDrive;

            switch(app.FWMode)
                case hdlbsp.util.SDUpdater.FWModes.FAT32_ZIP
                    app.writeSDCard_FAT32;
                case hdlbsp.util.SDUpdater.FWModes.DISK_IMAGE
                    app.writeSDCard_IMG;
                otherwise
                    error('Invalid setting for FWMode');
            end
        end
    end

    methods (Access = protected)
        function addVariant(app, BoardName, VariantName, ZipFile)
            if app.BoardVariantTable.isKey(BoardName)
                tblObj = app.BoardVariantTable(BoardName);
            else
                tblObj.BoardName = BoardName;
                tblObj.Variants = containers.Map();
            end
            tblObj.Variants(VariantName) = ZipFile;
            % if it's the first one, set it to active
            if isempty(app.privActiveBoard)
                app.privActiveBoard = BoardName;
                app.privActiveVariant = VariantName;
            end
            app.BoardVariantTable(BoardName) = tblObj;
        end

        function launchApp(app)
            % Load some parameters
            if isempty(app.FWUpdater)
                switch(app.DeviceType)
                    case hdlbsp.util.SDUpdater.DeviceTypes.ZYNQ
                        app.FWUpdater = zynq.setup.ZynqFirmwareUpdate;
                    case hdlbsp.util.SDUpdater.DeviceTypes.ALTERA_SOC
                        app.FWUpdater = codertarget.alterasoc.setup.AlteraSoCFirmwareUpdate;
                end
            end

            % Update the SD Card List
            driveList = app.updateSDList;

            % Create the GUI
            if app.GUIEnabled
                app.createGUI();
                % Update the list of SD cards
                set(app.handles.menuSDCard, 'String', driveList)

                % Setup the board list
                set(app.handles.menuBoard, 'String', app.BoardVariantTable.keys);
                app.setSelectedItem(app.handles.menuBoard, app.ActiveBoard);

                % Setup the variant list
                set(app.handles.menuVariant, 'String', app.ActiveVariantList);
                app.setSelectedItem(app.handles.menuVariant, app.ActiveVariant);
            end
        end

        % Update a FAT32 SD Card
        function writeSDCard_FAT32(app)
            hWaitbar = waitbar(0,['SD Card Update...' ...
                char(ones(1,30, 'like', uint8(0)).*uint8(' '))]);
            cleanup = onCleanup(@()delete(hWaitbar));

            if app.WipeSD
                waitbar(0.25,hWaitbar,'Formatting SD Card...');
                app.FWUpdater.formatFATWindows;
            end

            % Program the SD Card
            waitbar(0.50,hWaitbar,'Updating SD Card Contents...');
            unzip(app.ActiveZip, [app.ActiveSDDrive filesep]);
            waitbar(0.75,hWaitbar,'Applying final updates...');
        end

        % Update a image SD Card
        function writeSDCard_IMG(app)
            fwDir = tempname;
            mkdir(fwDir);
            cleanup = onCleanup(@()rmdir(fwDir, 's'));
            % Get the info on the image file
            [~, fwFileName, ~] = fileparts(app.ActiveZip);

            % populate the FW Updater
            app.FWUpdater.DownloadFolder = fwDir;
            app.FWUpdater.FirmwareName = fwFileName;
            gunzip(app.ActiveZip, fwDir);

            % Determine the file size
            fileInfo = dir(fullfile(fwDir, fwFileName));
            app.FWUpdater.FirmwareSize = fileInfo.bytes;
            % Write the image
            app.FWUpdater.writeImage;
        end
    end

    %% Helper Functions
    methods (Access = protected)
        function driveList = updateSDList(app)
            driveList = app.DriveList;
            if ~ismember(app.ActiveSDDrive, driveList)
                app.privActiveSDDrive = driveList{1};
            end
        end
        function enableUI(app, en)
            elements = fields(app.handles);
            for ii = 1:numel(elements)
                handle = app.handles.(elements{ii});
                if ~isfield(get(handle), 'Enable')
                    continue;
                end
                if en
                    set(handle, 'Enable', 'on');
                else
                    set(handle, 'Enable', 'off');
                end
            end
        end
    end

    %% Static Methods
    methods (Hidden, Static)
        function item = getSelectedItem(menu)
            menuList = get(menu, 'String');
            item = menuList{get(menu, 'Value')};
        end
        function setSelectedItem(menu, item)
            menuList = get(menu, 'String');
            [is,Idx] = ismember(item, menuList);
            if ~is
                error('Could not find %s in list %s', item, menuList);
            end
            set(menu, 'Value', Idx);
        end
    end
end

author: analogdevicesinc


text:
function [vDistances_] = ex_MapDist3D(mCatalog, mPolygon, vPolygon)
% function [caNodeIndices] = ex_CreateIndexCatalog(mCatalog, mPolygon, bMap, nGriddingMode,
%                                                  nNumberEvents, fRadius, fSizeRectHorizontal, fSizeRectDepth)
% -------------------------------------------------------------------------------------------------------------
% Creates a cell-array with subcatalogs for every grid node defined by mPolygon. These subcatalogs
%   contain only indices to the earthquake "rows" in mCatalog.
%
% Input parameters:
%   mCatalog              Earthquake catalog
%   mPolygon              Polygon (defined by ex_selectgrid)
%   vPolygonDepth         Polygon depth
%
% Output parameters:
%   caNodeIndices         Cell-array with index-catalogs per grid node of mPolygon
%
% Thomas van Stiphout
% July 28, 2008

global bDebug;
if bDebug
  report_this_filefun(mfilename('fullpath'));
end

% Create the catalogs for each node with pointers to the overall catalog
nNumberNodes_ = length(mPolygon(:,1));

% vResolution give the radius (for nGriddingMode = 0) and no. of events
% (for nGriddingMode = 1)
vResolution_(:,1)=ones(nNumberNodes_,1)*nan;

% Loop over all points of the polygon
for nNode_ = 1:nNumberNodes_
  % Get the grid node coordinates
  fX_ = mPolygon(nNode_, 1);
  fY_ = mPolygon(nNode_, 2);
  fZ_ = vPolygon(nNode_, 1);

  % Calculate distance from center point
  vDistances_ = sqrt(((mCatalog(:,1)-fX_)*cos(pi/180*fY_)*111).^2 + ((mCatalog(:,2)-fY_)*111).^2 + (mCatalog(:,7)-fZ_).^2);

  % Fixed radius
  % Use all events within fRadius
%       caNodeIndices{nNode_} = find(vDistances_ <= fRadius);
%       vResolution_(nNode_) = length(find(vDistances_ <= fRadius));
end; % of for nNode_

author: CelsoReyes


text:
function [equalities,redundant] = mpt_detect_fixed_rows(Matrices)

A = [ Matrices.G -Matrices.E];
b = Matrices.W;
lb = Matrices.lb;
ub = Matrices.ub;

fixed = find(lb == ub);
if ~isempty(fixed)
    b = b-A(:,fixed)*lb(fixed);
    A(:,fixed) = [];
    ub(fixed)=[];
    lb(fixed)=[];
end

[AA, BB, AAeq, BBeq,indeq] = mpt_ineq2eq(A, b);
if ~isempty(indeq)
    i = indeq(:,1);
    j = indeq(:,2);
    equalities = unique(i);
    redundant = unique([i;j]);
    equalities = setdiff(equalities,intersect(i,j));  % The case [43 51;51 89;89 97]
    AL0A  = (A>0).*A;
    AG0A  = (A<0).*A;
    bi_up = AL0A*ub+AG0A*lb;
    redundant = unique([redundant;find([bi_up<b])]);
else
    AL0A  = (A>0).*A;
    AG0A  = (A<0).*A;
    bi_up = AL0A*ub+AG0A*lb;
    redundant = unique(find(bi_up<b));
    equalities = [];
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function z = dClambdadx(lambda)
z = log2(lambda)./(1 + lambda).^2;
end

author: YuYongRun

repo: PolarCodeDecodersInMatlab