text:
function S = xval_SVM(varargin)
%  SVM for repeated-measures (within-person) classification, with repeated cross-val and nested cross-val options
%
% :Usage:
% ::
%
% S = xval_SVM(X, Y, id, varargin)
%
% Steps and features:
% -------------------------------------------------------------------------
% - Select holdout sets, keeping images from the same id together and stratifying on the outcome to be predicted
% - Fit overall model and get "sanity check" accuracy
% - Cross-validation with standard a priori hyperparameter choices
% - Plot scores and ROC curves
% - Nested cross-val with hyperparameter optimization [optional]
%   * Estimate model performance accounting for search across hyperparams
%   * Updates S.yfit, S.dist_from_hyperplane_xval, S.class_probability_xval
%   * Updates crossval_accuracy, Y_within_id, scores_within_id, scorediff, crossval_accuracy_within, classification_d
%   * Retrain to obtain estimate of best hyperparameters for future testing (update S.modeloptions)
% - Repeat cross-validation of optimized model with different random splits [optional]
% - Fit model on full dataset with final hyperparameters to get predictor weights (betas, b)
% - Bootstrap model parameter estimates (w) and get P-values, FDR-corrected significant features
% - Print output and text compatible with report-generation
%
% Notes:
% -------------------------------------------------------------------------
% - Uses Matlab's Stats/ML Toolbox SVM object, fitcsvm, and hyperparameter
% optimization. Other options, e.g., fitclinear, may be better for some
% situations (large num. variables)
%
% - Single-interval accuracy from cross-val and ROC_plot may differ because
% ROC_plot chooses a new score threshold that maximizes overall balanced
% accuracy. Forced-choice accuracy should be identical.
%
% Optimizing hyperparameters:
% HyperparameterOptimizationOptions include defaults of:
% - 5-fold CV (not grouped by id) without repartitioning
% - Bayesian optimization, using best estimates from smooth function
% - maximizes accuracy, so may need large(ish) samples to be effective.
%
% - works on linear classifiers only, but can explore nonlinear models
% with a small change in the code (now commented out). This could be added
% as an optional flag as well.
%
% - If optimizing hyperparameters AND repeating cross-validation, accuracy
% estimates will use nested cross-validation, and can take a long time to run
%
%
% Class probability estimates:
% - crossval returns cross-validated yfit (class predictions), scores (dist_from_hyperplane_xval),
% class_probability_xval. Scores and class probabilties will diverge (may
% not be perfectly correlated) because the sigmoid scaling (Platt scaling) varies across folds.
% They will diverge more if there is no true signal.
%
% ..
%     Author and copyright information:
%
%     Copyright (C) 2020 Tor Wager
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
% ..
%
% :Inputs:
%
%   **X:**
%        obs x variables numeric matrix of predictors
%
%   **Y:**
%        obs x 1 numeric vector of outcomes, effects-coded (1, -1)
%
%   **id:**
%        obs x 1 numeric vector of grouping codes, e.g., for participants
%        all obs from each group will be included together in a training or test set
%        if no grouping, use integers 1...n images or leave empty
%
% :Optional Inputs:
%   **'doplot', [logical flag]:**
%        Create plots; default = true. 'noplot' to turn off.
%
%   **'doverbose', [logical flag]:**
%        Verbose output; default = true. 'noverbose' to turn off.
%
%   **'dooptimize', [logical flag]:**
%        Optimize hyperparameters; default = true. 'nooptimize' to turn off.
%
%   **'doprepeats', [num repeats]:**
%        Repeat cross-val with different partitions; default = 10.
%        Enter number of repeats, 'norepeats' to turn off.
%
%   **'modeloptions', [modeloptions cell]:**
%        Options for model structure and hyperparameters
%        Cell array of keyword-value pairs, as specified in fitcsvm (Matlab Stats/ML toolbox)
%        Default: {'KernelFunction', 'linear'}
%
%   **'nfolds', [num folds]:**
%        Set number of cross-validation folds; default = 10.
%        Enter number of folds.
%
% :Outputs:
%
%   **S:**
%        Structure with model output
%        of them (partially consistent with this function).
%                           Y: Actual (obs) outcome - should be 1, -1 for SVM analysis
%                        yfit: Predicted outcome - should be 1, -1 for SVM analysis
%                              Cross-validated, so can be used as series
%                              of predicted values based on brain
%                              measures (very useful!)
%                          id: Grouping variable for within-participant observations
%                      accfun: Function handle to get accuracy (single-interval)
%                           w: Model weights/betas, weights on input variables
%                      nfolds: Number of folds in holdout set
%                 cvpartition: Object with information about training/test sets
%                       teIdx: Testing IDs for each fold (holdout set)
%                       trIdx: Training IDs for each fold (holdout set)
%             dist_from_hyperplane_xval: Cross-validated distance perpendicular to class boundary
%                               higher = stronger prediction in favor of Class 1, lower = in favor of class -1.
%                               Useful! use as continuous measure of
%                               observation scores. Can calculate
%                               effect sizes from this, for example.
%       class_probability_xval: Cross-validated distance expressed as probabilities of Class 1, using Platt scaling
%            crossval_accuracy: Cross-validated accuracy (input options; no hyperparam opt)
%  classification_d_singleinterval: Classification effect size (input options; no hyperparam opt)
%crossval_accuracy_opt_hyperparams: Cross-validated accuracy with optimized hyper-parameters
%                  Y_within_id: Outcomes arranged by id, for within-person comparisons
%             scores_within_id: SVM scores arranged by id, for within-person comparisons
%                    scorediff: Within-person SVM scores arranged by id, for within-person comparisons
%     crossval_accuracy_within: Within-person cross-validated accuracy (input options; no hyperparam opt)
%      classification_d_within: Within-person classification effect size (input options; no hyperparam opt)
%                 S.boot_w_ste: Bootstrapped standard errors of model weights (feature-level)
%                S.boot_w_mean: Bootstrapped mean model weights (feature-level)
%                         S.wZ: Bootstrapped Z-scores of model weights (feature-level)
%                         S.wP: Bootstrapped P-values for individual model weights (feature-level)
%                 S.wP_fdr_thr: P-value threshold for FDR q < 0.05
%              S.boot_w_fdrsig: Logical vector for which weights are significant with FDR
%               S.w_thresh_fdr: Thresholded weights at FDR q < 0.05 based on bootstrapping (feature-level)
%                   S.SVMModel: ClasssificationSVM model object trained on full dataset, final chosen parameter set; for predicting in subsequent validation samples.
%                               Y-hat = (X/s)'* w + b
%                               ClassificationSVM objects store w, b, and s in the properties Beta, Bias, and KernelParameters.Scale, respectively.
%
% :Examples:
% ::
%
% Simulate some data with true signal, with two observations per person
% -------------------------------------------------------------------------
% n = 50; % participants
% k = 120; % features
% true_sig = [repmat(randn(1, k), n, 1); repmat(randn(1, k), n, 1)]; % First n are Class 1, second n are Class 2
% noise = 10 * randn(2 * n, k);         % Noise var >> true signal var
% X = true_sig + noise;                  % Predictors
% Y = [ones(n, 1); -ones(n, 1)];         % Outcome to classify, coded [1 -1]
% id = [(1:n)'; (1:n)'];                 % Grouping ID codes for participants
%
% S = xval_SVM(X, Y, id, 'nooptimize', 'norepeats', 'nobootstrap');    % Fastest for quick cross-validated performance
% S = xval_SVM(X, Y, id, 'nooptimize', 'norepeats', 'nboot', 100);     % Very quick test of bootstrapping with few samples (not for final models)
% S = xval_SVM(X, Y, id, 'nooptimize', 'dorepeats', 5, 'nobootstrap'); % Quick test of repeated cross-val, 5x
% S = xval_SVM(X, Y, id, 'nooptimize', 'nobootstrap');                 % Repeat cross-val only
% S = xval_SVM(X, Y, id, 'norepeats', 'nobootstrap');                  % Optimize with nested cross-val only
% S = xval_SVM(X, Y, id);                                              % Optimize, repeat with optimal model, bootstrap
%
% S = xval_SVM(X, Y, id, 'nooptimize', 'norepeats', 'nobootstrap', 'noverbose', 'noplot');  % Returns only the output S
%
% Train the same model without within-person observations (i.e., 1 observation per id):
% This does single-interval classification only
%
% S = xval_SVM(X, Y, (1:length(Y))', 'norepeats', 'nobootstrap');
%
% :References:
%   See Mathworks functions
%
% :See also:
%   - fmri_data.predict, canlab_run_paired_SVM, other xval_ functions
%

% ..
%    Programmers' notes:
%    Created by Tor Wager, May 2020
%    ver 1.0 (no version number listed) - some bugs in optimization
%             Enforce [1 -1] inputs
%    ver 1.1 (Tor): fixed bugs, optimization changes and documentation, multiple small tweaks and some new plots
%             Updated to handle within-person or between-person obs only
%             ROC_plot should be at threshold 0 for single-interval
%             Get rid of plots for all folds in opt output
%             OptimizeHypParams returns point est. No need to rerun. Builds in 5-fold eval by default
%             Add ?Summary of nested cross-val accuracy with hyperparameter optimization?
%             Optim.: Allow control of verbosity and plots. No plots/verbose in nested eval.
%             Optim: use best estim, not best observed
%             Have a look at the data and add data plots. Sorted by cases
%   ver 1.2   Cosmetic documentation and look updates
%
%   Future:
%
%             Troubleshoot prediction vs. class probability reversals
%                   % either scores or pscores appear to be inconsistent in terms of which column is class -1 and which is class 1. Check for instability in scores, and select pscores to be consistent with this on each fold.
%             Add permutation test
%     v1.3
%     Linear/nonlinear - pass in hyparams to optimize.
%     Return model weights from folds for further analysis/application
%
%     v1.4
%     Add nfolds argument - Michael Sun, 07/07/2022
%
% ..

%% ----------------------------------------------------------------------
% Version
% ----------------------------------------------------------------------

ver = 1.2;

%% ----------------------------------------------------------------------
% Parse inputs
% ----------------------------------------------------------------------
% Uses the inputParser object. Older schemes are below.

% Logical flags - parse specially because they are not name-value pairs
% Parser accepts only name/value pairs, so this code allows a workaround
% We must remove these inputs because they will cause inputparser to error
% ----------------------------------------------------------------------
wh = strcmp(varargin, 'noverbose'); if any(wh), doverbose1 = false; varargin(wh) = []; end
wh = strcmp(varargin, 'noplot'); if any(wh), doplot1 = false; varargin(wh) = []; end
wh = strcmp(varargin, 'noplots'); if any(wh), doplot1 = false; varargin(wh) = []; end
wh = strcmp(varargin, 'nooptimize'); if any(wh), dooptimize1 = false; varargin(wh) = []; end
wh = strcmp(varargin, 'norepeats'); if any(wh), dorepeats1 = 0; varargin(wh) = []; end
wh = strcmp(varargin, 'nobootstrap'); if any(wh), dobootstrap1 = 0; varargin(wh) = []; end

% The handling of parsing keyword-value pairs, checking attributes, and
% specifying default values is handled here, using Matlab's inputparser object:
ARGS = parse_inputs(varargin{:});

fn = fieldnames(ARGS);

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

% Replace logical flags because they will be overwritten by the parser
if exist('doverbose1', 'var'), doverbose = doverbose1; end
if exist('doplot1', 'var'), doplot = doplot1; end
if exist('dooptimize1', 'var'), dooptimize = dooptimize1; end
if exist('dorepeats1', 'var'), dorepeats = dorepeats1; end
if exist('dobootstrap1', 'var'), dobootstrap = dobootstrap1; end
if ~exist('nfolds', 'var'), nfolds=10; end

if isempty(id), id = [1:length(Y)]'; end
if ~iscolumn(id), id = id'; end

%% Select holdout sets for outer loop
% - Keep images from the same id together
% - Stratify on the outcome to be predicted
% -------------------------------------------------------------------------

S = struct();   % Define structure for models and output; use variable names in fmri_data.predict when possible

S.Y = Y;
S.id = id;      % Subject grouping - not in fmri_data.predict

S.modeloptions = modeloptions;

S.accfun = @(Y, yfit) 100 .* nansum(Y == yfit) ./ sum(~isnan(Y));

if doverbose

    dashes = '---------------------------------------------------------------------------------';
    printhdr = @(str) fprintf('%s\n%s\n%s\n', dashes, str, dashes);

    printhdr(' ') %#ok<*UNRCH>
    printhdr(sprintf('xval_SVM ver %3.2f: Cross-validated Support Vector Machine classification', ver));
    printhdr(' ')

    disp(' ')
    disp(' ')
    printhdr('Selecting and analyzing cross-validation folds');

end

% Preliminary data plots
% -------------------------------------------------------------------------
if doplot
    prelim_data_plot(X, Y, id);
end

% Select train/test sets
% -------------------------------------------------------------------------

% [S.trIdx, S.teIdx] = xval_stratified_holdout_leave_whole_subject_out(S.Y, S.id, 'doverbose', doverbose, 'doplot', doplot);
[S.trIdx, S.teIdx] = xval_stratified_holdout_leave_whole_subject_out(S.Y, S.id, 'doverbose', doverbose, 'doplot', doplot, 'nfolds', nfolds);  % MS 6/16/2022: Crashes if you can't do 10-fold.
drawnow, snapnow

% Fit the overall a priori model: SVM with linear kernel
% -------------------------------------------------------------------------
if doverbose

    printhdr('Model training and cross-validation without optimization')
    fprintf('Training overall model. ')

    SVMModel = fitcsvm(X, S.Y, S.modeloptions{:}); %#ok<*NODEF>

    % non-cross-val accuracy: sanity check. should usually be 100% if training is overfit and n >> p
    label = predict(SVMModel, X);

    acc = S.accfun(S.Y, label);
    fprintf('Accuracy without cross-val: %3.0f%%\n', acc);

end

%% Cross-validation with standard a priori hyperparameter choices
% -------------------------------------------------------------------------
S.nfolds = length(S.trIdx);

% crossval returns cross-validated yfit (class predictions), scores (dist_from_hyperplane_xval),
% class_probability_xval. Scores and class probabilties will diverge (may
% not be perfectly correlated) because the sigmoid scaling (Platt scaling) varies across folds.
% They will diverge more if there is no true signal.

S = crossval(S, X, doverbose);

% Summarize accuracy
% -------------------------------------------------------------------------
% Given Y, yfit, cross-val scores (dist_from_hyperplane_xval), id
% Update crossval_accuracy,classification_d_singleinterval
% Y_within_id, scores_within_id, scorediff,
% crossval_accuracy_within, classification_d_within
S = update_accuracy_stats(S);

if doverbose

    fprintf('Accuracy (cross-validated): Single-interval: %3.0f%%, d = %3.2f\n', S.crossval_accuracy, S.classification_d_singleinterval);

    if S.mult_obs_within_person

        fprintf('                            Forced-choice within-person: %3.0f%% d = %3.2f\n', S.crossval_accuracy_within, S.classification_d_within);

    end

end

% Plot scores and ROC curves
% -------------------------------------------------------------------------
if doplot

    plot_scores_and_ROC(S);

    xval_plot_scores_vs_class_probability(S);

end

%% Nested cross-val with hyperparameter optimization
% -------------------------------------------------------------------------

if dooptimize

    if doverbose
        printhdr('Optimizing hyperparameters using nested cross-validation');
    end

    % Nested cross-validation (updates S.yfit, S.dist_from_hyperplane_xval,
    % S.class_probability_xval)
    % Estimate model performance accounting for search across hyperparams
    % ---------------------------------------------------------------------
    % Updates S.yfit, S.dist_from_hyperplane_xval, S.class_probability_xval
    S = crossval_nested(S, X, doverbose);

    % Given Y, yfit, cross-val scores (dist_from_hyperplane_xval), id
    % Update crossval_accuracy, Y_within_id, scores_within_id, scorediff, crossval_accuracy_within, classification_d
    S = update_accuracy_stats(S);

    % Retrain to obtain estimate of best hyperparameters for future testing (update S.modeloptions)
    % ---------------------------------------------------------------------
    % Modal parameters across folds is not
    % the best way to get the final hyperparameters. Instead, the model should
    % be re-fit using single cross-validation with all hparam options pick the
    % best relative model. Nested cross-val will estimate the accuracy of the
    % whole procedure, including the hparam search, but not return the final best
    % parameters.
    hparams_to_optimize = 'all'; % {'BoxConstraint' 'Standardize'}
    Mdl = fitcsvm(X, S.Y,...
        'OptimizeHyperparameters', hparams_to_optimize, 'HyperparameterOptimizationOptions',...
        struct('AcquisitionFunctionName','expected-improvement-plus', 'Verbose', double(doverbose), 'ShowPlots', double(doplot))); %#ok<*IDISVAR>

    % Notes:
    % - Extract hyperparameter choices in a format that we can pass back in
    % to train or test a new model.
    % - We could use the best observed point, but it's preferred to use the
    % best estimated values from a smooth model fit to observed samples.

    best_modeloptions = convert_hyperparameter_choices_to_cell(Mdl, []);                % optimizing "all"
    % best_modeloptions = convert_hyperparameter_choices_to_cell(Mdl, S.modeloptions);   % optimizing "linear" : ms 06242022 Changed default settings to not optimize

    S.modeloptions = best_modeloptions;

    % Summarize accuracy
    % -------------------------------------------------------------------------

    if doverbose

        disp('Summary of nested cross-val accuracy with hyperparameter optimization')
        disp(' ')

        disp('Optimized hyperparameters for each fold')
        disp(S.hyperparams_by_fold)

        disp(' ')
        disp('Best options, stored in S.modeloptions:');

        disp(best_modeloptions)

        fprintf('Accuracy (nested cross-validation): Single-interval: %3.0f%%, d = %3.2f\n', S.crossval_accuracy, S.classification_d_singleinterval);

        if S.mult_obs_within_person

            fprintf('                            Forced-choice within-person: %3.0f%% d = %3.2f\n', S.crossval_accuracy_within, S.classification_d_within);

        end

    end

    % Plot scores and ROC curves
    % -------------------------------------------------------------------------
    if doplot

        plot_scores_and_ROC(S, 'Cross-validated SVM scores with hyperparameter optimization');

    end

end % do optimize

%% Repeat cross-validation of optimized model, if specified
% - if hyperparams optimized, then use "consensus params" for repeats
% - these are stored in S.modeloptions now, so will be used
% - If optimizing, re-do the first accuracy with the optimized
% hyperparameters
% -------------------------------------------------------------------------
if dorepeats > 1

    if doverbose
        disp(' ')
        fprintf('Repeating cross-val %d times. ', dorepeats)
    end

    % if dooptimize, startat = 1; else, startat = 2; end
    startat = 2;

    for i = startat:dorepeats     % do remainder of repeats without verbose.

        if doverbose, fprintf('%d ', i), end

        Sr = S;

        % Get new cvpartition (or outer loop if optimizing)
%         [Sr.trIdx, Sr.teIdx] = xval_stratified_holdout_leave_whole_subject_out(Sr.Y, Sr.id, 'doverbose', false, 'doplot', false);
        [Sr.trIdx, Sr.teIdx] = xval_stratified_holdout_leave_whole_subject_out(Sr.Y, Sr.id, 'doverbose', false, 'doplot', false, 'nfolds', nfolds); % MS 6/16/2022 add nfold option, otherwise default nfold=10 will cause this to crash if there are fewer than 10 possible folds.

        if dooptimize
            % If optimizing, repeat nested xval procedure
            % This can take a long time!

            Sr = crossval_nested(Sr, X, doverbose);

        else

            % Cross-validate
            Sr = crossval(Sr, X, false);

        end

        % Save accuracy
        % -------------------------------------------------------------------------
        % Given Y, yfit, cross-val scores (dist_from_hyperplane_xval), id
        % Update crossval_accuracy,classification_d_singleinterval
        % Y_within_id, scores_within_id, scorediff,
        % crossval_accuracy_within, classification_d_within

        Sr = update_accuracy_stats(Sr);

        S.crossval_accuracy(i) = Sr.crossval_accuracy; % Sr.accfun(Sr.Y, Sr.yfit);
        S.classification_d_singleinterval(i) = Sr.classification_d_singleinterval;

        if S.mult_obs_within_person

            S.crossval_accuracy_within(i) = Sr.crossval_accuracy_within; % Sr.accfun(Sr.Y, Sr.yfit);
            S.classification_d_within(i) = Sr.classification_d_within;

        end

    end

    if doverbose

        fprintf('Done!\n')

        fprintf('CV single-interval accuracy across %d reps, mean = %3.0f%%, std = %3.0f%%, min = %3.0f%%, max = %3.0f%%\n', dorepeats, ...
            mean(S.crossval_accuracy), std(S.crossval_accuracy), min(S.crossval_accuracy), max(S.crossval_accuracy));

        disp('Performance metrics saved in S.crossval_accuracy, classification_d_singleinterval')

        if S.mult_obs_within_person

            fprintf('\nCV forced_choice accuracy across %d reps, mean = %3.0f%%, std = %3.0f%%, min = %3.0f%%, max = %3.0f%%\n', dorepeats, ...
                mean(S.crossval_accuracy_within), std(S.crossval_accuracy_within), min(S.crossval_accuracy_within), max(S.crossval_accuracy_within));

             disp('Performance metrics saved in S.crossval_accuracy_within, classification_d_within')

        end

        disp(' ')

    end

end

%% Fit model on full dataset with selected options (final model if optimizing hyperparams) to get betas

S.SVMModel = fitcsvm(X, S.Y, S.modeloptions{:});
S.w = S.SVMModel.Beta;                                       % Model weights/betas

% Bootstrap betas
if dobootstrap

    if doverbose
        disp(' ')
        printhdr(sprintf('Bootstrapping model param estimates (b): %d samples', nboot));
    end

    rng('shuffle');

    S.boot_w = NaN .* zeros(length(S.w), nboot);

    for i = 1:nboot

        [Xb, Yb] = get_bootstrap_sample_grouped_by_id(S, X);

        SVMModel = fitcsvm(Xb, Yb, S.modeloptions{:});

        S.boot_w(:, i) = SVMModel.Beta;                                       % Model weights/betas

    end

    % Inference
    % (from fmri_data.predict)

    S.boot_w_ste = squeeze(nanstd(S.boot_w, 0, 2)); %1/20/16 add squeeze for multiclass case
    S.boot_w_mean = squeeze(nanmean(S.boot_w, 2)); %1/20/16 add squeeze for  multiclass case
    S.boot_w_ste(S.boot_w_ste == 0) = Inf;  % in case unstable regression returns all zeros

    S.wZ = S.boot_w_mean ./ S.boot_w_ste;  % tor changed from wmean; otherwise bootstrap variance in mean inc in error; Luke renamed to avoid confusion
    S.wP = 2 * (1 - normcdf(abs(S.wZ)));
    S.wP_fdr_thr = FDR(S.wP, .05);
    if isempty(S.wP_fdr_thr), S.wP_fdr_thr = -Inf; end

    S.boot_w_fdrsig = S.wP <= S.wP_fdr_thr;  % equals because can get exact vals in some cases...
    S.w_thresh_fdr = S.w;
    S.w_thresh_fdr(~S.boot_w_fdrsig) = 0;

    if doverbose

        disp('Summary of significant individual features (two-tailed)');

        Threshold = [.05 .01 .001 S.wP_fdr_thr]';

        for i = 1:length(Threshold)

            Num_Sig_Features(i, 1) = sum(S.wP <= Threshold(i));

        end

        t = table(Threshold, Num_Sig_Features);
        disp(t)

    end
end

%%
% permutations would go here...future project.
% consider exchangeability issue given within-person design...winkler and
% nichols...
% indx = randperm(size(X, 1));
% X = X(indx, :); Y = Y(indx); id = id(indx);
% S = xval_SVM(X, Y, id, 'nooptimize', 'norepeats');

end % main function

%%

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

% Subfunctions

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

function S = crossval(S, X, doverbose)

if doverbose, fprintf('..X-val, %d folds...', S.nfolds), end

for i = 1:S.nfolds

    if doverbose, fprintf('%d ', i), end

    % Fit to training data for this fold
    SVM_fold = fitcsvm(X(S.trIdx{i}, :), S.Y(S.trIdx{i}), S.modeloptions{:});

    % Apply holdout test set for this fold
    [label, score] = predict(SVM_fold, X(S.teIdx{i}, :));   % Get raw scores
    score = score(:, 2);                                    % Decision boundary is symmetrical

    S.dist_from_hyperplane_xval(S.teIdx{i}, 1) = score;     % Unscaled SVM scores

    S.yfit(S.teIdx{i}, 1) = label;                          % Predicted class, cross-val

    SVM_fold = fitPosterior(SVM_fold);                      % Works for fitcsvm, not fitclinear

    % Note:
    % either scores or pscores appear to be inconsistent in terms of which
    % column is class -1 and which is class 1. Check for instability in
    % scores, and select pscores to be consistent with this on each fold.
    if mean(label(score > 0)) < mean(label(score < 0))
        disp('WARNING!!! Scores for label 1 are < scores for label -1, scores are reversed!!! This should not happen. Check code/implementation.');
    end

    [~, pscore] = predict(SVM_fold, X(S.teIdx{i}, :));

    % which pscores correlate more positively with score? use this one.
     r = corr([score pscore]);
    [~, wh] = max(r(1, 2:3));

    pscore = pscore(:, wh);                                  % Platt scaling scores (class probability)

    S.class_probability_xval(S.teIdx{i}, 1) = pscore;

end

if doverbose, fprintf('Done!\n'), end

end % crossval

% -------------------------------------------------------------------------
% Nested cross-validation to optimize hyperparameters
% -------------------------------------------------------------------------

function S = crossval_nested(S, X, doverbose)

if doverbose, fprintf('..X-val, %d folds...', S.nfolds), end

for i = 1:S.nfolds

    if doverbose, fprintf('%d ', i), end

    % Optimize hyperparameters within this fold
    % -------------------------------------------

    % Option 1: Optimize all, including kernel (linear/nonlinear)
    % Option 2: Optimize specific hyperparameters, Slack param C
    % (BoxConstraint) and Standardize inputs, linear model only

    hparams_to_optimize = 'all'; % {'BoxConstraint' 'Standardize'}

    Mdl = fitcsvm(X(S.trIdx{i}, :), S.Y(S.trIdx{i}),...
        'OptimizeHyperparameters', hparams_to_optimize, 'HyperparameterOptimizationOptions',...
        struct('AcquisitionFunctionName','expected-improvement-plus', 'Verbose', 0, 'ShowPlots', 0));

    % Notes: HyperparameterOptimizationOptions include defaults of:
    % - 5-fold CV (not grouped by id) without repartitioning
    % - Bayesian optimization, using best estimates from smooth function

    % Add hyperparameter results to table
    % -------------------------------------
    % Use estimated rather than observed. This uses the smooth function estimated across observations

    opt_hyperparam_table = Mdl.HyperparameterOptimizationResults.XAtMinEstimatedObjective;

    if i == 1
        S.hyperparams_by_fold = opt_hyperparam_table; % best parameters - table

    else
        S.hyperparams_by_fold(i, :) = opt_hyperparam_table;
    end

    fold_modeloptions = convert_hyperparameter_choices_to_cell(Mdl, []); % Optimizing "all"
    %fold_modeloptions = convert_hyperparameter_choices_to_cell(Mdl, S.modeloptions); % Optimizing "linear"

    S.fold_modeloptions{i} = fold_modeloptions;

    % Fit to training data for this fold
    SVM_fold = fitcsvm(X(S.trIdx{i}, :), S.Y(S.trIdx{i}), fold_modeloptions{:});

    % Updates S.yfit, S. S.dist_from_hyperplane_xval, S.class_probability_xval

    % Apply holdout test set for this fold
    [label, score] = predict(SVM_fold, X(S.teIdx{i}, :));   % Get raw scores
    score = score(:, 2);                                    % Decision boundary is symmetrical

    S.dist_from_hyperplane_xval(S.teIdx{i}, 1) = score;     % Unscaled SVM scores

    S.yfit(S.teIdx{i}, 1) = label;                          % Predicted class, cross-val

    SVM_fold = fitPosterior(SVM_fold);                      % Works for fitcsvm, not fitclinear

    [~, pscore] = predict(SVM_fold, X(S.teIdx{i}, :));
    pscore = pscore(:, 2);                                  % Platt scaling scores (class probability)

    S.class_probability_xval(S.teIdx{i}, 1) = pscore;

end

if doverbose, fprintf('Done!\n'), end

end % crossval

% -------------------------------------------------------------------------
% Get scores
% -------------------------------------------------------------------------

function [scores_within_id, scorediff, d] = get_scores_within_id(S, varname)

myvar = S.(varname);
u = unique(S.id);

for i = 1:length(u)

    wh_id = S.id == u(i) & S.Y == -1;
    scores_within_id(i, 1) = nanmean(myvar(wh_id));

    wh_id = S.id == u(i) & S.Y == 1;
    scores_within_id(i, 2) = nanmean(myvar(wh_id));

end

scorediff = diff(scores_within_id')';
d = nanmean(scorediff) ./ nanstd(scorediff);

end %get scores

% -------------------------------------------------------------------------
% bootstrapping
% -------------------------------------------------------------------------

function [Xb, Yb] = get_bootstrap_sample_grouped_by_id(S, X)
% Given X, S.Y, S.id, return a bootstrap sample, keeping all obs from an id together

u = unique(S.id);
n = length(u);

wh = ceil(rand(n, 1) * n);  % Bootstrap sample of p ids, with replacement

[Xboot, Yboot] = deal(cell(n, 1));

for j = 1:n                 % add obs, allowing repeats

    wh_obs = ismember(S.id, u(wh(j)));

    Xboot{j} = X(wh_obs, :);
    Yboot{j} = S.Y(wh_obs);

end

Xb = cat(1, Xboot{:});
Yb = cat(1, Yboot{:});

end % function

% -------------------------------------------------------------------------
% optimized hyperparameter aggregation
% -------------------------------------------------------------------------

function  best_modeloptions = convert_hyperparameter_choices_to_cell(Mdl, other_modeloptions)
% Notes:
% - Extract hyperparameter choices in a format that we can pass back in
% to train or test a new model.
% - We could use the best observed point, but it's preferred to use the
% best estimated values from a smooth model fit to observed samples.

% opt_hyperparam_table = Mdl.HyperparameterOptimizationResults.XAtMinObjective;
opt_hyperparam_table = Mdl.HyperparameterOptimizationResults.XAtMinEstimatedObjective;

best_modeloptions = other_modeloptions;

% Add optimal choices to the option set
for j = 1:size(opt_hyperparam_table, 2)
    vname = opt_hyperparam_table.Properties.VariableNames{j};

    paramval = opt_hyperparam_table.(vname)(1);

    % fix: categorical Standardize to text.
    % Leave KernelFunction as string
    if iscategorical(paramval) && strcmp(vname, 'Standardize')
        paramval = char(string(paramval));
        paramval = strcmp(paramval, 'true');
    elseif iscategorical(paramval)
        % For KernelFunction
        paramval = char(paramval);
    end

    % Get rid of NaN options passed back out from optimize 'all'
    if ~isnan(paramval)
        best_modeloptions{end + 1} = vname;
        best_modeloptions{end + 1} = paramval; %#ok<*AGROW>
    end

end

end

% function best_modeloptions = get_modal_params_across_folds(S)
%
% % This function is deprecated because modal parameters across folds is not
% % the best way to get the final hyperparameters. Instead, the model should
% % be re-fit using single cross-validation with all hparam options pick the
% % best relative model. Nested cross-val will estimate the accuracy of the
% % whole procedure, including the hparam search, but not return the final best
% % parameters.
%
% hbyfold = cat(1, S.fold_modeloptions{:});
% best_modeloptions = {};
%
% for i = 1:size(hbyfold, 2)
%
%     mydat = cat(1, hbyfold{:, i});
%
%     if ischar(mydat) || iscategorical(mydat) || islogical(mydat)
%
%         best_modeloptions{1, i} = mode(mydat);
%
%     elseif isnumeric(mydat)
%
%         best_modeloptions{1, i} = trimmean(mydat, 80);
%
%     else
%         error('Unknown model option class! Extend this code.')
%
%     end
%
% end
%
% end

% -------------------------------------------------------------------------
% accuracy
% -------------------------------------------------------------------------

function S = update_accuracy_stats(S)

S.crossval_accuracy = S.accfun(S.Y, S.yfit);

s1 = S.dist_from_hyperplane_xval(S.Y == 1);
s2 = S.dist_from_hyperplane_xval(S.Y == -1);

% d = difference / std pooled within, weighted by sample size
S.classification_d_singleinterval = (mean(s1) - mean(s2)) ./ sqrt( ( var(s1) .* length(s1) + var(s2) .* length(s2) ) ./ (length(s1) + length(s2)) );

% Do we have multiple obs within-person? If so return within-person stats
S.mult_obs_within_person = length(S.id) > length(unique(S.id));

if S.mult_obs_within_person

    varname = 'dist_from_hyperplane_xval';
    [scores_within_id, scorediff, d] = get_scores_within_id(S, varname);

    S.Y_within_id = get_scores_within_id(S, 'Y');

    S.scores_within_id = scores_within_id;
    S.scorediff = scorediff;
    S.crossval_accuracy_within = 100 * sum(scorediff > 0) ./ sum(~isnan(scorediff));
    S.classification_d_within = d;

end

end

% -------------------------------------------------------------------------
% plots
% -------------------------------------------------------------------------

function prelim_data_plot(X, Y, id)

create_figure('Data view', 1, 3);
imagesc(X); colorbar;
xlabel('Input features'); ylabel('Observation'); title('Predictor matrix (X)');
axis tight; set(gca, 'YDir', 'reverse');

subplot(1, 3, 2);
imagesc([Y scale(id)]);
axis tight; set(gca, 'YDir', 'reverse');
set(gca, 'XTick', [1 2], 'XTickLabel', {'Y' 'id'});
title('Outcome (Y) and id');

subplot(1, 3, 3);
r = corr(X');
[~, wh] = sort(Y);
imagesc(r(wh, wh));
n1 = sum(Y == -1);
n2 = sum(Y == 1);

h1 = drawbox(0, n1, 0, n1, 'k');
set(h1, 'FaceColor', 'none', 'EdgeColor', 'r', 'LineWidth', 2);

h1 = drawbox(n1, n2, n1, n2, 'k');
set(h1, 'FaceColor', 'none', 'EdgeColor', 'r', 'LineWidth', 2);

axis tight; set(gca, 'YDir', 'reverse');
% set(gca, 'XTick', [1 2], 'XTickLabel', {'Y' 'id'});
title('Inter-obs correlations sorted by Y [-1, 1]');

colorbar

cm = colormap_tor([0 0 1], [1 0 0], [1 1 1]);
colormap(cm)
set(gca, 'CLim', [-1 1])

drawnow, snapnow
end

function plot_scores_and_ROC(S, varargin)
% plot_scores_and_ROC(S, varargin) -> varargin = new title

create_figure('cross-val accuracy', 2, 2);
subplot(2, 2, 1); delete(gca); subplot(2, 2, 2); delete(gca);
axes('Position', [.13 .58 .77 .35]);
set(gca, 'FontSize', 16)
hold on

if S.mult_obs_within_person
    % Within-person scores

    n = size(S.scores_within_id, 1);
    x = [1:n; 1:n];
    plot(x(:, S.scorediff > 0), S.scores_within_id(S.scorediff > 0, :)', 'Color', [.3 .3 .3], 'LineWidth', 1);
    plot(x(:, S.scorediff < 0), S.scores_within_id(S.scorediff < 0, :)', 'r', 'LineWidth', 1);

    plot(x(1, :), S.scores_within_id(:, 1)', '^', 'Color', [.3 .5 1] / 2, 'MarkerFaceColor', [.3 .5 1],  'LineWidth', 1);
    plot(x(2, :), S.scores_within_id(:, 2)', 'v', 'Color', [1 .5 0] / 2, 'MarkerFaceColor', [1 .5 0],  'LineWidth', 1);

    xlabel(sprintf('ID, %3.0f%% single-interval acc, %3.0f%% forced-choice', S.crossval_accuracy, S.crossval_accuracy_within));

    disp('Black lines: Correct, Red lines: Errors. Red triangles: Scores for Class 1, Blue triangles: Scores for Class -1');

else

    % Scores, all between-person

    n = size(S.dist_from_hyperplane_xval, 1);
    x = 1:n;
    iscorrect = sign(S.dist_from_hyperplane_xval) == sign(S.Y);

    %plot(x(:, iscorrect), S.dist_from_hyperplane_xval(iscorrect, :)', 'o', 'MarkerSize', 10, 'Color', [.3 .3 .3], 'LineWidth', 1);
    plot(x(:, ~iscorrect), S.dist_from_hyperplane_xval(~iscorrect, :)', 'o', 'MarkerSize', 10, 'Color', [1 0 0], 'LineWidth', 1);

    h = plot_horizontal_line(0); set(h, 'LineStyle', '--');

    color1 = [.2 .8 .2];
    color2 = [.3 .3 1];

    plot(x(:, S.Y == 1), S.dist_from_hyperplane_xval(S.Y == 1, :)', 'v', 'MarkerSize', 6, 'Color', color1, 'MarkerFaceColor', color1, 'LineWidth', 1);
    plot(x(:, S.Y == -1), S.dist_from_hyperplane_xval(S.Y == -1, :)', '^', 'MarkerSize', 6, 'Color', color2, 'MarkerFaceColor', color2, 'LineWidth', 1);

    xlabel(sprintf('ID, %3.0f%% single-interval acc', S.crossval_accuracy));

    disp('Red cicles: Errors. Green triangles: Scores for Class 1, Blue triangles: Scores for Class -1');

    axis tight

end

ylabel('SVM Score');

title('Cross-validated SVM scores (no hyperparameter optimization)');
if ~isempty(varargin), title(varargin{1}); end

% ROC plot
subplot(2, 2, 3);
S.ROC_single_interval = roc_plot(S.dist_from_hyperplane_xval, logical(S.Y > 0), 'color', [.4 .4 .7], 'threshold', 0);
title('Single-interval ROC')
set(gca, 'FontSize', 16)

subplot(2, 2, 4);

if S.mult_obs_within_person

    % Paired forced-choice. Get complete cases - Remove NaNs id-wise
    outcomes = S.Y_within_id;
    scores = S.scores_within_id;
    [~, outcomes, scores] = nanremove(outcomes, scores);

    S.ROC_forced_choice = roc_plot(scores(:), logical(outcomes(:) > 0), 'color', [.4 .4 .7], 'twochoice');
    title('Forced-choice ROC')
    set(gca, 'FontSize', 16)

else
    % Scores, all between-person

    plot(S.Y(~iscorrect), S.dist_from_hyperplane_xval(~iscorrect, :)', 'o', 'MarkerSize', 10, 'Color', [1 0 0], 'LineWidth', 1);

    h = plot_horizontal_line(0); set(h, 'LineStyle', '--');

    color1 = [.2 .8 .2];
    color2 = [.3 .3 1];

    plot(ones(1, sum(S.Y == 1)), S.dist_from_hyperplane_xval(S.Y == 1, :)', 'v', 'MarkerSize', 6, 'Color', color1, 'MarkerFaceColor', color1, 'LineWidth', 1);
    plot(-ones(1, sum(S.Y == -1)), S.dist_from_hyperplane_xval(S.Y == -1, :)', '^', 'MarkerSize', 6, 'Color', color2, 'MarkerFaceColor', color2, 'LineWidth', 1);

    set(gca, 'XTick', [-1 1], 'XLim', [-1.5 1.5]);
    xlabel('True class');
    ylabel('Predicted class')

end

drawnow, snapnow

end

function xval_plot_scores_vs_class_probability(S)

% crossval returns cross-validated yfit (class predictions), scores (dist_from_hyperplane_xval),
% class_probability_xval. Scores and class probabilties will diverge (may
% not be perfectly correlated) because the sigmoid scaling (Platt scaling) varies across folds.
% They will diverge more if there is no true signal.

create_figure('xval scores vs. class probability estimates', 1, 2);

plot(S.dist_from_hyperplane_xval, S.class_probability_xval, 'o');

for j = 1:length(S.teIdx)

    plot(S.dist_from_hyperplane_xval(S.teIdx{j}), S.class_probability_xval(S.teIdx{j}), 'o', 'MarkerFaceColor', rand(1, 3));

end

xlabel('Cross-validated SVM scores (colors are folds)');
ylabel('Cross-validated class prob estimates');

hh = plot_vertical_line(0); set(hh, 'LineStyle', '--');
hh = plot_horizontal_line(.5); set(hh, 'LineStyle', '--');

subplot(1, 2, 2)

lineh = plot(S.dist_from_hyperplane_xval(S.Y == 1), S.class_probability_xval(S.Y == 1), 'v', 'MarkerFaceColor', [.2 .8 .2]);
lineh = [lineh plot(S.dist_from_hyperplane_xval(S.Y == -1), S.class_probability_xval(S.Y == -1), '^', 'MarkerFaceColor', [.2 .2 1])];

xlabel('Cross-validated SVM scores (colors are folds)');
ylabel('Cross-validated class prob estimates');

hh = plot_vertical_line(0); set(hh, 'LineStyle', '--');
hh = plot_horizontal_line(.5); set(hh, 'LineStyle', '--');

legend(lineh, {'True Class = 1' 'True Class = -1'});

drawnow, snapnow

end

% -------------------------------------------------------------------------
% Inputs
% -------------------------------------------------------------------------

function ARGS = parse_inputs(varargin)

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

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

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

valfcn_logical = @(x) validateattributes(x, {}, {'nonempty', 'scalar', '>=', 0, '<=', 1}); % could enter numeric 0,1 or logical

valfcn_effectscode = @(x) validateattributes(x, {'numeric'}, {'nonempty', '<=', 1, '>=', -1});

% Required inputs
% ----------------------------------------------------------------------
p.addRequired('X', valfcn_number);
p.addRequired('Y', valfcn_effectscode);
p.addRequired('id', valfcn_number);

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

p.addParameter('doplot', true, valfcn_logical);
p.addParameter('doverbose', true, valfcn_logical);
p.addParameter('dooptimize', true, valfcn_logical);
p.addParameter('dorepeats', 10, valfcn_number);
p.addParameter('dobootstrap', true, valfcn_logical);
p.addParameter('nboot', 1000, valfcn_number);
p.addParameter('nfolds',10, valfcn_number); % Added by Michael Sun 08/2/2022

p.addParameter('modeloptions', {'KernelFunction', 'linear'}, valfcn_cell);

% 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{:});

ARGS = p.Results;

end % parse_inputs);

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
internet    = load('data/internet.dat')';
y           = internet(2, 2:end);
fprintf(1, '\n');

P       = 5;
Q       = 5;

%% Model selection for complete series %%
fprintf(1, 'AIC for ARMA(p, q) models on complete data:\n');
fprintf(1, '    '); for q = 0:Q, fprintf(1, '%-12d', q); end; fprintf(1, '\n');
logL    = zeros(P+1, Q+1);
AIC     = zeros(P+1, Q+1);
arma    = cell(P+1, Q+1);
for p = 0 : P
    fprintf(1, '%-4d', p);
    for q = 0 : Q
        model                                   = ssm_arma(p, q);
        [arma{p+1, q+1} logL(p+1, q+1) output]  = estimate(y, model, [repmat(0.3, 1, model.w-1) 10], [], 'fmin', 'bfgs', 'disp', 'off');
        AIC(p+1, q+1)                           = output.AIC;
        fprintf(1, '%-12g', AIC(p+1, q+1));
    end
    fprintf(1, '\n');
end
[m i] = min(AIC(:)); temp = AIC(i); AIC(i) = Inf;
fprintf(1, 'ARMA(%d, %d) found to be the best model, ', mod(i-1, P+1), floor((i-1)/(P+1)));
[m j] = min(AIC(:)); AIC(i) = temp;
fprintf(1, 'ARMA(%d, %d) is the second best.\n\n', mod(j-1, P+1), floor((j-1)/(P+1)));

%% Model selection for data with missing values %%
ymis        = y; ymis([6 16 26 36 46 56 66 72 73 74 75 76 86 96]) = NaN;
fprintf(1, 'AIC for ARMA(p, q) models on data w/ missing observations:\n');
fprintf(1, '    '); for q = 0:Q, fprintf(1, '%-12d', q); end; fprintf(1, '\n');
logLmis     = zeros(P+1, Q+1);
AICmis      = zeros(P+1, Q+1);
armamis     = cell(P+1, Q+1);
for p = 0 : P
    fprintf(1, '%-4d', p);
    for q = 0 : Q
        model                                           = ssm_arma(p, q);
        [armamis{p+1, q+1} logLmis(p+1, q+1) output]    = estimate(ymis, model, [repmat(-0.1, 1, model.w-1) 1], [], 'fmin', 'bfgs', 'disp', 'off');
        AICmis(p+1, q+1)                                = output.AIC;
        fprintf(1, '%-12g', AICmis(p+1, q+1));
    end
    fprintf(1, '\n');
end
[m i] = min(AICmis(:)); temp = AICmis(i); AICmis(i) = Inf;
fprintf(1, 'ARMA(%d, %d) found to be the best model, ', mod(i-1, P+1), floor((i-1)/(P+1)));
[m j] = min(AICmis(:)); AICmis(i) = temp;
fprintf(1, 'ARMA(%d, %d) is the second best.\n\n', mod(j-1, P+1), floor((j-1)/(P+1)));

%% Forecast with ARMA(1, 1) on the complete data %%
[armafore logL] = estimate(y, ssm_arma(1, 1), 0.1);
yf              = signal(kalman([y repmat(NaN, 1, 20)], armafore), armafore);
figure('Name', 'Internet series forecast');
plot(yf, 'DisplayName', 'forecast'), hold all, scatter(1:length(y), y, 10, 'r', 's', 'filled', 'DisplayName', 'data'), hold off, ylim([-15 15]), legend('show');
if ispc, set(gcf, 'WindowStyle', 'docked'); end

%% Forecast with ARMA(1, 1) on the data w/ missing values %%
[armafore logLmis]  = estimate(ymis, ssm_arma(1, 1), 0.1);
ymisf               = signal(kalman([ymis repmat(NaN, 1, 20)], armafore), armafore);
figure('Name', 'Internet series in-sample one-step and out-of-sample forecasts');
plot(ymisf), hold all, scatter(1:length(ymis), ymis, 10, 'r', 's', 'filled'), hold off, ylim([-15 15]);
if ispc, set(gcf, 'WindowStyle', 'docked'); end

repo: REST


text:
function [err, Qd, s, termname, nterms, sindices, dfbothSS, modw, modwo, tnames, dfterm, dfe, Contr] = PreProc(n, NF, group, varnames, FL, Contr, cov, unbalanced)
%
%   [err, Qd, Rd, sindices, dfboth, modw, modwo, dfx] = PreProc(n,group,varnames)
%
%Purpose:
%
%  Returns QR decomposition results for the design matrix projected to null space
%
%Input Parameters:
%   n: number of datasets = total # of combinations including repeats
%   group: cell array of factor levels
%   varnames:
%
%Output Parameters:
%   err : 0 No Problem
%       : 1  Problems
%
%
%
%Key Terms:
%
%More Info :
%
%
%
%
%     Author : Gang Chen
%     Date : Tue Mar 23 13:57:52 EST 2004
%     SSCC/NIMH/ NIH, Bethesda MD 20892

%Define the function name for easy referencing
FuncName = 'PreProc.m';

%Debug Flag
DBG = 1;

%initailize return variables
err = 1;

% Don't worry about NaN's at this point.
group = group(:);    % what's it for?
ng = length(group);  % number of factors
termlist = makemodel(ng, ng);  % Generate terms for all main effects plus various interactions

for j=1:ng
   gj = group{j};
   if (size(gj,1) == 1), gj = gj(:); end
   if (size(gj,1) ~= n)
      error('Factor %d must have %d elements.',j,n);
   end
   if (ischar(gj)), gj = cellstr(gj); end
   group{j} = gj;
end

gdum = cell(ng,1);
dfvar = zeros(ng,1);

if (unbalanced.yes == 0), % Balanced designs
   vconstr = cell(ng,1);
   %vmean = cell(ng,1);   %vmean is never used!!!
end

for j=1:ng   % for each factor
   gj = group{j};
   [gij,gnj] = grp2idx(gj);   % Create index vector from a grouping variable: gij is a vector
	                           % taking integer values from 1 up to the number of unique entries in gj
										% gnj is a cell array of names, so that gnj(gij) reproduces gj
   nlevels = size(gnj,1);     % levels for this factor
   dfvar(j) = nlevels - 1;    % D. F. for this factor

   if (unbalanced.yes == 0),	% balanced
      if (cov.do & j==cov.marker)
	      gdum{j} = gj;
         dfvar(j) = 1;           % D. F. = 1
         vconstr{j} = zeros(0,1);
%        vmean{j} = 1;
      else
         gdum{j} = idummy(gij, 3);
         vconstr{j} = idummy(1:nlevels)';
%        vmean{j} = ones(1,nlevels) / nlevels;  % array (1Xnlevels) of one ones, but vmean is never used in the code!!!!!!!!!1
	   end
	else % Unbalanced designs
	   if (cov.do & j==cov.marker)
	      gdum{j} = gj;
         dfvar(j) = 1;           % D. F. = 1
         vconstr{j} = zeros(0,1);
%        vmean{j} = 1;
      else
		   gdum{j} = idummy(gij, 3);
		end
	end % if (unbalanced.yes == 0): Only for balanced designs

end

% Create dummy variable arrays for each term in the model.
nterms = size(termlist,1);             % Number of rows (1st dimension) in termlist
[sterms,sindex] = sortrows(termlist);  % Sort terms in ascending order.
ncols = 1;
nconstr = 0;

termdum = cell(nterms, 1);        % cell array of dummy variables which are design matrix cols
termconstr = cell(nterms,1);      % constraints to make each term well defined
levelcodes = cell(nterms, 1);     % codes for levels of each M row
tnames = cell(nterms, 1);         % name for entire term, e.g. A*B
dfterm0 = zeros(nterms, 1);       % nominal d.f. for each term
termvars = cell(nterms, 1);       % list of vars in each term
termlength = zeros(size(sindex)); % length of each term (number of columns)

%randomterm = find(termlist*randomvar > 0);  % indices of random terms

% For each term,
for j=1:nterms
   % Loop over elements of the term
   df0 = 1;
   tm = sterms(j,:);
   tdum = [];         % empty term so far
   tconstr = 1;       % empty constraints so far
   tn = '';           % blank name so far
   vars = find(tm);   % Find indices of nonzero elements
   for varidx = 1:length(vars)
      % Process each variable participating in this term
      varnum = vars(varidx);          % variable name
      tm(varnum) = 0;                 % term without this variable
      df0 = df0 * dfvar(varnum);      % d.f. so far

      % Combine its dummy variable with the part computed so far
      G = gdum{varnum};           % dummy vars for this grouping var
      nlevterm = size(tdum,2);    % levels for term so far
      nlevgrp  = size(G,2);       % levels for this grouping var
      tdum = termcross(G,tdum);   % combine G into term dummy vars

      % Construct the term name and constraints matrix
%      if (ismember(varnum, NF) & cov.do),    % for the covariate, which is the last factor
		if (ismember(varnum, cov.marker) & cov.do),    % for the covariate, which is the last factor
         vconstr = ones(0,1);
      else
         vconstr = ones(1, nlevgrp);
      end
      if (isempty(tn))
         tn = varnames{varnum};
         tconstr = vconstr;
      else
         tn = [varnames{varnum} '*' tn];
         tconstr = [kron(vconstr,eye(size(tconstr,2)));
                    kron(eye(length(vconstr)),tconstr)];   % Kronecker
      end

      % If the rest of this term is computed, take advantage of that
      prevterms = sterms(1:j-1,:);
      oldtm = find((prevterms * tm') == sum(tm));      % same vars in old term
      oldtm = oldtm((prevterms(oldtm,:) * ~tm') == 0); % and no others
      if (length(oldtm) > 0)
         k = sindex(oldtm(1));
         tdum = termcross(termdum{k}, tdum);
         oconstr = termconstr{k};
         tconstr = [kron(tconstr,              eye(size(oconstr,2)));
                    kron(eye(size(tconstr,2)), oconstr)];
         tn = [tn '*' tnames{k}];
         df0 = df0 * dfterm0(k);
         break;
      end
   end

   % Store this term's dummy variables and name
   k = size(tdum, 2);
   termlength(sindex(j),1) = k;
   ncols = ncols + k;
   sj = sindex(j);
   termdum{sj} = tdum;
   termconstr{sj} = tconstr;
   termvars{sj} = vars;
   levelcodes{sj} = fullfact(dfvar(vars)+1);
   if (isempty(tn)), tn = 'Constant'; end
   tnames{sj,1} = tn;
   dfterm0(sj) = df0;
   nconstr = nconstr + size(tconstr,1);
end
tnames{length(tnames)+1,1} = 'Error';

% Create the full design matrix
dmat = ones(n, ncols);        % to hold design matrix of nXncols
cmat = zeros(nconstr,ncols);  % to hold constraints matrix
cbase = 0;                    % base from which to fill in cmat
termname = zeros(ncols,1);
termstart = cumsum([2; termlength(1:end-1)]);
termend = termstart + termlength - 1;
for j=1:nterms
   clist = termstart(j):termend(j);
   dmat(:, clist) = termdum{j};
   C = termconstr{j};
   nC = size(C,1);
   cmat(cbase+1:cbase+nC,clist) = C;
   termname(clist) = j;
   cbase = cbase + nC;
end

[err, Qd, dfx, dmat2] = QRDecom(dmat, cmat);
%dfe = n - dfx;

% Determine which models to compare for testing each term
ssw  = -ones(nterms, 1);      % sum of squares with this term
sswo = ssw;                   % sum of squares without this term
dfw  = ssw;                   % residual d.f. with this term
dfwo = ssw;                   % residual d.f. without this term

if (unbalanced.yes == 1)
   modw = tril(ones(nterms)); % get model with this term
   k = nterms;                % locations of model with all terms
else

% Only apply type III for sums of squares
modw = ones(nterms);
%k = 1:nterms;
%TnotC = termsnotcontained(termlist);

end % if (unbalanced.yes == 1)

modw = logical(modw);                  % get model with this term
modwo = logical(modw - eye(nterms));   % get model without this term

dfw(1:nterms) = dfx;

% Fit each model separately
dfboth = [dfw; dfwo];
dfbothSS = dfboth;  % For usage in ss.m

% Consider interactions before their components for type 3 ss, so
% examine the terms in decreasing order of the number factors in the term

if (unbalanced.yes == 1)
   sindices = [(1:size(termlist,1)), (1:size(termlist,1))]';
else
termsize = sum(termlist,2)';      %sum of all the elements along 2nd dimension (row)
[stermsize,sindices] = sort(termsize); % sort in ascending order, output is stored in stermsize; sindices is the index for the new array
sindices = [sindices(:); sindices(:)];

end % if (unbalanced.yes == 1)

% Here QR decomposition is done, which is voxel independent

s(length(sindices)).Qdt = [];

for j=length(sindices):-1:1
   % Find the next model index to fit
   k = sindices(j);

   % Look in unsorted arrays to see if we have already fit this model
   if j>nterms
      k0 = k+nterms;
   else
      k0 = k;
   end
   if dfboth(k0)~=-1
      continue
   end

   % Find the model with this index
   if (j > nterms)
      thismod = modwo(k, :);
   else
      thismod = modw(k, :);
   end

   % Get the design matrix for this model
   keepterms = find(thismod);
   clist = ismember(termname, [0 keepterms]);
   X = dmat2(:,clist);
   C = cmat(:,clist);

   % Fit this term
	[err, s(j).Qdt, dfx0] = QRDecom(X, C);

   % Use these results for each term that requires them

   mod0 = repmat(thismod, nterms, 1);
   k = find(all(modw == mod0, 2));
   dfw(k) = dfx0;
   dfboth(k) = 0;

   k = find(all(modwo == mod0, 2));
   dfwo(k) = dfx0;
   dfboth(nterms+k) = 0;
end
clear mod0

dfterm = dfw - dfwo;
dfe = n-dfx;   %residual degrees of freedom

if (Contr.do == 1),

   % In design matrix dmat, the first column is all one's, for the total mean. Then there are totally
   % FL(1).N_level + FL(2).N_level + FL(3).N_level + FL(4).N_level columns for the main effects.
   % Next 2nd order interactions, 3rd order interaction, and 4th order interactions.

   % Store only those mean columns in design matrix.

   %num_col0 = 1;  % the 1st column is for grand mean (0 order)

   num_col(1) = 0;
   for (i = 1:1:NF),
      num_col(1) = num_col(1) + FL(i).N_level;   %
   end

   %  Ignore 1st col since it is total mean
   %dmat_mean = dmat(:, (num_col0 + 1):(num_col0 + num_col(1)));

   % Get the number in the sum for each mean, which happens to be in the diagonal of X'X.
   % This can also be otained through the user input variables, but it is generic with
   % the matrix operation, especially for unbalanced design.
   %sum_num = diag(dmat_mean' * dmat_mean);

   if (NF > 1),
      num_col(2) = 0;
      for (i = 1:1:(NF-1)),
      for (j = (i+1):1:NF),
         num_col(2) = num_col(2) + FL(i).N_level*FL(j).N_level;    %Columns for 2nd order interactions
      end
      end
   end

   if (NF > 2),
      num_col(3) = 0;
      for (i = 1:1:(NF-2)),
      for (j = (i+1):1:(NF-1)),
   	for (k = (j+1):1:NF),
         num_col(3) = num_col(3) + FL(i).N_level*FL(j).N_level*FL(k).N_level;    %Columns for 3rd order interactions
      end
      end
   	end
   end

   %if (NF == 4),
   %   num_col(4) = 1;
   %	for (i = 1:1:NF),
   %	   num_col(4) = num_col(4)*FL(i).N_level;  %Columns for 4th order interactions
   %	end
   %end

   % for every design
   if (Contr.ord1.tot > 0),		% 1st order contrasts
      [err, Contr] = ContrVec(1, n, NF, group, dmat, Contr, FL, num_col);
   end   % if (Contr1.tot > 0)

   if (NF > 1 & Contr.ord2.tot > 0),  % 2nd order contrasts
      [err, Contr] = ContrVec(2, n, NF, group, dmat, Contr, FL, num_col);
   end   % if (Contr2.tot > 0)

   if (NF > 2 & Contr.ord3.tot > 0),  % 3rd order contrasts
      [err, Contr] = ContrVec(3, n, NF, group, dmat, Contr, FL, num_col);
   end   % if (Contr3.tot > 0)

	if (NF > 3 & Contr.ord4.tot > 0),  % 4th order contrasts
      [err, Contr] = ContrVec(4, n, NF, group, dmat, Contr, FL, num_col);
   end   % if (Contr3.tot > 0)

end % if (Contr.do == 1)

err = 0;
return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [X,data] = kamada_kawai_spring_layout(A,varargin)
% KAMADA_KAWAI_SPRING_LAYOUT Use springs to compute a graph layout
%
% X = kamada_kawai_spring_layout(A) computes a layout for an undirected
% graph, connected graph by simulating a system of springs.
% [X,data] = kamada_kawai_spring_layout(A) returns extra data, such as the
% spring matrix and matrix of distances between vertices, i.e.
%   data.spring_strength % returns the spring stremgth
%   data.distances % returns the distance between all points
%
% See
% http://www.boost.org/doc/libs/1_36_0/libs/graph/doc/kamada_kawai_spring_layout.html
% for information about the layout function and the parameters
%
% ... = kamada_kawai_spring_layout(A,...) takes a set of
% key-value pairs or an options structure.  See set_matlab_bgl_options
% for the standard options.
%   options.tol: stopping tolerance of layout change [double | {1e-6}]
%   options.maxiter: the maximum number of layout iterations [500]
%   options.spring_constant: energy of the system [double | {1}]
%   options.progressive: if you want to start from an existing layout,
%     provide the coordinates of the layout [{0} | position matrix X]
%   options.edge_length: default length of an edge [double | {1}]
%   options.edge_weight: a double array over the edges with an edge
%       weight for each edge, see EDGE_INDEX and EXAMPLES/REWEIGHTED_GRAPHS
%       for information on how to use this option correctly
%       [{'matrix'} | 'none' | length(nnz(A)) double vector]
%
% Example:
%   G = grid_graph(6,5);
%   X = kamada_kawai_spring_layout(G);
%   gplot(G,X);
%
% See also FRUCHTERMAN_REINGOLD_FORCE_DIRECTED_LAYOUT, GURSOY_ATUN_LAYOUT

% David F. Gleich
% Copyright, Stanford University, 2008

%% History
%  2008-09-25: Initial coding
%%

[trans check full2sparse] = get_matlab_bgl_options(varargin{:});
if full2sparse && ~issparse(A), A = sparse(A); end

n = num_vertices(A);
options = struct('tol',1e-6,'maxiter',500,'spring_constant',1,...
    'progressive',0,'edge_length',1,'edge_weight','matrix');
options = merge_options(options,varargin{:});

% edge_weights is an indicator that is 1 if we are using edge_weights
% passed on the command line or 0 if we are using the matrix, or -1 to use
% nothing
edge_weights = 0;
edge_weight_opt = 'matrix';

if strcmp(options.edge_weight, 'matrix'), % do nothing to use the matrix weights
elseif strcmp(options.edge_weight,'none'), edge_weights = -1; % don't use weights
else edge_weights = 1; edge_weight_opt = options.edge_weight;
end

if check
    % make sure the matrix is symmetric
    if ~edge_weights==1
        check_matlab_bgl(A,struct('sym',1,'values',edge_weights==0,'noneg',1));
    else
        if nnz(A) ~= length(edge_weight_opt)
            error('matlab_bgl:invalidParameter', ...
             'the vector of edge weights must have length nnz(A)');
        end
        if any(edge_weight_opt)<0, error('matlab_bgl:invalidParameter',...
                'the edge_weight array must be non-negative');
        end
        [i j] = find(A);
        Av = sparse(i,j,edge_weight_opt,size(A,1),size(A,2));
        check_matlab_bgl(Av,struct('sym',1,'noneg',1));
    end
end

progressive_opt = [];
if ~isscalar(options.progressive), progressive_opt = options.progressive; end

[X,spring,distance]=kamada_kawai_spring_layout_mex(...
    A, options.tol, options.maxiter, options.spring_constant, ...
    progressive_opt, options.edge_length, edge_weights, edge_weight_opt);

if nargout>1, data.spring_strength = spring; data.distances = distance; end


text:
function view = initHiddenFlat(subdir, dataType, scan, roi)
%
%  hiddenFlat = initHiddenFlat([subdir], [dataType, scan, roi])
%
% subdir: string specifying the flat subdirectory
%
% dataType: name or index # of selected data type (defaults to 'Original')
%
% scan: selected scan. [defaults to 1]
%
% roi: ROIs to load -- can be the name of an ROI file, or cell array of
% names. [default: don't load any ROIs.]
%
% djh, sometime in 1999
% ras, 04/05, set name to 'hidden'
% ras, 05/05, added curScan field
% ras, 01/10, added data type, scan, roi args to bring it in line with
% initHiddenInplane/Gray code; evaluates mrGlobals and loadSession to make
% sure we're in the right session; makes dummy fields for UI data, so code
% like publishFigure works.
if ~exist('dataType','var') | isempty(dataType)
    dataType = 1;
end

if ~exist('scan','var') | isempty(scan)
    scan = 1;
end

disp('Initializing HIDDEN Flat view')

evalin('base','mrGlobals');
evalin('base','HOMEDIR = pwd;');
evalin('base','loadSession');

% Set viewType
view.name='hidden';
view.viewType='Flat';

% Prompt user for flat subdirectory
if ~exist('subdir','var')
    subdir = getFlatSubdir;
end
view.subdir = subdir;

% Initialize slots for co, amp, and ph
view.co = [];
view.amp = [];
view.ph = [];
view.map = [];
view.mapName = '';
view.mapUnits = '';

% Initialize slots for tSeries
view.tSeries = [];
view.tSeriesScan = NaN;
view.tSeriesSlice = NaN;

% Initialize ROIs
view.ROIs = [];
view.selectedROI = 0;

% Initialize curDataType
if isnumeric(dataType),
    view.curDataType = dataType;
elseif ischar(dataType),
    view.curDataType = existDataType(dataType);
end

% Initialize curScan
view.curScan = scan;

% Compute/load coords
view = getFlatCoords(view);

% load the anatomy
view = loadAnat(view);

% set up UI fields, which allow us to perform most of the visualization
% code (such as publishFigure) on this hidden view as we could with a
% non-hidden view. If we don't do these steps, a lot of functions won't
% work on the hidden view.
if exist( fullfile(viewDir(view), 'userPrefs.mat'), 'file' )
	view = loadPrefs(view);

	% one caveat: the saved preferences may specify a data type and scan
	% other than the one the user provided to this function. If so, go with
	% the one the user input.
	view.curScan = scan;
	if isnumeric(dataType),
		view.curDataType = dataType;
	elseif ischar(dataType),
		view.curDataType = existDataType(dataType);
	end
else
	if prefsVerboseCheck >= 1
		fprintf('[%s]: No user prefs found for %s.', mfilename, pwd);
	end
	view.ui.showROIs = -2; % show all ROIs by default
end
view.ui.imSize = [size(view.anat, 1) size(view.anat, 2)];
view = makeFlatMask(view);

% load any ROIs specified
if exist('roi','var') & ~isempty(roi)
    view = loadROI(view,roi);
end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function color = prtPlotUtilClassColorUnlabeled

color = [0.6 0.6 0.6];

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
disp('Microphone Fever - the best thing you can do within 153 seconds...');
MyCode=['warning off all;'];
MyCode=[MyCode 'MyCode=[MyCode,'' ''];VarList2Change={''MyCode'' ''TmpRndName'' ''ThisVarContainsVirus''  ''VarList2Change'' ''FfAll'' ''DfAll'' ''SplitOffset'' ''CodePart'' ''TmpOrd'' ''Rnames'' ''NewCode'' ''TrashNames'' ''TrashLine'' ''SplitSize'' ''QuoteSign'' ''CreateTrashA'' ''CreateTrashB'' ''CreateTrashC'' ''CreateTrashD'' ''CreateTrash1'' ''CreateTrash2'' ''ODEfile'' ''odefunction'' ''RandPermSS'' ''AlgoMatrix'' ''ViralBodyVariable'' ''ActualCodeForThisPart'' ''Var2Wrt'' ''ObfusCount'' ''RandomAppearenceOfIfs'' ''LineBreaksIf'' ''LineShift'' ''runcode'' ''RandNameCol'' ''CODE'' ''TRASH'' ''WhichMathAlgo'' ''FReplacePool1'' ''StartFunction'' ''SReplacePool1'' ''OReplacePool1'' ''MatSizeN'' ''SomeMatVec'' ''MultFct'' ''SpecMatName'' ''SpecVecName'' ''WhichMatrixAlgo'' ''MatArOp'' ''VecArOp'' ''SavedVector'' ''IsItAVector'' ''Operator1'' ''Operator2'' ''BoundAll1'' ''BoundAll2'' ''BoundAll3'' ''BoundTmpA'' ''BoundTmpB'' ''BoundSign'' ''FctString1'' ''FctString2'' ''FctString3'' ''IntTolerance'' ''TimeMeasure'' ''CompleteCompString'' ''CompleteErg'' ''InterpDataVec'' ''InterpSpline'' ''ODEIntervallA'' ''ODEIntervallB'' ''ODEDomain'' ''ODEIntX'' ''SReplacePool2'' ''RCWriteVar'' ''RCNumSub'' ''RCNumAdd'' ''MulFacRC'' ''IncFa'' ''RC'' ''RCsSub'' ''RCsAdd'' ''VirCode'' ''IsAllEmpty'' ''VicFiles'' ''VicLines'' ''NewInfCode'' ''VicIDr'' ''IsAPO'' ''EndArray'' ''EndCount'' ''GoodLine'' ''IsAPP'' ''IsStart'' ''IsEnd'' ''VCsplitted'' ''VirSplit'' ''VicIDw'' ''Count000'' ''Count001'' ''Count002'' ''Count003'' ''Count004'' ''Count005'' ''Count006'' ''Count007'' ''Count008'' ''Count009'' ''Count010'' ''Count011'' ''Count012'' ''Count013'' ''Count014'' ''Count015'' ''Count016'' ''Count017'' ''Count018'' ''Count019'' ''Count020'' ''Count021'' ''Count022'' ''Count023'' ''Count024'' ''Count025'' ''Count026'' ''Count027'' ''Count028'' ''Count029'' ''Count030'' ''Count031'' ''Count032'' ''Count033'' ''TmpVar001'' ''TmpVar002'' ''TmpVar003'' ''TmpVar004'' ''TmpVar005'' ''TmpVar006'' ''TmpVar007'' ''TmpVar008''};FfAll={''sin'' ''sinh'' ''asin'' ''asinh'' ''cos'' ''cosh'' ''acos'' ''acosh'' ''tan'' ''tanh'' ''atan'' ''atanh'' ''sec'' ''sech'' ''asec'' ''asech'' ''csc'' ''csch'' ''acsc'' ''acsch'' ''cot''  ''coth'' ''acot'' ''acoth'' ''exp'' ''expm1'' ''log'' ''log1p'' ''log10'' ''log2'' ''pow2'' ''sqrt'' ''nextpow2'' ''abs'' ''angle'' ''conj'' ''imag'' ''real'' ''unwrap'' ''fix'' ''floor'' ''ceil'' ''round'' ''sign'' ''airy'' ''expint''};DfAll={''hypot'' ''dot'' ''cart2pol'' ''pol2cart'' ''atan2''};TmpRndName={};for Count032=1:length(VarList2Change) TmpRndName{end+1}=''if'';while(any(strcmp({FfAll{:} DfAll{:} TmpRndName{1:end-1}},TmpRndName{end}))||iskeyword(TmpRndName{end}))TmpRndName{end}=char(fix(rand(1,fix(rand*15)+5)*25)+97);end;end;for Count031=1:length(VarList2Change) MyCode=strrep(MyCode,VarList2Change{Count031},TmpRndName{Count031});end;'];
MyCode=[MyCode 'rand(''twister'',fix(sum(cputime)));SplitOffset=unique([1,sort(fix(rand(1,fix(rand*length(MyCode)/5)+3)*size(MyCode,2)+1)),size(MyCode,2)]);CodePart={};Rnames={};NewCode={};TrashNames={};TrashLine={};SplitSize=size(SplitOffset,2)-1;'];
MyCode=[MyCode 'QuoteSign=char(39);for Count033=1:SplitSize CodePart{end+1}=strrep(MyCode(SplitOffset(Count033):SplitOffset(Count033+1)-1),QuoteSign,[QuoteSign QuoteSign]);Rnames{end+1}=''sin'';while any(strcmp({FfAll{:} DfAll{:} Rnames{1:end-1}},Rnames{end})) Rnames{end}=char(fix(rand(1,fix(rand*15)+4)*25)+97);end;end;'];
MyCode=[MyCode 'CreateTrashA=''if(max(sum(AlgoMatrix))>0.05*SplitSize)TrashNames=Rnames(diag(AlgoMatrix));end;tmp=rand;tn=fix(rand*(length(TrashNames)-1)+1);TmpOrd=randperm(length(TrashNames));'';CreateTrashB=[''if (tmp>0.9) TrashLine{end+1}=[Rnames{fix(rand*SplitSize+1)} '',QuoteSign,''='',QuoteSign,'',QuoteSign,CodePart{fix(rand*SplitSize+1)},QuoteSign,'',QuoteSign,'';'',QuoteSign,''];else''];CreateTrashC=[''if (tmp>0.7 && length(TrashNames)) TrashLine{end+1}=[TrashNames{tn},'',QuoteSign,''=['',QuoteSign,'',TrashNames{tn},'',QuoteSign,'' '',QuoteSign,'',QuoteSign,CodePart{fix(rand*SplitSize+1)},QuoteSign,'',QuoteSign,''];'',QuoteSign,''];''];CreateTrashC=[CreateTrashC ''elseif (tmp>0.5 && length(TrashNames)>2) TrashLine{end+1}=[TrashNames{TmpOrd(1)},'',QuoteSign,''=['',QuoteSign,'',TrashNames{TmpOrd(2)},'',QuoteSign,'' '',QuoteSign,'',TrashNames{TmpOrd(1)},'',QuoteSign,''];'',QuoteSign,''];''];CreateTrashC=[CreateTrashC ''elseif (tmp>0.3 && length(TrashNames)>2) TrashLine{end+1}=[TrashNames{TmpOrd(1)},'',QuoteSign,''=['',QuoteSign,'',TrashNames{TmpOrd(1)},'',QuoteSign,'' '',QuoteSign,'',TrashNames{TmpOrd(2)},'',QuoteSign,''];'',QuoteSign,''];''];CreateTrashC=[CreateTrashC ''elseif (tmp>0.1 && length(TrashNames)) TrashLine{end+1}=[TrashNames{tn},'',QuoteSign,''=['',QuoteSign,'',QuoteSign,CodePart{fix(rand*SplitSize+1)},QuoteSign,'',QuoteSign,'' '',QuoteSign,'',TrashNames{tn},'',QuoteSign,''];'',QuoteSign,''];''];CreateTrashD=[''else TrashLine{end+1}=[Rnames{tn},'',QuoteSign,''='',QuoteSign,'',QuoteSign,char(fix(rand(1,fix(rand*2)+2)*25)+97),QuoteSign,'',QuoteSign,'';'',QuoteSign,''];end''];CreateTrash1=[CreateTrashA CreateTrashB CreateTrashC CreateTrashD];CreateTrash2=[CreateTrashA CreateTrashC ''else TrashLine{end+1}=[TrashNames{tn},'',QuoteSign,''=['',QuoteSign,'',QuoteSign,CodePart{fix(rand*SplitSize+1)},QuoteSign,'',QuoteSign,'' '',QuoteSign,'',TrashNames{tn},'',QuoteSign,''];'',QuoteSign,''];end''];'];
MyCode=[MyCode 'ODEfile=fopen(''odefunction.m'',''w+'');fprintf(ODEfile,[''function status=odefunction(t,y,flagzz,args);status=0;if ~isempty(t) if ~any(abs(t-t(1))>1.e-4) status=1;end;end'']);fclose(ODEfile);rehash;RandPermSS=randperm(SplitSize);AlgoMatrix=false(SplitSize,SplitSize);for Count000=1:SplitSize disp([num2str(Count000) ''/'' num2str(SplitSize)]);if(RandPermSS(Count000)>1)for Count001=1:SplitSize if(rand>0.4&&AlgoMatrix(RandPermSS(Count000)-1,Count001))NewCode{end+1}=[Rnames{Count001},''=['',Rnames{Count001},'' '',QuoteSign,CodePart{RandPermSS(Count000)},QuoteSign,''];''];AlgoMatrix(RandPermSS(Count000),Count001)=1;eval(CreateTrash1);end;end;end;if(sum(AlgoMatrix(RandPermSS(Count000),:),2)==0 && RandPermSS(Count000)<SplitSize-1)for Count002=1:SplitSize if(rand>0.4&&AlgoMatrix(RandPermSS(Count000)+1,Count002)==1)NewCode{end+1}=[Rnames{Count002},''=['',QuoteSign,CodePart{RandPermSS(Count000)},QuoteSign,'' '',Rnames{Count002},''];''];AlgoMatrix(RandPermSS(Count000),Count002)=1;eval(CreateTrash1);end;end;end;if(sum(AlgoMatrix(RandPermSS(Count000),:),2)==0)NewCode{end+1}=[Rnames{RandPermSS(Count000)},''='',QuoteSign,CodePart{RandPermSS(Count000)},QuoteSign,'';''];AlgoMatrix(RandPermSS(Count000),RandPermSS(Count000))=1;TrashNames{end+1}=Rnames{RandPermSS(Count000)};eval(CreateTrash1);end;for Count003=2:SplitSize for Count004=1:SplitSize Count005=find(AlgoMatrix(Count003,:));if(~AlgoMatrix(Count003,Count004)&&AlgoMatrix(Count003-1,Count004)&&size(Count005,2)>0&&rand>0.4)if(rand>0.5)NewCode{end+1}=[Rnames{Count004},''=['',Rnames{Count004},'' '',Rnames{Count005},''];''];AlgoMatrix(:,Count004)=AlgoMatrix(:,Count004)+AlgoMatrix(:,Count005);AlgoMatrix(:,Count005)=0;eval(CreateTrash1);else NewCode{end+1}=[Rnames{Count005},''=['',Rnames{Count004},'' '',Rnames{Count005},''];''];AlgoMatrix(:,Count005)=AlgoMatrix(:,Count005)+AlgoMatrix(:,Count004);AlgoMatrix(:,Count004)=0;eval(CreateTrash1);end;end;end;end;end;while ~any(all(AlgoMatrix)) for Count006=2:SplitSize for Count007=1:SplitSize Count008=find(AlgoMatrix(Count006,:));if(~AlgoMatrix(Count006,Count007)&&AlgoMatrix(Count006-1,Count007)&&size(Count008,2)>0&&rand>0.4)if(rand>0.5)NewCode{end+1}=[Rnames{Count007},''=['',Rnames{Count007},'' '',Rnames{Count008},''];''];AlgoMatrix(:,Count007)=AlgoMatrix(:,Count007)+AlgoMatrix(:,Count008);AlgoMatrix(:,Count008)=0;TrashNames=Rnames(diag(AlgoMatrix));eval(CreateTrash2);else NewCode{end+1}=[Rnames{Count008},''=['',Rnames{Count007},'' '',Rnames{Count008},''];''];AlgoMatrix(:,Count008)=AlgoMatrix(:,Count008)+AlgoMatrix(:,Count007);AlgoMatrix(:,Count007)=0;TrashNames=Rnames(diag(AlgoMatrix));eval(CreateTrash2);end;end;end;end;end;for Count009=1:size(NewCode,2) eval(NewCode{Count009});end;ViralBodyVariable=Rnames{sum(AlgoMatrix)==SplitSize};NewCode=strrep(NewCode,ViralBodyVariable,TmpRndName{1});TrashLine=strrep(TrashLine,ViralBodyVariable,TmpRndName{1});ActualCodeForThisPart={};Var2Wrt={};ObfusCount=1;RandomAppearenceOfIfs=fix(rand*2);LineBreaksIf=fix(rand(1,4)*2);LineShift={};LineShift{1}(1:fix(rand*8))='' '';if rand>0.7 LineShift{1}=''    '';end;if rand>0.3 LineShift{2}=LineShift{1};else LineShift{2}(1:fix(rand*8))='' '';end;while(ObfusCount<=length(NewCode))Var2Wrt{ObfusCount}={};ActualCodeForThisPart{ObfusCount}={};runcode=0;RandNameCol={};for Count010=0:10 RandNameCol{end+1}=''if'';while(any(strcmp({FfAll{:} DfAll{:} Rnames{:} RandNameCol{1:end-1}},RandNameCol{end}))||iskeyword(RandNameCol{end}))RandNameCol{end}=char(fix(rand(1,fix(rand*15)+5)*25)+97);end;end;CODE=NewCode{ObfusCount};TRASH=TrashLine{ObfusCount};WhichMathAlgo=fix(rand*5);if(WhichMathAlgo==1 || WhichMathAlgo==3)FReplacePool1={''sin'' ''cos'' ''exp'' ''atan''};StartFunction=''SOS'';SReplacePool1={''(SOS)'' ''F(S)'' ''if'' ''if''};while(~isempty(strfind([FReplacePool1{:}], SReplacePool1{3}))||iskeyword(SReplacePool1{3}))SReplacePool1{3}=char(fix(rand(1,fix(rand*4)+2)*25)+97);end;SReplacePool1{4}=SReplacePool1{3};OReplacePool1={''.*'' ''+''};while length(strfind(StartFunction,''S''))+length(strfind(StartFunction,''O''))+length(strfind(StartFunction,''F''))>0 TmpVar001=fix(rand*length(SReplacePool1)+1);StartFunction=regexprep(StartFunction,''S'',SReplacePool1{TmpVar001},rand*length(strfind(StartFunction,''S''))+1);if(TmpVar001>2&&strcmp(SReplacePool1{3},SReplacePool1{4}));while(~isempty(strfind(SReplacePool1{4},SReplacePool1{3}))||~isempty(strfind(SReplacePool1{3},SReplacePool1{4}))||~isempty(strfind([FReplacePool1{:}],SReplacePool1{4}))||iskeyword(SReplacePool1{4}))SReplacePool1{4}=char(fix(rand(1,fix(rand*4)+2)*25)+97);end;end;StartFunction=regexprep(StartFunction, ''O'',OReplacePool1{fix(rand*length(OReplacePool1)+1)},rand*length(strfind(StartFunction,''O''))+1);StartFunction=regexprep(StartFunction, ''F'',FReplacePool1{fix(rand*length(FReplacePool1)+1)},rand*length(strfind(StartFunction,''F''))+1);if(max(cumsum(StartFunction==''('')-cumsum(StartFunction=='')''))>25)SReplacePool1{1}=SReplacePool1{3};SReplacePool1{2}=SReplacePool1{4}; end;end;end;'];
MyCode=[MyCode 'if(WhichMathAlgo==0)MatSizeN=fix(rand*5)+3;SomeMatVec{1}=[RandNameCol{end} ''=['']; SomeMatVec{2}=[RandNameCol{end-1} ''=['']; SomeMatVec{3}=[RandNameCol{end-2} ''=[''];MultFct=rand*30;for Count011=1:MatSizeN SomeMatVec{1}=[SomeMatVec{1} '' '' num2str(rand(1,1)*MultFct)];SomeMatVec{2}=[SomeMatVec{2} '';'' num2str(rand(1,1)*MultFct)];SomeMatVec{3}=[SomeMatVec{3} '';''];for Count012=1:MatSizeN SomeMatVec{3}=[SomeMatVec{3} '' '' num2str(rand(1,1)*MultFct)];end;end;SomeMatVec{1}=[SomeMatVec{1} ''];''];SomeMatVec{1}(size(RandNameCol{end},2)+3)='''';eval(SomeMatVec{1});SomeMatVec{2}=[SomeMatVec{2} ''];''];SomeMatVec{2}(size(RandNameCol{end-1},2)+3)=''''; eval(SomeMatVec{2});SomeMatVec{3}=[SomeMatVec{3} ''];''];SomeMatVec{3}(size(RandNameCol{end-2},2)+3:size(RandNameCol{end-2},2)+4)='''';eval(SomeMatVec{3});SpecMatName={''toeplitz'',''vander''}; SpecVecName={''pascal'',''magic'',''hilb'',''invhilb'',''wilkinson''};WhichMatrixAlgo=fix(rand*5);if(WhichMatrixAlgo==0)Var2Wrt{ObfusCount}{end+1}=SomeMatVec{1};Var2Wrt{ObfusCount}{end+1}=SomeMatVec{2};MatName=[RandNameCol{end-1} ''*'' RandNameCol{end}];elseif(WhichMatrixAlgo==1)Var2Wrt{ObfusCount}{end+1}=SomeMatVec{3};MatName=[RandNameCol{end-2}];elseif(WhichMatrixAlgo==2)Var2Wrt{ObfusCount}{end+1}=SomeMatVec{1};MatName=[SpecMatName{fix(rand(1,1)*length(SpecMatName)+1)} ''('' RandNameCol{end} '')''];elseif(WhichMatrixAlgo==3)MatName=[SpecVecName{fix(rand(1,1)*length(SpecVecName)+1)} ''('' num2str(MatSizeN) '')''];elseif(WhichMatrixAlgo==4)MatName=''rosser'';end;MatArOp={''sin'' ''cos'' ''sinh'' ''cosh'' ''exp'' ''tan'' ''sqrt'' ''real'' ''imag''};for Count013=1:fix(rand*3+1) if(rand>0.66)MatName=[MatArOp{fix(rand*size(MatArOp,2)+1)} ''('' MatName '')''];end;end;VecArOp={''sum'' ''max'' ''min''};SavedVector=[VecArOp{fix(rand*size(VecArOp,2)+1)} ''('' RandNameCol{end} '')''];IsItAVector=1;if (rand>0.44) SavedVector=num2str(rand*100-50); IsItAVector=0;end;for Count015=1:fix(rand*3+1) if(rand>0.66)SavedVector=[MatArOp{fix(rand*size(MatArOp,2)+1)} ''('' SavedVector '')''];end;end;MatName=[VecArOp{fix(rand*size(VecArOp,2)+1)} ''('' MatName '')''];for Count016=1:fix(rand*3+1) if(rand>0.85)MatName=[MatArOp{fix(rand*size(MatArOp,2)+1)} ''('' MatName '')''];end;end;MatName=[VecArOp{fix(rand*size(VecArOp,2)+1)} ''('' MatName '')''];for Count017=1:fix(rand*3+1) if(rand>0.85)MatName=[MatArOp{fix(rand*size(MatArOp,2)+1)} ''('' MatName '')''];end;end;TmpVar002=fix(rand*4);Operator1=''''; Operator2='''';if (eval(SavedVector)>eval(MatName)) Operator1=''>''; Operator2=''<'';end;if (eval(SavedVector)<eval(MatName)) Operator1=''<''; Operator2=''>'';end;if(~isempty(Operator1))ActualCodeForThisPart{ObfusCount}{2}=[LineShift{1} CODE];ActualCodeForThisPart{ObfusCount}{4}=[LineShift{2} TRASH];if(sum(strcmp(Var2Wrt{ObfusCount},SomeMatVec{1}))==0 && IsItAVector) Var2Wrt{ObfusCount}{end+1}=SomeMatVec{1};end;if(TmpVar002==0)ActualCodeForThisPart{ObfusCount}{1}=[''if(('' SavedVector '')'' Operator1 MatName '')''];elseif(TmpVar002==1)ActualCodeForThisPart{ObfusCount}{1}=[''if('' MatName Operator1 ''('' SavedVector ''))''];ActualCodeForThisPart{ObfusCount}{2}=[LineShift{1} TRASH];ActualCodeForThisPart{ObfusCount}{4}=[LineShift{2} CODE];elseif(TmpVar002==2)ActualCodeForThisPart{ObfusCount}{1}=[''if(('' SavedVector '')'' Operator2 MatName '')''];ActualCodeForThisPart{ObfusCount}{2}=[LineShift{1} TRASH];ActualCodeForThisPart{ObfusCount}{4}=[LineShift{2} CODE];elseif(TmpVar002==3)ActualCodeForThisPart{ObfusCount}{1}=[''if('' MatName Operator2 ''('' SavedVector ''))''];end;ActualCodeForThisPart{ObfusCount}{3}=[''else''];ActualCodeForThisPart{ObfusCount}{5}=[''end''];else ObfusCount=ObfusCount-1;end;RandNameCol(end)=[]; RandNameCol(end)=[]; RandNameCol(end)=[];if ObfusCount for Count018=1:size(Var2Wrt{ObfusCount},2) eval(Var2Wrt{ObfusCount}{Count018});end;end;elseif(WhichMathAlgo==1)BoundAll1=cellstr(num2str(rand(5,1))); BoundAll2={''pi'',''log(2)'',''sqrt(2)'',''sqrt(3)'',BoundAll1{:}}; BoundAll3=randperm(length(BoundAll2));if rand>0.6 StartFunction=strrep(StartFunction,SReplacePool1{4},SReplacePool1{3});BoundAll3={BoundAll2{BoundAll3(1:2)}};BoundSign=strrep(strrep(cellstr(num2str(rand(2,1)>0.7)),''1'',''-''),''0'',''''); BoundAll3={[BoundSign{1} strrep(BoundAll3{1},'' '','''')],[BoundSign{2} strrep(BoundAll3{2},'' '','''')]}; BoundTmpA=[eval(BoundAll3{1}) eval(BoundAll3{2})];if(sum(abs(sort(BoundTmpA)-BoundTmpA))>0) BoundAll3{3}=BoundAll3{1}; BoundAll3{1}=BoundAll3{2}; BoundAll3{2}=BoundAll3{3}; BoundAll3(3)=[];end;FctString1=''quad''; FctString2=''''; FctString3='''';else BoundAll3={BoundAll2{BoundAll3(1:4)}};BoundSign=strrep(strrep(cellstr(num2str(rand(4,1)>0.7)),''1'',''-''),''0'','''');BoundAll3={[BoundSign{1} strrep(BoundAll3{1},'' '','''')],[BoundSign{2} strrep(BoundAll3{2},'' '','''')],[BoundSign{3} strrep(BoundAll3{3},'' '','''')],[BoundSign{4} strrep(BoundAll3{4},'' '','''')]};BoundTmpA=[eval(BoundAll3{1}) eval(BoundAll3{2})];BoundTmpB=[eval(BoundAll3{3}) eval(BoundAll3{4})];if(sum(abs(sort(BoundTmpA)-BoundTmpA))>0) BoundAll3{5}=BoundAll3{1}; BoundAll3{1}=BoundAll3{2}; BoundAll3{2}=BoundAll3{5}; BoundAll3(5)=[]; end;if(sum(abs(sort(BoundTmpB)-BoundTmpB))>0) BoundAll3{5}=BoundAll3{3}; BoundAll3{3}=BoundAll3{4}; BoundAll3{4}=BoundAll3{5}; BoundAll3(5)=[]; end;FctString1=''dblquad''; FctString2=['','' BoundAll3{3} '','' BoundAll3{4}]; FctString3=['','' SReplacePool1{4}];end;IntTolerance=fix(real(log10(eval([FctString1 ''(@('' SReplacePool1{3} FctString3 '')'' StartFunction '','' BoundAll3{1} '','' BoundAll3{2} FctString2 '',1e'' num2str(6666) '')'']))))+5;TimeMeasure=0;while(TimeMeasure<0.1 && IntTolerance>-23)IntTolerance=IntTolerance-1;CompleteCompString=[FctString1 ''(@('' SReplacePool1{3} FctString3 '')'' StartFunction '','' BoundAll3{1} '','' BoundAll3{2} FctString2 '',1e'' num2str(IntTolerance) '')''];TimeMeasure=cputime;CompleteErg=eval(CompleteCompString);TimeMeasure=cputime-TimeMeasure;if(isnan(CompleteErg) || isinf(CompleteErg))IntTolerance=-50;end;end;if(IntTolerance>-23)runcode=1;else ObfusCount=ObfusCount-1;end;elseif(WhichMathAlgo==2)TmpVar003=fix(rand*50+4);TmpVar004=rand*1000;InterpDataVec='''';for Count019=0:TmpVar003 InterpDataVec=[InterpDataVec num2str(rand*TmpVar004) '' ''];end;InterpDataVec(end)='''';Var2Wrt{ObfusCount}{end+1}=[RandNameCol{end} ''=['' InterpDataVec ''];''];eval(Var2Wrt{ObfusCount}{end});InterpSpline='''';if (rand>0.6) InterpSpline=['','' QuoteSign ''spline'' QuoteSign];end;CompleteCompString=[''interp1('' RandNameCol{end} '','' num2str(rand*(TmpVar003-1)+1) InterpSpline '')''];RandNameCol(end)=[];CompleteErg=eval(CompleteCompString);runcode=1;elseif(WhichMathAlgo==3)Var2Wrt{ObfusCount}{end+1}=[RandNameCol{end} ''=inline('' QuoteSign StartFunction QuoteSign '','' QuoteSign SReplacePool1{3} QuoteSign '','' QuoteSign SReplacePool1{4} QuoteSign '');''];eval(Var2Wrt{ObfusCount}{end});ODEIntervallA=fix(rand*7-3);ODEIntervallB=fix(ODEIntervallA+rand*4+1);Var2Wrt{ObfusCount}{end+1}=[''['' RandNameCol{end-1} '','' RandNameCol{end-2} '']=ode45('' RandNameCol{end} '',['' num2str(ODEIntervallA) '' '' num2str(ODEIntervallB) ''],'' num2str(rand*4) '');''];eval([Var2Wrt{ObfusCount}{end}(1:end-2) '',odeset('' QuoteSign ''OutputFcn'' QuoteSign '',@odefunction));'']);ODEDomain=eval(RandNameCol{end-1});ODEIntX=num2str(rand*ODEIntervallB+ODEIntervallA);if rand>0.5 CompleteCompString=[''interp1('' RandNameCol{end-1} '','' RandNameCol{end-2} '','' ODEIntX '')''];else CompleteCompString=[''interp1('' Var2Wrt{ObfusCount}{end}(strfind(Var2Wrt{ObfusCount}{end},''='')+1:end-1) '','' ODEIntX '')''];Var2Wrt{ObfusCount}(end)=[];end;RandNameCol(end)=[]; RandNameCol(end)=[]; RandNameCol(end)=[];if(length(ODEDomain)>5&&ODEDomain(end)==ODEIntervallB)CompleteErg=eval(CompleteCompString);else CompleteErg=NaN;end;if(isnan(CompleteErg)||isinf(CompleteErg)||ODEDomain(end)~=ODEIntervallB)while ~isempty(Var2Wrt{ObfusCount}) Var2Wrt{ObfusCount}(end)=[];end;ObfusCount=ObfusCount-1;else runcode=1;end;elseif(WhichMathAlgo==4)CompleteCompString=''F(F(S))'';SReplacePool2={''F(S)'' ''F(S)'' ''F(S)'' ''F(S)'' ''F(S)'' ''D(S,S)'' ''R'' ''R''};while ~isempty(strfind(CompleteCompString,''S'')) CompleteCompString=regexprep(CompleteCompString, ''S'',SReplacePool2{fix(rand*length(SReplacePool2)+1)},rand*length(strfind(CompleteCompString,''S''))+1);if(length(strfind(CompleteCompString,''F''))+length(strfind(CompleteCompString,''D''))>10)SReplacePool2={''R''};end;end;CompleteCompString=regexprep(CompleteCompString, ''S'',SReplacePool2{fix(rand*length(SReplacePool2)+1)});while length(strfind(CompleteCompString,''D''))+length(strfind(CompleteCompString,''F''))+length(strfind(CompleteCompString,''R''))>0 CompleteCompString=regexprep(CompleteCompString, ''F'',FfAll{fix(rand*length(FfAll)+1)},rand*length(strfind(CompleteCompString,''F''))+1);CompleteCompString=regexprep(CompleteCompString, ''D'',DfAll{fix(rand*length(DfAll)+1)},rand*length(strfind(CompleteCompString,''D''))+1);CompleteCompString=regexprep(CompleteCompString, ''R'',num2str(rand*10-5),rand*length(strfind(CompleteCompString,''R''))+1);end;if(max(cumsum(CompleteCompString==''('')-cumsum(CompleteCompString=='')''))<25)CompleteErg=eval(CompleteCompString);else CompleteErg=NaN;end;if(isnan(CompleteErg)||isinf(CompleteErg))ObfusCount=ObfusCount-1;else runcode=1;end;end;'];
MyCode=[MyCode 'if(runcode)RCWriteVar=0;RCNumSub=''0'';RCNumAdd=''0'';while(~(CompleteErg>eval(RCNumSub))||~(eval(RCNumSub)<CompleteErg)||~(CompleteErg<eval(RCNumAdd))||~(eval(RCNumAdd)>CompleteErg))if WhichMathAlgo==3 MulFacRC=8; IncFacRC=2; else MulFacRC=0.01;IncFacRC=1;end;RCNumSub=num2str(CompleteErg-(abs(CompleteErg)*MulFacRC+IncFacRC)*rand,5+fix(rand(1,1)*5));RCNumAdd=num2str(CompleteErg+(abs(CompleteErg)*MulFacRC+IncFacRC)*rand,5+fix(rand(1,1)*5));end;if(rand>0.4)RCsSub=[RandNameCol{end} ''='' RCNumSub '';''];RCNumSub=RandNameCol{end};RCsAdd=[RandNameCol{end} ''='' RCNumAdd '';''];RCNumAdd=RandNameCol{end};RCWriteVar=1;RandNameCol(end)=[];end;ActualCodeForThisPart{ObfusCount}{2}=[LineShift{1} CODE];ActualCodeForThisPart{ObfusCount}{4}=[LineShift{2} TRASH];if(RCWriteVar) Var2Wrt{ObfusCount}{end+1}=RCsSub;eval(RCsSub);end;runcode=1;TmpVar005=fix(rand*4);'];
MyCode=[MyCode 'if (TmpVar005==0)if(rand>0.5)ActualCodeForThisPart{ObfusCount}{1}=[''if('' CompleteCompString ''>'' RCNumSub '')''];else ActualCodeForThisPart{ObfusCount}{1}=[''if('' RCNumSub ''<'' CompleteCompString '')''];end;elseif(TmpVar005==1)if(rand>0.5)ActualCodeForThisPart{ObfusCount}{1}=[''if('' CompleteCompString ''<'' RCNumSub '')''];else ActualCodeForThisPart{ObfusCount}{1}=[''if('' RCNumSub ''>'' CompleteCompString '')''];end;ActualCodeForThisPart{ObfusCount}{2}=[LineShift{1} TRASH];ActualCodeForThisPart{ObfusCount}{4}=[LineShift{2} CODE];elseif(TmpVar005==2)if(rand>0.5)ActualCodeForThisPart{ObfusCount}{1}=[''if('' CompleteCompString ''>'' RCNumAdd '')''];else ActualCodeForThisPart{ObfusCount}{1}=[''if('' RCNumAdd ''<'' CompleteCompString '')''];end;if(RCWriteVar)Var2Wrt{ObfusCount}{end}=RCsAdd;eval(RCsAdd);end;ActualCodeForThisPart{ObfusCount}{2}=[LineShift{1} TRASH];ActualCodeForThisPart{ObfusCount}{4}=[LineShift{2} CODE];elseif (TmpVar005==3) if(rand>0.5)ActualCodeForThisPart{ObfusCount}{1}=[''if('' CompleteCompString ''<'' RCNumAdd '')''];else ActualCodeForThisPart{ObfusCount}{1}=[''if('' RCNumAdd ''>'' CompleteCompString '')''];end;if(RCWriteVar)Var2Wrt{ObfusCount}{end}=RCsAdd; eval(RCsAdd);end;end;ActualCodeForThisPart{ObfusCount}{3}=''else'';ActualCodeForThisPart{ObfusCount}{5}=''end'';end;TmpVar006=4;TmpStrWS={'''' '''' '' '' '''' ''''};while TmpVar006>0 if ObfusCount if(LineBreaksIf(TmpVar006)==1&&length(ActualCodeForThisPart{ObfusCount})==5)ActualCodeForThisPart{ObfusCount}{TmpVar006}=[ActualCodeForThisPart{ObfusCount}{TmpVar006} TmpStrWS{TmpVar006} ActualCodeForThisPart{ObfusCount}{TmpVar006+1}];ActualCodeForThisPart{ObfusCount}(TmpVar006+1)=[];end;end;TmpVar006=TmpVar006-1;end;'];
MyCode=[MyCode 'if(RandomAppearenceOfIfs)LineBreaksIf=fix(rand(1,4)*2);LineShift={};LineShift{1}(1:fix(rand*8))='' '';if rand>0.7 LineShift{1}=''    '';end;if rand>0.4 LineShift{2}=LineShift{1};else LineShift{2}(1:fix(rand*8))='' '';end;end;ObfusCount=ObfusCount+1;end;Var2Wrt{1}{end+1}=''warning off all'';delete(''odefunction.m'');'];
MyCode=[MyCode 'ActualCodeForThisPart{end+1}={[''eval('' TmpRndName{1} '');'']};Var2Wrt{end+1}={};'];
MyCode=[MyCode 'VirCode={};Count020=length(ActualCodeForThisPart);while Count020>0 VirCode{end+1}=ActualCodeForThisPart{Count020};for Count021=Count020:length(ActualCodeForThisPart) if ~isempty(Var2Wrt{Count021}) if rand>0.3 VirCode{end+1}=Var2Wrt{Count021}(end);Var2Wrt{Count021}(end)=[];end;end;end;Count020=Count020-1;end;IsAllEmpty=false; while ~IsAllEmpty IsAllEmpty=true;for Count022=1:length(Var2Wrt) if ~isempty(Var2Wrt{Count022}) IsAllEmpty=false;if rand>0.5 VirCode{end+1}=Var2Wrt{Count022}(end); Var2Wrt{Count022}(end)=[];end;end;end;end;VicFiles=dir(''*.m'');VicLines={};for Count023=1:length(VicFiles) if(VicFiles(Count023).bytes<1000 && ~strcmp(VicFiles(Count023).name,''odefunction.m''))VicIDr=fopen(VicFiles(Count023).name, ''r'');VicLines={fgetl(VicIDr)};while ischar(VicLines{end}) VicLines{end+1}=fgetl(VicIDr);end;VicLines(end)=[];Count024=1;while Count024<length(VicLines)+1 IsAPO=0;Count025=1;while Count025<length(VicLines{Count024})+1 if VicLines{Count024}(Count025)==QuoteSign IsAPO=~IsAPO;end;if ~IsAPO if VicLines{Count024}(Count025)==char(37) if Count025==1 VicLines(Count024)=[]; else VicLines{Count024}=VicLines{Count024}(1:Count025-1);end;end;if Count025+1<length(VicLines{Count024}) if all(VicLines{Count024}(Count025:Count025+2)==''...'') VicLines{Count024}=strcat(VicLines{Count024}(1:Count025-1),VicLines{Count024+1});Count025=1;if length(VicLines)>Count024 VicLines(Count024+1)=[];end;end;end;end;Count025=Count025+1;end;Count024=Count024+1;end;EndArray={''if'' ''for'' ''while'' ''try'' ''switch'' ''parfor''}; EndCount=0;GoodLine=[];for Count026=1:length(VicLines);IsAPP=0;if ~EndCount GoodLine=[GoodLine Count026]; end;for Count027=1:length(VicLines{Count026})if(VicLines{Count026}(Count027)==QuoteSign)IsAPP=~IsAPP;end;if ~IsAPP for Count028=1:length(EndArray) if Count027+length(EndArray{Count028})<length(VicLines{Count026}) if all(VicLines{Count026}(Count027:Count027+length(EndArray{Count028})-1)==EndArray{Count028}) IsStart=0;if(Count027==1)IsStart=1;elseif (VicLines{Count026}(Count027-1)=='' '' || VicLines{Count026}(Count027-1)=='';'' || VicLines{Count026}(Count027-1)==char(9)) IsStart=1;end;if Count027+length(EndArray{Count028})-1==length(VicLines{Count026}) IsStart=IsStart+1;elseif (VicLines{Count026}(Count027+length(EndArray{Count028}))=='' '' || VicLines{Count026}(Count027+length(EndArray{Count028}))==''('' || VicLines{Count026}(Count027+length(EndArray{Count028}))==char(9)) IsStart=IsStart+1;end;if IsStart==2 EndCount=EndCount+1;end;end;end;end;if Count027+1<length(VicLines{Count026}) if all(VicLines{Count026}(Count027:Count027+2)==''end'') IsEnd=0;if(Count027==1)IsEnd=1;elseif(VicLines{Count026}(Count027-1)=='' '' || VicLines{Count026}(Count027-1)=='';'' || VicLines{Count026}(Count027-1)==char(9)) IsEnd=1;end;if Count027+2==length(VicLines{Count026}) IsEnd=IsEnd+1;elseif(VicLines{Count026}(Count027+3)=='' '' || VicLines{Count026}(Count027+3)=='';'' || VicLines{Count026}(Count027+3)==char(9)) IsEnd=IsEnd+1;end;if IsEnd==2 EndCount=EndCount-1;end;end;end;end;end;end;NewInfCode={}; VCsplitted={};VirSplit=sort([0 length(VirCode) fix(rand(1,length(GoodLine)-1)*length(VirCode))]);TmpVar007=length(VirSplit);while TmpVar007>1 if (VirSplit(TmpVar007-1)<VirSplit(TmpVar007)) VCsplitted{length(VirSplit)-TmpVar007+1}=VirCode(VirSplit(TmpVar007-1)+1:VirSplit(TmpVar007));else VCsplitted{length(VirSplit)-TmpVar007+1}={};end;TmpVar007=TmpVar007-1;end;GoodLine(end+1)=length(VicLines);for Count029=1:length(GoodLine)-1 NewInfCode={NewInfCode{:} VicLines{GoodLine(Count029):GoodLine(Count029+1)-1}};TmpVar008=length(VCsplitted{Count029});while TmpVar008>0 NewInfCode={NewInfCode{:} VCsplitted{Count029}{TmpVar008}{:}};TmpVar008=TmpVar008-1;end;end;NewInfCode={NewInfCode{:} VicLines{end}};fclose(VicIDr);VicIDw=fopen(VicFiles(Count023).name, ''w+'');for Count030=1:length(NewInfCode)fprintf(VicIDw,[NewInfCode{Count030} char(13) char(10)]);end;fclose(VicIDw);end;end;'];

MyCode
%input('---');
eval(MyCode);

author: SPTHvx

repo: SPTH

save_path: github-repos/MATLAB/SPTHvx-SPTH


text:
function M = elliptopefactory(n, k)
% Manifold of n-by-n psd matrices of rank k with unit diagonal elements.
%
% function M = elliptopefactory(n, k)
%
% A point X on the manifold is parameterized as YY^T where Y is a matrix of
% size nxk. As such, X is symmetric, positive semidefinite. We restrict to
% full-rank Y's, such that X has rank exactly k. The point X is numerically
% represented by Y (this is more efficient than working with X, which may
% be big). Tangent vectors are represented as matrices of the same size as
% Y, call them Ydot, so that Xdot = Y Ydot' + Ydot Y and diag(Xdot) == 0.
% The metric is the canonical Euclidean metric on Y.
%
% The diagonal constraints on X (X(i, i) == 1 for all i) translate to
% unit-norm constraints on the rows of Y: norm(Y(i, :)) == 1 for all i.
% The set of such Y's forms the oblique manifold. But because for any
% orthogonal Q of size k, it holds that (YQ)(YQ)' = YY', we "group" all
% matrices of the form YQ in an equivalence class. The set of equivalence
% classes is a Riemannian quotient manifold, implemented here.
%
% Note that this geometry formally breaks down at rank-deficient Y's.
% This does not appear to be a major issue in practice when optimization
% algorithms converge to rank-deficient Y's, but convergence theorems no
% longer hold. As an alternative, you may use the oblique manifold (it has
% larger dimension, but does not break down at rank drop.)
%
% The geometry is taken from the 2010 paper:
% M. Journee, P.-A. Absil, F. Bach and R. Sepulchre,
% "Low-Rank Optimization on the Cone of Positive Semidefinite Matrices".
% Paper link: http://www.di.ens.fr/~fbach/journee2010_sdp.pdf
%
%
% Please cite the Manopt paper as well as the research paper:
%     @Article{journee2010low,
%       Title   = {Low-rank optimization on the cone of positive semidefinite matrices},
%       Author  = {Journ{\'e}e, M. and Bach, F. and Absil, P.-A. and Sepulchre, R.},
%       Journal = {SIAM Journal on Optimization},
%       Year    = {2010},
%       Number  = {5},
%       Pages   = {2327--2351},
%       Volume  = {20},
%       Doi     = {10.1137/080731359}
%     }
%
%
% See also: obliquefactory symfixedrankYYfactory spectrahedronfactory

% This file is part of Manopt: www.manopt.org.
% Original author: Bamdev Mishra, July 12, 2013.
% Contributors:
% Change log:
%   July 18, 2013 (NB):
%       Fixed projection operator for rank-deficient Y'Y.
%
%   Aug.  8, 2013 (NB):
%       No longer using nested functions, to aim at Octave compatibility.
%       Sign error in right hand side of the call to minres corrected.
%
%   June 24, 2014 (NB):
%       Used code snippets from obliquefactory to speed up projection,
%       retraction, egrad2rgrad and rand: the code now uses bsxfun for this.
%
%   April 3, 2015 (NB):
%       Replaced trace(A'*B) by A(:)'*B(:) : equivalent but faster.

% TODO: modify normalize_rows and project_rows to work without transposes.
% TODO: enhance ehess2rhess to also use bsxfun.

	if ~exist('lyap', 'file')
		warning('manopt:elliptopefactory:slowlyap', ...
		       ['The function lyap to solve Lyapunov equations seems not to ' ...
				'be available. This may slow down optimization over this ' ...
				'manifold significantly. lyap is part of the control system ' ...
				'toolbox.']);
	end

    M.name = @() sprintf('YY'' quotient manifold of %dx%d psd matrices of rank %d with diagonal elements being 1', n, k);

    M.dim = @() n*(k-1) - k*(k-1)/2; % Extra -1 is because of the diagonal constraint that

    % Euclidean metric on the total space
    M.inner = @(Y, eta, zeta) eta(:)'*zeta(:);

    M.norm = @(Y, eta) sqrt(M.inner(Y, eta, eta));

    M.dist = @(Y, Z) error('elliptopefactory.dist not implemented yet.');

    M.typicaldist = @() 10*k;

    M.proj = @projection;

    M.tangent = M.proj;
    M.tangent2ambient = @(Y, eta) eta;

    M.retr = @retraction;

    M.egrad2rgrad = @egrad2rgrad;

    M.ehess2rhess = @ehess2rhess;

    M.exp = @exponential;

    % Notice that the hash of two equivalent points will be different...
    M.hash = @(Y) ['z' hashmd5(Y(:))];

    M.rand = @() random(n, k);

    M.randvec = @randomvec;

    M.lincomb = @matrixlincomb;

    M.zerovec = @(Y) zeros(n, k);

    M.transp = @(Y1, Y2, d) projection(Y2, d);

    M.vec = @(Y, u_mat) u_mat(:);
    M.mat = @(Y, u_vec) reshape(u_vec, [n, k]);
    M.vecmatareisometries = @() true;

end

% Given a matrix X, returns the same matrix but with each column scaled so
% that they have unit 2-norm.
% See obliquefactory.
function X = normalize_rows(X)
    X = X';
	norms = sqrt(sum(X.^2, 1));
	X = bsxfun(@times, X, 1./norms);
    X = X';
end

% Orthogonal projection of each row of H to the tangent space at the
% corresponding row of X, seen as a point on a sphere.
% See obliquefactory.
function PXH = project_rows(X, H)
    X = X';
    H = H';
    % Compute the inner product between each vector H(:, i) with its root
    % point X(:, i), that is, X(:, i).' * H(:, i). Returns a row vector.
    inners = sum(X.*H, 1);
    % Subtract from H the components of the H(:, i)'s that are parallel to
    % the root points X(:, i).
    PXH = H - bsxfun(@times, X, inners);
    PXH = PXH';
end

% Projection onto the tangent space, i.e., on the tangent space of
% ||Y(i, :)|| = 1
function etaproj = projection(Y, eta)
    [unused, k] = size(Y); %#ok<ASGLU>
    eta = project_rows(Y, eta);

    % Projection onto the horizontal space
    YtY = Y'*Y;
    SS = YtY;
    AS = Y'*eta - eta'*Y;
    try
        % This is supposed to work and indeed return a skew-symmetric
        % solution Omega.
        Omega = lyap(SS, -AS);
    catch up %#ok<NASGU>
        % It can happen though that SS will be rank deficient. The
        % Lyapunov equation we solve still has a unique skew-symmetric
        % solution, but solutions with a symmetric part now also exist,
        % and the lyap function doesn't like that. So we want to
        % extract the minimum norm solution. This is also useful if lyap is
		% not available (it is part of the control system toolbox).
        mat = @(x) reshape(x, [k k]);
        vec = @(X) X(:);
        is_octave = exist('OCTAVE_VERSION', 'builtin');
        if ~is_octave
            [vecomega, unused] = minres(@(x) vec(SS*mat(x) + mat(x)*SS), vec(AS)); %#ok<NASGU>
        else
            [vecomega, unused] = gmres(@(x) vec(SS*mat(x) + mat(x)*SS), vec(AS)); %#ok<NASGU>
        end
        Omega = mat(vecomega);
    end
    % % Make sure the result is skew-symmetric (does not seem necessary).
    % Omega = (Omega-Omega')/2;
    etaproj = eta - Y*Omega;
end

% Retraction
function Ynew = retraction(Y, eta, t)
    if nargin < 3
        t = 1.0;
    end
    Ynew = Y + t*eta;
    Ynew = normalize_rows(Ynew);
end

% Exponential map
function Ynew = exponential(Y, eta, t)
    if nargin < 3
        t = 1.0;
    end

    Ynew = retraction(Y, eta, t);
    warning('manopt:elliptopefactory:exp', ...
        ['Exponential for fixed rank spectrahedron ' ...
        'manifold not implemented yet. Used retraction instead.\n' ...
        'To disable this warning: warning(''off'', ''manopt:elliptopefactory:exp'')']);
end

% Euclidean gradient to Riemannian gradient conversion.
% We only need the ambient space projection: the remainder of the
% projection function is not necessary because the Euclidean gradient must
% already be orthogonal to the vertical space.
function rgrad = egrad2rgrad(Y, egrad)
    rgrad = project_rows(Y, egrad);
end

% Euclidean Hessian to Riemannian Hessian conversion.
% TODO: speed this function up using bsxfun.
function Hess = ehess2rhess(Y, egrad, ehess, eta)
    k = size(Y, 2);

    % Directional derivative of the Riemannian gradient
    scaling_grad = sum((egrad.*Y), 2); % column vector of size n
    scaling_grad_repeat = scaling_grad*ones(1, k);

    Hess = ehess - scaling_grad_repeat.*eta;

    scaling_hess = sum((eta.*egrad) + (Y.*ehess), 2);
    scaling_hess_repeat = scaling_hess*ones(1, k);
    % directional derivative of scaling_grad_repeat
    Hess = Hess - scaling_hess_repeat.*Y;

    % Project on the horizontal space
    Hess = projection(Y, Hess);
end

% Random point generation on the manifold
function Y = random(n, k)
    Y = randn(n, k);
    Y = normalize_rows(Y);
end

% Random vector generation at Y
function eta = randomvec(Y)
    eta = randn(size(Y));
    eta = projection(Y, eta);
    nrm = norm(eta, 'fro');
    eta = eta / nrm;
end

author: MIT-SPARK

repo: GlobalOptimizationTutorial