text:
function save_model_path = rfcn_train(conf, imdb_train, roidb_train, varargin)
% save_model_path = rfcn_train(conf, imdb_train, roidb_train, varargin)
% --------------------------------------------------------
% R-FCN implementation
% Modified from MATLAB Faster R-CNN (https://github.com/shaoqingren/faster_rcnn)
% Copyright (c) 2016, Jifeng Dai
% Licensed under The MIT License [see LICENSE for details]
% --------------------------------------------------------

%% inputs
    ip = inputParser;
    ip.addRequired('conf',                              @isstruct);
    ip.addRequired('imdb_train',                        @iscell);
    ip.addRequired('roidb_train',                       @iscell);
    ip.addParamValue('do_val',          false,          @isscalar);
    ip.addParamValue('imdb_val',        struct(),       @isstruct);
    ip.addParamValue('roidb_val',       struct(),       @isstruct);
    ip.addParamValue('val_iters',       500,            @isscalar);
    ip.addParamValue('val_interval',    5000,           @isscalar);
    ip.addParamValue('snapshot_interval',...
                                        10000,          @isscalar);
    ip.addParamValue('solver_def_file', fullfile(pwd, 'models', 'rfcn_prototxts', 'ResNet-50L_res3a', 'solver_80k120k_lr1_3.prototxt'), ...
                                                        @isstr);
    ip.addParamValue('net_file',        fullfile(pwd, 'models', 'pre_trained_models', 'ResNet-50L', 'ResNet-50-model.caffemodel'), ...
                                                        @isstr);
    ip.addParamValue('cache_name',      'ResNet-50L_res3a', ...
                                                        @isstr);
    ip.addParamValue('caffe_version',   'Unkonwn',      @isstr);

    ip.parse(conf, imdb_train, roidb_train, varargin{:});
    opts = ip.Results;

%% try to find trained model
    imdbs_name = cell2mat(cellfun(@(x) x.name, imdb_train, 'UniformOutput', false));
    cache_dir = fullfile(pwd, 'output', 'rfcn_cachedir', opts.cache_name, imdbs_name);
    save_model_path = fullfile(cache_dir, 'final');
    if exist(save_model_path, 'file')
        return;
    end

%% init
    % set random seed
    prev_rng = seed_rand(conf.rng_seed);
    caffe.set_random_seed(conf.rng_seed);

    % init caffe solver
    mkdir_if_missing(cache_dir);
    caffe_log_file_base = fullfile(cache_dir, 'caffe_log');
    caffe.init_log(caffe_log_file_base);
    caffe_solver = caffe.Solver(opts.solver_def_file);
    caffe_solver.net.copy_from(opts.net_file);

    % init log
    timestamp = datestr(datevec(now()), 'yyyymmdd_HHMMSS');
    mkdir_if_missing(fullfile(cache_dir, 'log'));
    log_file = fullfile(cache_dir, 'log', ['train_', timestamp, '.txt']);
    diary(log_file);

    % set gpu/cpu
    if conf.use_gpu
        caffe.set_mode_gpu();
    else
        caffe.set_mode_cpu();
    end

    disp('conf:');
    disp(conf);
    disp('opts:');
    disp(opts);

%% making tran/val data
    fprintf('Preparing training data...');
    [image_roidb_train, bbox_means, bbox_stds] = rfcn_prepare_image_roidb(conf, opts.imdb_train, opts.roidb_train);
    fprintf('Done.\n');

    if opts.do_val
        fprintf('Preparing validation data...');
        [image_roidb_val] = rfcn_prepare_image_roidb(conf, opts.imdb_val, opts.roidb_val, bbox_means, bbox_stds);
        fprintf('Done.\n');

        % fix validation data
        shuffled_inds_val = generate_random_minibatch([], image_roidb_val, conf.ims_per_batch);
        shuffled_inds_val = shuffled_inds_val(randperm(length(shuffled_inds_val), opts.val_iters));
    end

%% training
    shuffled_inds = [];
    train_results = [];
    val_results = [];
    iter_ = caffe_solver.iter();
    max_iter = caffe_solver.max_iter();

    p = new_parpool(1);
    parfor i=1:1
        seed_rand(conf.rng_seed);
    end
    [shuffled_inds, sub_db_inds] = generate_random_minibatch(shuffled_inds, image_roidb_train, conf.ims_per_batch);
    parHandle = parfeval(p, @rfcn_get_minibatch, 1, conf, image_roidb_train(sub_db_inds));
    tic
    while (iter_ < max_iter)
        caffe_solver.net.set_phase('train');

        % generate minibatch training data
        % gather date
        [~, net_inputs] = fetchNext(parHandle);

        % generate minibatch training data
        % generate data asynchronously
        [shuffled_inds, sub_db_inds] = generate_random_minibatch(shuffled_inds, image_roidb_train, conf.ims_per_batch);
        parHandle = parfeval(p, @rfcn_get_minibatch, 1, conf, image_roidb_train(sub_db_inds));

        caffe_solver.net.reshape_as_input(net_inputs);

        % one iter SGD update
        caffe_solver.net.set_input_data(net_inputs);
        caffe_solver.step(1);

        rst = caffe_solver.net.get_output();
        train_results = parse_rst(train_results, rst);

        % do valdiation per val_interval iterations
        if ~mod(iter_, opts.val_interval)
            if opts.do_val
                caffe_solver.net.set_phase('test');
                for i = 1:length(shuffled_inds_val)
                    sub_db_inds = shuffled_inds_val{i};
                    net_inputs = rfcn_get_minibatch(conf, image_roidb_val(sub_db_inds));
                    caffe_solver.net.reshape_as_input(net_inputs);

                    caffe_solver.net.forward(net_inputs);

                    rst = caffe_solver.net.get_output();
                    val_results = parse_rst(val_results, rst);
                end
            end

            show_state(iter_, train_results, val_results);
            toc;tic;
            train_results = [];
            val_results = [];
            diary; diary; % flush diary
        end

        % snapshot
        if ~mod(iter_, opts.snapshot_interval)
            snapshot(caffe_solver, bbox_means, bbox_stds, cache_dir, sprintf('iter_%d', iter_));
        end

        iter_ = caffe_solver.iter();
    end

    % final snapshot
    snapshot(caffe_solver, bbox_means, bbox_stds, cache_dir, sprintf('iter_%d', iter_));
    save_model_path = snapshot(caffe_solver, bbox_means, bbox_stds, cache_dir, 'final');

    diary off;
    caffe.reset_all();
    rng(prev_rng);
end

function [shuffled_inds, sub_inds] = generate_random_minibatch(shuffled_inds, image_roidb_train, ims_per_batch)

    % shuffle training data per batch
    if isempty(shuffled_inds)
        % make sure each minibatch, only has horizontal images or vertical
        % images, to save gpu memory

        hori_image_inds = arrayfun(@(x) x.im_size(2) >= x.im_size(1), image_roidb_train, 'UniformOutput', true);
        vert_image_inds = ~hori_image_inds;
        hori_image_inds = find(hori_image_inds);
        vert_image_inds = find(vert_image_inds);

        % random perm
        lim = floor(length(hori_image_inds) / ims_per_batch) * ims_per_batch;
        hori_image_inds = hori_image_inds(randperm(length(hori_image_inds), lim));
        lim = floor(length(vert_image_inds) / ims_per_batch) * ims_per_batch;
        vert_image_inds = vert_image_inds(randperm(length(vert_image_inds), lim));

        % combine sample for each ims_per_batch
        hori_image_inds = reshape(hori_image_inds, ims_per_batch, []);
        vert_image_inds = reshape(vert_image_inds, ims_per_batch, []);

        shuffled_inds = [hori_image_inds, vert_image_inds];
        shuffled_inds = shuffled_inds(:, randperm(size(shuffled_inds, 2)));

        shuffled_inds = num2cell(shuffled_inds, 1);
    end

    if nargout > 1
        % generate minibatch training data
        sub_inds = shuffled_inds{1};
        assert(length(sub_inds) == ims_per_batch);
        shuffled_inds(1) = [];
    end
end

function model_path = snapshot(caffe_solver, bbox_means, bbox_stds, cache_dir, file_name)
    bbox_pred_layer_name = 'rfcn_bbox';
    weights = caffe_solver.net.params(bbox_pred_layer_name, 1).get_data();
    biase = caffe_solver.net.params(bbox_pred_layer_name, 2).get_data();
    weights_back = weights;
    biase_back = biase;

    rep_time = size(weights, 4)/length(bbox_means(:));

    bbox_stds_flatten = bbox_stds';
    bbox_stds_flatten = bbox_stds_flatten(:);
    bbox_stds_flatten = repmat(bbox_stds_flatten, [1,rep_time])';
    bbox_stds_flatten = bbox_stds_flatten(:);
    bbox_stds_flatten = permute(bbox_stds_flatten, [4,3,2,1]);

    bbox_means_flatten = bbox_means';
    bbox_means_flatten = bbox_means_flatten(:);
    bbox_means_flatten = repmat(bbox_means_flatten, [1,rep_time])';
    bbox_means_flatten = bbox_means_flatten(:);
    bbox_means_flatten = permute(bbox_means_flatten, [4,3,2,1]);

    % merge bbox_means, bbox_stds into the model
    weights = bsxfun(@times, weights, bbox_stds_flatten); % weights = weights * stds;
    biase = biase .* bbox_stds_flatten(:) + bbox_means_flatten(:); % bias = bias * stds + means;

    caffe_solver.net.set_params_data(bbox_pred_layer_name, 1, weights);
    caffe_solver.net.set_params_data(bbox_pred_layer_name, 2, biase);

    model_path = fullfile(cache_dir, file_name);
    caffe_solver.net.save(model_path);
    fprintf('Saved as %s\n', model_path);

    % restore net to original state
    caffe_solver.net.set_params_data(bbox_pred_layer_name, 1, weights_back);
    caffe_solver.net.set_params_data(bbox_pred_layer_name, 2, biase_back);
end

function show_state(iter, train_results, val_results)
    fprintf('\n------------------------- Iteration %d -------------------------\n', iter);
    fprintf('Training : accuracy %.3g, loss (cls %.3g, reg %.3g)\n', ...
        mean(train_results.accuarcy.data), ...
        mean(train_results.loss_cls.data), ...
        mean(train_results.loss_bbox.data));
    if exist('val_results', 'var') && ~isempty(val_results)
        fprintf('Testing  : accuracy %.3g, loss (cls %.3g, reg %.3g)\n', ...
            mean(val_results.accuarcy.data), ...
            mean(val_results.loss_cls.data), ...
            mean(val_results.loss_bbox.data));
    end
end

repo: R-FCN


text:
function [subband size_band] = complex_divisive_normalized(pyro,pind,Nsc,Nor,parent,neighbor,blSzX,blSzY)

Nband = size(pind,1)-1;

p = 1;
for scale=1:Nsc-1
    for orien=1:Nor
        nband = (scale-1)*Nor+orien+1; % except the ll
        %      if(prod(band-nband-1) ~=0)
        %            continue;
        %      end
        aux_c = pyrBand(pyro, pind, nband);
        aux = abs(aux_c);
        [Nsy,Nsx] = size(aux);

        prnt = parent & (nband <= Nband-(Nsc-1)*Nor);  % has the subband a parent?
%       define that only the finest scale has a parent
        BL = zeros(size(aux,1),size(aux,2),1 + prnt);
        BL(:,:,1) = aux;
        if prnt,
            auxp = pyrBand(pyro, pind, nband+Nor);
            %    if nband>Nor+1,     % resample 2x2 the parent if not in the high-pass oriented subbands.
            % 	   auxp = real(imenlarge2(auxp)); % this was uncommented
            auxp = abs(imresize(auxp,2)); %
            %   end
            %  fprintf('parent band and size is %d %d %d \n',nband+Nor,Nsy,Nsx)
            BL(:,:,2) = auxp(1:Nsy,1:Nsx);
        end
        y = BL;
        [nv,nh,nb] = size(y);
        block = [blSzX blSzY];

        nblv = nv-block(1)+1;	% Discard the outer coefficients
        nblh = nh-block(2)+1;   % for the reference (centrral) coefficients (to avoid boundary effects)
        nexp = nblv*nblh;			% number of coefficients considered
        N = prod(block) + prnt; % size of the neighborhood

        Ly = (block(1)-1)/2;		% block(1) and block(2) must be odd!
        Lx = (block(2)-1)/2;
        if (Ly~=floor(Ly))|(Lx~=floor(Lx)),
            error('Spatial dimensions of neighborhood must be odd!');
        end
        Y = zeros(nexp,N);		% It will be the observed signal (rearranged in nexp neighborhoods)
        % Rearrange observed samples in 'nexp' neighborhoods
        n = 0;
        for ny=-Ly:Ly,	% spatial neighbors
            for nx=-Lx:Lx,
                n = n + 1;
                foo = shift(y(:,:,1),[ny nx]);
                foo = foo(Ly+1:Ly+nblv,Lx+1:Lx+nblh);
                Y(:,n) = (foo(:));
            end
        end

        if prnt,	% parent
            n = n + 1;
            foo = y(:,:,2);
            foo = foo(Ly+1:Ly+nblv,Lx+1:Lx+nblh);
            Y(:,n) = (foo(:));
        end

        %      including neighbor
        if neighbor,
            for neib=1:Nor
                if neib == orien
                    continue;
                end
                n=n+1;
                nband1 = (scale-1)*Nor+neib+1; % except the ll
                aux1 = abs(pyrBand(pyro, pind, nband1));
                aux1 = aux1(Ly+1:Ly+nblv,Lx+1:Lx+nblh);
                Y(:,n) = (aux1(:));
            end
        end

        C_x = innerProd(Y)/nexp;
        % C_x is positive definete covariance matrix
        [Q,L] = eig(C_x);
        % correct possible negative eigenvalues, without changing the overall variance
        L = diag(diag(L).*(diag(L)>0))*sum(diag(L))/(sum(diag(L).*(diag(L)>0))+(sum(diag(L).*(diag(L)>0))==0));
        C_x = Q*L*Q';

        o_c = aux_c(Ly+1:Ly+nblv,Lx+1:Lx+nblh);
        o_c = (o_c(:));
%         o_c_r = real(o_c);
%         o_c_i = imag(o_c);
%         o_c_r = o_c_r - mean(o_c_r);
%         o_c_i = o_c_i - mean(o_c_i);
        o_c = o_c - mean(o_c);

        tempY = (Y*inv(C_x)).*Y/N;
        z = sqrt(sum(tempY,2));
        ind = find(z~=0);

        g_c = o_c(ind)./z(ind);
        size_band(p,1) = nblv;
        size_band(p,2) = nblh;
        g_c = g_c - mean(g_c);
        g_c_m = reshape(g_c, nblv, nblh);
        subband{p} = g_c_m;
        p = p+1;

    end
end
return;

author: vztu

repo: VIDEVAL

save_path: github-repos/MATLAB/vztu-VIDEVAL


text:
function [A, B, C, D, tau, kappa, epsilon] = theta2abcd(theta, idx, R, L )
% Transform DCM parameters from vectorized to structured format
%
% This is a protected method of the tapas_Huge class. It cannot be called
% from outside the class.
%

% Author: Yu Yao (yao@biomed.ee.ethz.ch)
% Copyright (C) 2019 Translational Neuromodeling Unit
%                    Institute for Biomedical Engineering,
%                    University of Zurich and ETH Zurich.
%
% This file is part of TAPAS, which is released under the terms of the GNU
% General Public Licence (GPL), version 3. For further details, see
% <https://www.gnu.org/licenses/>.
%
% This software is provided "as is", without warranty of any kind, express
% or implied, including, but not limited to the warranties of
% merchantability, fitness for a particular purpose and non-infringement.
%
% This software is intended for research only. Do not use for clinical
% purpose. Please note that this toolbox is under active development.
% Considerable changes may occur in future releases. For support please
% refer to:
% https://github.com/translationalneuromodeling/tapas/issues
%

tmp = zeros(idx.P_f, 1);
tmp([idx.clustering; idx.homogenous]) = theta;
% hemodynamic parameters
iEd = numel(tmp);
epsilon = tmp(iEd); % ratio of intra- and extravascular signal
iEd = iEd - 1;
kappa = tmp(iEd-R+1:iEd); % decay of vasodilatory signal
iEd = iEd - R;
tau = tmp(iEd-R+1:iEd); % transit time
iEd = iEd - R;

% connectivity parameters
iSt = 0;
A = reshape(tmp(iSt+1:iSt+R^2), R, R);
iSt = iSt + R^2;
B = reshape(tmp(iSt+1:iSt+R^2*L), R, R, L);
iSt = iSt + R^2*L;
C = reshape(tmp(iSt+1:iSt+R*L), R, L);
iSt = iSt + R*L;
% nonlinear connections
if iEd - iSt == R^3
    D = reshape(tmp(iSt+1:iEd), R, R, R);
else
    D = zeros(R, R, 0);
end

end

author: translationalneuromodeling

repo: tapas


text:
    % design_matrix: data class for creating a design matrix to be used with a linear model including fmri data.
    %
    %
    % Inputs:
    % ---------------------------------------------------------------------
    % dat                       : M x N numeric matrix containing Observations and Variables
    %                           If dat is a file name then will try to
    %                           import data into design_matrix object.
    %                           Make sure there is only text in header and
    %                           not anywhere else.
    %
    % varname                   : Cell array containing variable names.  Must
    %                             match number of column in data matrix
    %
    % Examples:
    % ---------------------------------------------------------------------
    % DM = design_matrix([ones(10,1), (1:10)', (1:10).^2'],{'Intercept','X','X2'})
    %
    % Original version: Copyright Luke Chang 2/2014

    % Notes:
    % Need to add these Methods:
    % -create regressor from stim times
    % -pca

classdef design_matrix < handle
    properties
        dat = [];
        varname = {};
        fname = '';
    end

    methods
        function obj = design_matrix(dat, varname)
            class constructor % Initialize instance of design_matrix

            if(nargin > 1)
                try
                    if(~ismatrix(dat) || ~isnumeric(dat) || iscell(dat))
                        error('Make sure input data is a matrix')
                    end
                    if(length(varname) ~= size(dat,2) || ~iscell(varname));
                        error('Make sure the number of variable names corresponds to number of data columns.')
                    end
                    obj.dat = dat;
                    obj.varname = varname;
                catch err
                    error('Make sure input variable names are in a cell array with length equal to number of data columns and data is a matrix.')
                end
                obj.varname = varname;

            elseif(nargin > 0)

                %import data into design_matrix if dat is a filename
                if ischar(dat)
                    ftest = exist(dat,'file'); %check if valid file
                    if ftest == 2
                        data = importdata(dat);
                        obj.dat = data.data;
                        obj.varname = data.colheaders;
                    end

                elseif(~ismatrix(dat) || ~isnumeric(dat) || iscell(dat))
                    error('Make sure input data is a matrix')
                else

                    obj.dat = dat;
                end

            else % if nothing initialize empty object
                return
            end
        end

        function names(obj)
            % names(obj)
            %
            % List variable names for each regressor
            display(obj.varname)
        end

        function dim = size(obj, varargin)
            % dim = size(obj, varargin)
            %
            % Return dimensions of design matrix
            % Optional Input: Indicate Dimension(row = 1 or column = 2)

            if nargin > 1
                dim = size(obj.dat, varargin{1});
            else
                dim = size(obj.dat);
            end
        end

        function plot(obj)
            % plot(obj)
            %
            % Plot design matrix
            imagesc(obj.dat)
        end

        function save(obj, fname)
            % save(obj, fname)
            %
            % Save Design Matrix to file

            save(fname, obj)
        end

        function obj = addvariable(obj, x, varargin)
            % obj = addvariable(obj, x, varargin)
            %
            % Add regressor to design matrix
            %
            % optional inputs
            % -------------------------------------------------------------------
            % 'Name'        : 'Name' followed by Variable name
            %                  Default is 'newVx'
            %
            % 'Order'       : 'Order' followed by where to insert new data columns location
            %                  Default is end

            % Check Inputs
            if size(x,1) ~= size(obj,1)
                error('Make sure new variable column is the same length as design matrix')
            end

            % Defaults
            for i = 1:size(x,2)
                newname{i} = ['newV' num2str(i)];
            end
            isnewname = 0;
            varorder = size(obj,2); %add to end

            % Parse inputs
            % -------------------------------------------------------------------
            for varg = 1:length(varargin)
                if ischar(varargin{varg})
                    % reserved keywords
                    if strcmpi('name',varargin{varg})
                        if length(varargin{varg + 1}) == size(x,2)
                            newname = varargin{varg + 1};
                            isnewname = 1;
                        else
                            error('Make sure ''Name'' is follwed by a valid variable name')
                        end
                        varargin{varg} = {}; varargin{varg + 1} = {};
                    end

                    if strcmpi('order',varargin{varg})
                        if isnumeric(varargin{varg + 1}) && varargin{varg + 1} <= size(obj,2)
                            varorder = varargin{varg + 1};
                        else
                            error('Make sure ''Order'' is follwed by a valid column number')
                        end
                        varargin{varg} = {}; varargin{varg + 1} = {};
                    end
                end
            end
            % -------------------------------------------------------------------

            % Update Name vector if not empty
            if ~isempty(obj.varname)
                if varorder == 1 %begin
                    obj.varname = [newname, obj.varname];
                elseif varorder == size(obj,2) %end
                    obj.varname = [obj.varname, newname];
                else %Somewhere inbetween
                    obj.varname = [obj.varname(1:varorder), newname, obj.varname(varorder + 1 : end)];
                end
            end

            % Add new variables to design Matrix
            if varorder == 1 %begin
                obj.dat = [x, obj.dat];
            elseif varorder == size(obj,2) %end
                obj.dat = [obj.dat, x];
            else %Somewhere inbetween
                obj.dat = [obj.dat(:,1:varorder), x, obj.dat(:,varorder + 1 : end)];
            end
        end

        function obj = zscore(obj, varargin)
            % obj = zscore(obj, varargin)
            %
            % Standardize columns of design matrix
            %
            % optional inputs
            % -------------------------------------------------------------------
            % 'center'        : Only remove mean, don't standardize

            % Defaults
            center = 0; %Only remove mean

            % optional inputs
            % -------------------------------------------------------------------
            for varg = 1:length(varargin)
                if ischar(varargin{varg})
                    % reserved keywords
                    if strcmpi('center',varargin{varg})
                        center=1;
                        varargin{varg} = {};
                    end
                end
            end
            % -------------------------------------------------------------------

            if center
                obj.dat = obj.dat - repmat(mean(obj.dat),size(obj,1),1);
            else
                obj.dat = zscore(obj.dat);
            end
        end

        function obj = removevariable(obj, x)
            % obj = removevariable(obj, x)
            %
            % Remove columns from design_matrix
            %
            % Inputs
            % -------------------------------------------------------------------
            % x             : Input vector of columns to remove

            obj.dat(:,x) = [];
            obj.varname(x) = [];
        end

        function obj = addintercept(obj)
            % obj = addintercept(obj)
            %
            % Add intercept to design matrix

            %Check if Intercept exists
            if sum(strcmpi(obj.varname,'intercept')) > 0
                error('Intercept Name already included in obj.varname')
            end

            %Check if any variable only includes ones
            for i = 1:size(obj,2)
                if sum(obj.dat(:,i)==1) == size(obj,1)
                    error('There is already a column of ones that resembles an intercept')
                end
            end

            obj.dat = [obj.dat, ones(size(obj,1),1)];
            obj.varname = [obj.varname, 'Intercept'];
        end

        function obj = removeintercept(obj)
            %  obj = removeintercept(obj)
            %
            % Remove intercept from design matrix

            %Check if any variable only includes ones or if intercept is in varname
            for i = 1:size(obj,2)
                whereint(i) = sum(obj.dat(:,i)==1) == size(obj,1);
            end
            if sum(whereint) == 0
                error('There does not appear to be any column of ones that resembles an intercept')
            elseif sum(strcmpi('intercept',obj.varname)) == 0
                error('Intercept Name is not included in obj.varname')
            end

            % now remove intercept
            obj.dat(:,whereint) = [];
            obj.varname(whereint) = [];
        end

        function vif = vif(obj, varargin)
            % vif = vif(obj, varargin)
            %
            % Check for multicollinearity by getting variance inflation factors
            %
            % See original getvif.m in canlab repository - OptimizeDesign11/core_functions/getvif.m
            %
            % optional inputs
            % -------------------------------------------------------------------
            % 'nointercept'       : Remove intercept (turned off by default)

            % Defaults
            noint = 0;
            if strcmpi('nointercept',varargin)
                noint = 1;
            end

            %Remove intercept if asked
            if noint
                obj = removeintercept(obj);
            end

            %Calculate VIF
            for i = 1:size(obj.dat,2)
                X = obj.dat;
                y = X(:,i);
                X(:,i) = [];
                b = X\y;
                fits = X * b;
                rsquare = var(fits) / var(y);

                if rsquare == 1,rsquare = .9999999;end

                vif(i) = 1 / (1 - rsquare);
            end
        end

        function r = corr(obj)
            % Calculate pairwise correlation of regressors in design_matrix

            r = corr(obj.dat);
        end

        function obj = normalizedrank(obj, varargin)
            % obj = normalizedrank(obj, varargin)
            %
            % Rank each regressor and normalize between [0,1]
            %
            % See normalizedrank.m for optional inputs

            obj.dat = normalizedrank(obj.dat, varargin);
        end

        function obj = conv_hrf(obj, varargin)
            % obj = conv_hrf(obj, varargin)
            %
            % Convolve each regressors with hemodynamic response function
            % Uses spm_hrf.m
            %
            % optional inputs
            % -------------------------------------------------------------------
            % 'tr'           : Input TR to use for creating HRF
            %                  (e.g., 'tr', 3)
            %
            % 'select'       : Select Input vector of regressors to convolve
            %                  (e.g., 'select', [2,4])
            %
            % 'custom_hrf'   : Use custom HRF (e.g., 'custom_hrf', [1.00, 0.59, 0.39, 0.27]

            % Defaults
            include = (1:size(obj,2)); %Convolve entire Design Matrix by default
            tr = 2;

            % Check if spm_hrf is on path
            checkspm =  which('spm_hrf.m');
            if isempty(checkspm), error('Make sure spm is in matlab path'); end

            % Parse inputs
            % -------------------------------------------------------------------
            for varg = 1:length(varargin)
                if ischar(varargin{varg})
                    % reserved keywords
                    if strcmpi('tr',varargin{varg})
                        if isnumeric(varargin{varg + 1})
                            tr = varargin{varg + 1};
                            crf = spm_hrf(tr);
                        else
                            error('Make sure ''tr'' is followed by valid number')
                        end
                        varargin{varg} = {}; varargin{varg + 1} = {};
                    end

                    if strcmpi('select',varargin{varg})
                        if isnumeric(varargin{varg + 1}) && all(varargin{varg + 1} <= size(obj,2))
                            include = varargin{varg + 1};
                        else
                            error('Make sure ''select'' is followed by a valid column number')
                        end
                        varargin{varg} = {}; varargin{varg + 1} = {};
                    end

                    if strcmpi('custom_hrf',varargin{varg})
                        if isnumeric(varargin{varg + 1}) && all(varargin{varg + 1} <= size(obj,2))
                            crf = varargin{varg + 1};
                        else
                            error('Make sure ''custom_hrf'' is followed by a vector of a canonical response function')
                        end
                        varargin{varg} = {}; varargin{varg + 1} = {};
                    end
                end
            end
            % -------------------------------------------------------------------

            %Convolution of task
            for i = 1:length(include)
                convdat = conv(obj.dat(:,include(i)),crf);

                %Cut off extra data from convolution
                obj.dat(:,include(i)) = convdat(1:size(obj,1));
            end
        end

        function obj = hpfilter(obj, varargin)
            % obj = hpfilter(obj, varargin)
            %
            % Add High pass filter design matrix using spm's discrete
            % cosine Transform (spm_filter.m)
            %
            % optional inputs
            % -------------------------------------------------------------------
            % 'tr'           : Input TR to use for creating HRF
            %                  (e.g., 'tr', 3) Default = 2;
            %
            % 'duration'     : Duration of high pass filter in seconds
            %                  (e.g., 'duration', 100) Default = 180;

            % Defaults
            tr = 2;
            filterlength = 180;

            % Check if spm_hrf is on path
            checkspm =  which('spm_filter.m');
            if isempty(checkspm), error('Make sure spm is in matlab path'); end

            % Parse inputs
            % -------------------------------------------------------------------
            for varg = 1:length(varargin)
                if ischar(varargin{varg})
                    % reserved keywords
                    if strcmpi('tr',varargin{varg})
                        if isnumeric(varargin{varg + 1})
                            tr = varargin{varg + 1};
                        else
                            error('Make sure ''tr'' is followed by valid number')
                        end
                        varargin{varg} = {}; varargin{varg + 1} = {};
                    end
                    if strcmpi('duration',varargin{varg})
                        if isnumeric(varargin{varg + 1})
                            filterlength = varargin{varg + 1};
                        else
                            error('Make sure ''duration'' is followed by valid number')
                        end
                        varargin{varg} = {}; varargin{varg + 1} = {};
                    end
                end
            end
            % -------------------------------------------------------------------

            %create high pass filter
            K.RT = tr;
            K.row = 1:size(obj,1);
            K.HParam = filterlength;
            nK = spm_filter(K);
            if isempty(nK.X0), error('Check if filter duration is too long'); end

            %Add filter to design_matrix
            obj.dat = [obj.dat, nK.X0];

            %Add variable names
            for i = 1:size(nK.X0,2)
                filtname{i} = ['hpfilter' num2str(i)];
            end
            obj.varname = [obj.varname, filtname];
        end

        function stats = regress(obj, Y, varargin)
            % [B,BINT,R,RINT,STATS] = regress(obj, Y)
            %
            % Regress design matrix on vector Y
            % Uses matlab's regress function
            %
            % optional inputs
            % -------------------------------------------------------------------
            % 'robust'           : use robust regression

            % Defaults
            doRobust = 0;
            for i = 1:length(varargin)
                if strcmpi(varargin(i),'robust')
                    doRobust = 1;
                    find_robust = exist('robustfit');
                    if find_robust ~= 2
                        error('Make sure robustfit is on your path, requires the stats toolbox')
                    end
                    varargin(i) = [];
                end
            end

            if ~doRobust
                [stats.B, stats.BINT, stats.R, stats.RINT, stats.STATS] = regress(Y, obj.dat);
            else
                [stats.B, stats.STATS] = robustfit(obj.dat, Y,[],[],'off');
            end
        end

        function obj = onsettimes(obj, onset, names, tr, timing )
            % obj = onsettimes(obj, onset, names, tr, timing )
            %
            %Create stimulus regressor from onset times
            %
            % Inputs
            % -------------------------------------------------------------------
            % onset        : Input cell array of onset times for each
            %                   regressor in FSL's 3 column format (e.g., onset in sec, duration in sec, weight).
            %
            % names        : Cell array of variable names corresponding
            %                to each onset cell (e.g., {'BlueOn','RedOn'})
            %
            % tr           : Repetition time (e.g., 2)
            %
            % timing       : Timing converstion from onset array to design matrix
            %                  (e.g., 'sec2tr','tr2sec','sec2sec',or
            %                  'tr2tr'). Need to know which format each
            %                  array is in.

            %Convert Onset Times Into Boxcar Regressors
            r = zeros(size(obj,1),length(onset));
            for i = 1:length(onset)
                for j = 1:size(onset{i},1)
                    switch timing
                        case 'sec2tr'
                            if floor(onset{i}(j,1)/tr) == 0
                                r(1 : 1 + ceil(onset{i}(j,2)), i) = onset{i}(j,3);
                            else
                                r(floor(onset{i}(j,1) / tr) : floor(onset{i}(j,1) / tr) + ceil(onset{i}(j,2) / tr) - 1, i) = onset{i}(j,3);
                            end
                        case 'tr2sec'
                            r(floor(onset{i}(j,1) * tr) : floor(onset{i}(j,1) * tr) + ceil(onset{i}(j,2) * tr) - 1, i) = onset{i}(j,3);
                        case {'sec2sec', 'tr2tr'}
                            r(floor(onset{i}(j,1)) : floor(onset{i}(j,1)) + ceil(onset{i}(j,2)) - 1, i) = onset{i}(j,3);
                    end
                end
            end
            obj.dat = [obj.dat, r];

            %Add Variable names
            obj.varname = [obj.varname, names];
        end

        function obj = write(obj, varargin)
            % obj = write(obj, varargin)
            %
            % -------------------------------------------------------------------
            % write design_matrix object into csv file.  Will use obj.fname
            % or can specify optional name
            %
            % -------------------------------------------------------------------
            % Optional Inputs
            % -------------------------------------------------------------------
            % fname        : path and file name of csv file.
            % -------------------------------------------------------------------

            if nargin > 1 %use supplied file name
                if ischar(varargin{1})
                    hdr = sprintf('%s,',obj.varname{:});
                    hdr(end) = '';
                    dlmwrite(varargin{1}, hdr,'') %Write Header 1st
                    dlmwrite(varargin{1}, obj.dat, 'delimiter',',','-append','precision',10) %Append data
                end

            elseif ~isempty(obj.fname) %use obj.fname
                hdr = sprintf('%s,',obj.varname{:});
                hdr(end) = '';
                dlmwrite(obj.fname, hdr, ''); %Write Header 1st
                dlmwrite(obj.fname, obj.dat, 'delimiter',',','-append','precision',10) %Append data
            else
                error('Please supply valid file name with path to save.')
            end

        end

        function c = horzcat(varargin)
            % function c = horzcat(varargin)
            % -------------------------------------------------------------------
            % Implements the horzcat ([a b]) operator on design_matrix objects across variables.
            % Requires that each object has an equal number of rows
            % -------------------------------------------------------------------
            % Examples:
            % c = [dm1 dm2];
            % -------------------------------------------------------------------

            %check if number of rows is the same
                        %check if varnames are the same
            nrow = [];
            for i = 1:nargin
                nrow(i) = size(varargin{i},1);
            end
            for i = 1:nargin
                for j = 1:nargin
                    if nrow(i)~=nrow(j)
                        error('objects have a different number of rows')
                    end
                end
            end

            dat = [];
            varname = [];
            for i = 1:nargin
                    %Check if design_matrix object
                if ~isa(varargin{i}, 'design_matrix')
                    error('Input Data is not an design_matrix object')
                end

                dat = [dat, varargin{i}.dat];
                varname = [varname, varargin{i}.varname];
            end

            c = varargin{1};
            c.dat = dat;
            c.varname = varname;
        end

        function c = vertcat(varargin)
            % function c = vertcat(varargin)
            % -------------------------------------------------------------------
            % Implements the vertcat ([a b]) operator on design_matrix objects across rows.
            % Requires that each object has an equal number of columns and
            % that each varname is identical
            % -------------------------------------------------------------------
            % Examples:
            % c = [dm1; dm2];
            % -------------------------------------------------------------------

            %check if varnames are the same
            varname = {};
            for i = 1:nargin
                varname{i} = varargin{i}.varname;
            end
            for i = 1:nargin
                for j = 1:nargin
                    if ~strcmpi(varname{i},varname{j})
                        error('variable names do not match')
                    end
                end
            end
            dat = [];
            for i = 1:nargin
                %Check if design_matrix object
                if ~isa(varargin{i}, 'design_matrix')
                    error('Input Data is not an design_matrix object')
                end
                dat = [dat; varargin{i}.dat];
            end

            c = varargin{1};
            c.dat = dat;
        end

    end %methods
end %class

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [maxFreqDev,omega,freqs] = calcFreqDev(stimList,conditions,freqConditions)
% function [maxFreqDev,omega,freqs] = calcFreqDev(stimList,conditions,freqConditions)
%
%
% conditions is integer vector of conditions of interest
% freqConditions is vector of propotions of all conditions

			nc = size(conditions,2);

            for i = 1:nc
         		freqs(i) = sum(stimList == conditions(i));
      	    end

			% not necessary if forcing sum to 1
			% freqs =  freqs ./ size(stimList,1);

			% force frequencies to sum to 1 b/c otherwise rest intervals will influence frequencies
            freqs = freqs ./ sum(freqs);

			maxFreqDev = max(abs(freqConditions(1:nc) - freqs));

			omega = 1 - maxFreqDev;

return


text:
function y = apply_multiple_ouput(f, g, x)

% apply_multiple_ouput - compose multiple arguments functions
%
%   y = apply_multiple_ouput(f, g, x);
%
% compute y=f(g(x)) where f takes several input and g computes several outputs.
%
%   Copyright (c) 2010 Gabriel Peyre

if nargin(f)>=0 & nargout(g)>=0 & not(nargout(g)==nargin(f))
    error('one must have nargout(g)==nargin(f)');
end

switch nargin(f)
    case 1
        [a{1}] = g(x);
        y = f(a{1});
    case 2
        [a{1}, a{2}] = g(x);
        y = f(a{1}, a{2});
    case 3
        [a{1}, a{2}, a{3}] = g(x);
        y = f(a{1}, a{2}, a{3});
    case 4
        [a{1}, a{2}, a{3}, a{4}] = g(x);
        y = f(a{1}, a{2}, a{3}, a{4});
    otherwise
        error('Not implemented');
end

author: gpeyre

repo: matlab-toolboxes


text:
function cirbva()
    %   This subroutine "circle"  selects the Ni closest earthquakes
    %   around a interactively selected point.  Resets ZG.newcat and ZG.newt2
    %   Operates on "primeCatalog".
    %   changes newt2
    %
    % axis: h1
    % plots to: plos1 as xk
    % inCatalog: a
    % outCatalog: newt2, newcat
    % mouse controlled
    % closest events OR radius
    % calls: bdiff
    %
    % turned into function by Celso G Reyes 2017

    ZG=ZmapGlobal.Data; % used by get_zmap_globals

    %  Input Ni:
    %
    report_this_filefun();
    ZG=ZmapGlobal.Data;

    delete(findobj('Tag','plos1'));
    disp(ZmapGlobal.Data.hold_state)

    axes(h1)
    %zoom off

    titStr ='Selecting EQ in Circles                         ';
    messtext= ...
        ['                                                '
        '  Please use the LEFT mouse button              '
        ' to select the center point.                    '
        ' The "ni" events nearest to this point          '
        ' will be selected and displayed in the map.     '];

    msg.dbdisp(messtext, titStr);

    % Input center of circle with mouse
    %
    [xa0,ya0]  = ginput(1);

    stri1 = [ 'Circle: ' num2str(xa0,5) '; ' num2str(ya0,4)];
    stri = stri1;
    pause(0.1)
    %  calculate distance for each earthquake from center point
    %  and sort by distance
    %
    l = ZG.primeCatalog.epicentralDistanceTo(ya0,xa0);
    ZG.newt2=ZG.primeCatalog; % points to same thing
    if met == 'ni'
        % take first ni and sort by time
        [ZG.newt2, R2] = ZG.primeCatalog.selectClosestEvents(ni);
    elseif  met == 'ra'
        ZG.newt2 = ZG.primeCatalog.selectRadius(ra);
        R2=ra;
    elseif met == 'ti'
        global t1 t2 t3 t4

        lt =  ZG.newt2.Date >= t1 &  ZG.newt2.Date <t2 ;
        bdiff(ZG.newt2.subset(lt));
        ZG.hold_state=true;
        lt =  ZG.newt2.Date >= t3 &  ZG.newt2.Date <t4 ;
        bdiff(ZG.newt2.subset(lt));

    end
    ZG.newt2.sort('Date');
    R2 = ra;

    %
    % plot Ni clostest events on map as 'x':

    set(gca,'NextPlot','add')
    plot(ZG.newt2.Longitude,ZG.newt2.Latitude,'xk','Tag','plos1');

    % plot circle containing events as circle
    x = -pi-0.1:0.1:pi;
    pl = plot(xa0+sin(x)*R2/(cosd(ya0)*111), ya0+cos(x)*R2/(cosd(ya0)*111),'k')
    %plot(xa0+sin(x)*l(ni)/111, ya0+cos(x)*l(ni)/111,'k')

    set(gcf,'Pointer','arrow')

    %
    newcat = ZG.newt2;                   % resets ZG.newcat and ZG.newt2

    bdiff(ZG.newt2,ZmapGlobal.Data.hold_state)

end

author: CelsoReyes


text:
%   AUTHORSHIP
%   Primary Developer: Stephen Meehan <swmeehan@stanford.edu>
%   Math Lead & Secondary Developer:  Connor Meehan <connor.gw.meehan@gmail.com>
%   Bioinformatics Lead:  Wayne Moore <wmoore@stanford.edu>
%   Provided by the Herzenberg Lab at Stanford University
%   License: BSD 3 clause
%
classdef ToolBar < handle
    properties
        ht=[];
        comboBoxes;
        bag=[];
        firstItems={};
        jToolbar=[];
    end
    methods
        function setEnabled(tb, yes)
            N=tb.jToolbar.getComponentCount;
            for i=0:N-1
                btn=tb.jToolbar.getComponentAtIndex(i);
                btn.setRequestFocusEnabled(false);
                if isa(btn, 'javax.swing.JPanel')
                    btn.setVisible(yes);
                else
                    btn.setEnabled(yes);
                end
            end
        end
        function this=ToolBar(ht)
            this.comboBoxes=[];
            this.bag=struct();
            this.jToolbar=ToolBarMethods.getJ(ht);
            this.ht=ht;
        end
    end

    methods(Static)

        function this=New(fig, first, removeZoom, removeRotate3D, removeEdit)
            if nargin<5
                removeEdit=true;
                if nargin<4
                    removeRotate3D=true;
                    if nargin<3
                        removeZoom=true;
                    end
                end
            end
            if first
                set(fig, 'Toolbar', 'figure');
                Gui.removeToolbarExcess(fig, removeRotate3D, removeEdit);
                if nargin>2 && removeZoom
                    Gui.removeForTooltip(fig, 'Zoom In');
                    Gui.removeForTooltip(fig, 'Zoom Out');
                    Gui.removeForTooltip(fig, 'Pan');
                end
                drawnow;
                this=ToolBar(findall(fig,'tag','FigureToolBar'));
            else
                ht_=uitoolbar(fig);
                drawnow;
                this=ToolBar(ht_);
            end
        end
    end

end


text:
function y = sum_square_abs( varargin )

%SUM_SQUARE   Sum of squares.
%   For vectors, SUM_SQUARE_ABS(X) is the sum of the squares of the
%   absolute values of the elements of the vector; i.e., SUM(ABS(X).^2).
%
%   For matrices, SUM_SQUARE_ABS(X) is a row vector containing the
%   application of SUM_SQUARE_ABS to each column. For N-D arrays, the
%   SUM_SQUARE_ABS operation is applied to the first non-singleton
%   dimension of X.
%
%   SUM_SQUARE(X,DIM) takes the sum along the dimension DIM of X.
%
%   Disciplined convex programming information:
%       If X is real, then SUM_SQUARE(X,...) is convex and nonmonotonic in
%       X. If X is complex, then SUM_SQUARE(X,...) is neither convex nor
%       concave. Thus, when used in CVX expressions, X must be affine. DIM
%       must be constant.

try
    varargin{end+1:2} = [];
    y = sum_square( varargin{:}, true );
catch exc
	cvx_throw( exc );
end

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
% A simple tutorial file to interface with RF
% Options copied from http://cran.r-project.org/web/packages/randomForest/randomForest.pdf

%run plethora of tests
clc
close all

%compile everything
if strcmpi(computer,'PCWIN') |strcmpi(computer,'PCWIN64')
   compile_windows
else
   compile_linux
end

total_train_time=0;
total_test_time=0;

%diabetes
load data/diabetes

%modify so that training data is NxD and labels are Nx1, where N=#of
%examples, D=# of features

X = diabetes.x;
Y = diabetes.y;

[N D] =size(X);
%randomly split into 400 examples for training and 42 for testing
randvector = randperm(N);

X_trn = X(randvector(1:400),:);
Y_trn = Y(randvector(1:400));
X_tst = X(randvector(401:end),:);
Y_tst = Y(randvector(401:end));

% example 1:  simply use with the defaults
    model = regRF_train(X_trn,Y_trn);
    Y_hat = regRF_predict(X_tst,model);
    fprintf('\nexample 1: MSE rate %f\n',   sum((Y_hat-Y_tst).^2));

% % example 2:  set to 100 trees
%     model = regRF_train(X_trn,Y_trn, 100);
%     Y_hat = regRF_predict(X_tst,model);
%     fprintf('\nexample 2: MSE rate %f\n',   sum((Y_hat-Y_tst).^2));
%
% % example 3:  set to 100 trees, mtry = 2
%     model = regRF_train(X_trn,Y_trn, 100,2);
%     Y_hat = regRF_predict(X_tst,model);
%     fprintf('\nexample 3: MSE rate %f\n',   sum((Y_hat-Y_tst).^2));
%
% % example 4:  set to defaults trees and mtry by specifying values as 0
%     model = regRF_train(X_trn,Y_trn, 0, 0);
%     Y_hat = regRF_predict(X_tst,model);
%     fprintf('\nexample 4: MSE rate %f\n',   sum((Y_hat-Y_tst).^2));
%
% % % example 5: set sampling without replacement (default is with replacement)
%     extra_options.replace = 0 ;
%     model = regRF_train(X_trn,Y_trn, 100, 4, extra_options);
%     Y_hat = regRF_predict(X_tst,model);
%     fprintf('\nexample 5: MSE rate %f\n',   sum((Y_hat-Y_tst).^2));
%
% % example 6: sampsize example
%     %  extra_options.sampsize =  Size(s) of sample to draw. For classification,
%     %                   if sampsize is a vector of the length the number of strata, then sampling is stratified by strata,
%     %                   and the elements of sampsize indicate the numbers to be drawn from the strata.
%     clear extra_options
%     extra_options.sampsize = size(X_trn,1)*2/3;
%
%     model = regRF_train(X_trn,Y_trn, 100, 4, extra_options);
%     Y_hat = regRF_predict(X_tst,model);
%     fprintf('\nexample 6: MSE rate %f\n',   sum((Y_hat-Y_tst).^2));
%
% % example 7: nodesize
%     %  extra_options.nodesize = Minimum size of terminal nodes. Setting this number larger causes smaller trees
%     %                   to be grown (and thus take less time). Note that the default values are different
%     %                   for classification (1) and regression (5).
%     clear extra_options
%     extra_options.nodesize = 7;
%
%     model = regRF_train(X_trn,Y_trn, 100, 4, extra_options);
%     Y_hat = regRF_predict(X_tst,model);
%     fprintf('\nexample 7: MSE rate %f\n',   sum((Y_hat-Y_tst).^2));
%
%
% % example 8: calculating importance
%     clear extra_options
%     extra_options.importance = 1; %(0 = (Default) Don't, 1=calculate)
%
%     model = regRF_train(X_trn,Y_trn, 100, 4, extra_options);
%     Y_hat = regRF_predict(X_tst,model);
%     fprintf('\nexample 8: MSE rate %f\n',   sum((Y_hat-Y_tst).^2));
%
%     %model will have 3 variables for importance importanceSD and localImp
%     %importance = a matrix with nclass + 2 (for classification) or two (for regression) columns.
%     %           For classification, the first nclass columns are the class-specific measures
%     %           computed as mean decrease in accuracy. The nclass + 1st column is the
%     %           mean decrease in accuracy over all classes. The last column is the mean decrease
%     %           in Gini index. For Regression, the first column is the mean decrease in
%     %           accuracy and the second the mean decrease in MSE. If importance=FALSE,
%     %           the last measure is still returned as a vector.
%     figure('Name','Importance Plots')
%     subplot(3,1,1);
%     bar(model.importance(:,end-1));xlabel('feature');ylabel('magnitude');
%     title('Mean decrease in Accuracy');
%
%     subplot(3,1,2);
%     bar(model.importance(:,end));xlabel('feature');ylabel('magnitude');
%     title('Mean decrease in Gini index');
%
%
%     %importanceSD = The ?standard errors? of the permutation-based importance measure. For classification,
%     %           a D by nclass + 1 matrix corresponding to the first nclass + 1
%     %           columns of the importance matrix. For regression, a length p vector.
%     model.importanceSD
%     subplot(3,1,3);
%     bar(model.importanceSD);xlabel('feature');ylabel('magnitude');
%     title('Std. errors of importance measure');
%
% % example 9: calculating local importance
%     %  extra_options.localImp = Should casewise importance measure be computed? (Setting this to TRUE will
%     %                   override importance.)
%     %localImp  = a D by N matrix containing the casewise importance measures, the [i,j] element
%     %           of which is the importance of i-th variable on the j-th case. NULL if
%     %          localImp=FALSE.
%     clear extra_options
%     extra_options.localImp = 1; %(0 = (Default) Don't, 1=calculate)
%
%     model = regRF_train(X_trn,Y_trn, 100, 4, extra_options);
%     Y_hat = regRF_predict(X_tst,model);
%     fprintf('\nexample 9: MSE rate %f\n',   sum((Y_hat-Y_tst).^2));
%
%     model.localImp
%
% % example 10: calculating proximity
%     %  extra_options.proximity = Should proximity measure among the rows be calculated?
%     clear extra_options
%     extra_options.proximity = 1; %(0 = (Default) Don't, 1=calculate)
%
%     model = regRF_train(X_trn,Y_trn, 100, 4, extra_options);
%     Y_hat = regRF_predict(X_tst,model);
%     fprintf('\nexample 10: MSE rate %f\n',   sum((Y_hat-Y_tst).^2));
%
%     model.proximity
%
%
% % example 11: use only OOB for proximity
%     %  extra_options.oob_prox = Should proximity be calculated only on 'out-of-bag' data?
%     clear extra_options
%     extra_options.proximity = 1; %(0 = (Default) Don't, 1=calculate)
%     extra_options.oob_prox = 0; %(Default = 1 if proximity is enabled,  Don't 0)
%
%     model = regRF_train(X_trn,Y_trn, 100, 4, extra_options);
%     Y_hat = regRF_predict(X_tst,model);
%     fprintf('\nexample 11: MSE rate %f\n',   sum((Y_hat-Y_tst).^2));
%
%
% % example 12: to see what is going on behind the scenes
% %  extra_options.do_trace = If set to TRUE, give a more verbose output as randomForest is run. If set to
% %                   some integer, then running output is printed for every
% %                   do_trace trees.
%     clear extra_options
%     extra_options.do_trace = 1; %(Default = 0)
%
%     model = regRF_train(X_trn,Y_trn, 100, 4, extra_options);
%     Y_hat = regRF_predict(X_tst,model);
%     fprintf('\nexample 12: MSE rate %f\n',   sum((Y_hat-Y_tst).^2));
%
% % example 13: to see what is going on behind the scenes
% %  extra_options.keep_inbag Should an n by ntree matrix be returned that keeps track of which samples are
% %                   'in-bag' in which trees (but not how many times, if sampling with replacement)
%
%     clear extra_options
%     extra_options.keep_inbag = 1; %(Default = 0)
%
%     model = regRF_train(X_trn,Y_trn, 100, 4, extra_options);
%     Y_hat = regRF_predict(X_tst,model);
%     fprintf('\nexample 13: MSE rate %f\n',   sum((Y_hat-Y_tst).^2));
%
%     model.inbag
%
%
% % example 14: getting the OOB MSE rate. model will have mse field
%     model = regRF_train(X_trn,Y_trn);
%     Y_hat = regRF_predict(X_tst,model);
%     fprintf('\nexample 14: MSE rate %f\n',   sum((Y_hat-Y_tst).^2));
%
%     figure('Name','OOB error rate');
%     plot(model.mse); title('OOB MSE error rate');  xlabel('iteration (# trees)'); ylabel('OOB error rate');
%
% %
% % example 15: nPerm
% %               Number of times the OOB data are permuted per tree for assessing variable
% %               importance. Number larger than 1 gives slightly more stable estimate, but not
% %               very effective. Currently only implemented for regression.
%     clear extra_options
%     extra_options.importance=1;
%     extra_options.nPerm = 1; %(Default = 0)
%     model = regRF_train(X_trn,Y_trn,100,2,extra_options);
%     Y_hat = regRF_predict(X_tst,model);
%     fprintf('\nexample 15: MSE rate %f\n',   sum((Y_hat-Y_tst).^2));
%
%     figure('Name','Importance Plots nPerm=1')
%     subplot(2,1,1);
%     bar(model.importance(:,end-1));xlabel('feature');ylabel('magnitude');
%     title('Mean decrease in Accuracy');
%
%     subplot(2,1,2);
%     bar(model.importance(:,end));xlabel('feature');ylabel('magnitude');
%     title('Mean decrease in Gini index');
%
%     %let's now run with nPerm=3
%     clear extra_options
%     extra_options.importance=1;
%     extra_options.nPerm = 3; %(Default = 0)
%     model = regRF_train(X_trn,Y_trn,100,2,extra_options);
%     Y_hat = regRF_predict(X_tst,model);
%     fprintf('\nexample 15: MSE rate %f\n',   sum((Y_hat-Y_tst).^2));
%
%     figure('Name','Importance Plots nPerm=3')
%     subplot(2,1,1);
%     bar(model.importance(:,end-1));xlabel('feature');ylabel('magnitude');
%     title('Mean decrease in Accuracy');
%
%     subplot(2,1,2);
%     bar(model.importance(:,end));xlabel('feature');ylabel('magnitude');
%     title('Mean decrease in Gini index');
%
% % example 16: corr_bias (not recommended to use)
%     clear extra_options
%     extra_options.corr_bias=1;
%     model = regRF_train(X_trn,Y_trn,100,2,extra_options);
%     Y_hat = regRF_predict(X_tst,model);
%     fprintf('\nexample 16: MSE rate %f\n',   sum((Y_hat-Y_tst).^2));

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
% Fig. 5.9   Feedback Control of Dynamic Systems, 6e
%             Franklin, Powell, Emami
% Script to generate Figure 5.9

clf
n=1;
d=[1 8 32 0];
rlocus(n,d)
axis([-10 6 -6 6])
hold on
x=[ 0 -2 0 -2];
y=[0  2*sqrt(3)  0 -2*sqrt(3) ];
plot(x,y)
r=roots([1 4 16]);
plot(r,'*')
title('Fig. 5.9 Locus for L=1/[s(s+4)^2+16)]')
z=0:.1:.9;
 wn=1:1:10;
 sgrid(z, wn)
 hold off

author: Sable

repo: mcbench-benchmarks


text:
function [ rectsLTRB ] = RectLTWH2LTRB(rectsLTWH)
%rects (l, t, r, b) to (l, t, w, h)

rectsLTRB = [rectsLTWH(:, 1), rectsLTWH(:, 2), rectsLTWH(:, 1)+rectsLTWH(:,3)-1, rectsLTWH(:,2)+rectsLTWH(:,4)-1];
end

author: ShaoqingRen


text:
% Copyright 2019 Jonas Koenemann, Moritz Diehl, University of Freiburg
% Redistribution is permitted under the 3-Clause BSD License terms. Please
% ensure the above copyright notice is visible in any derived work.
%
function warning(msg)
  warning(msg)
  global oclHasWarnings
  oclHasWarnings = true;

author: OpenOCL

repo: OpenOCL

save_path: github-repos/MATLAB/OpenOCL-OpenOCL


text:
report_this_filefun(mfilename('fullpath'));

% drap a colormap onto topography
l = isnan(tmap);
tmap(l) = 1;

%l = tmap< 0.1;
%tmap(l) = 0;

[lat,lon] = meshgrat(tmap,tmapleg);
%[smap,smapleg] = country2mtx('switzerland',100);
%[lat0, lon0] = meshgrat(smap,smapleg);

%   % tmap = km2deg(tmap/1);
%   [X , Y]  = meshgrid(gx,gy);
%
%   %sw = interp2(lon0,lat0,smap,lon,lat);
%
%
%
%   ren = interp2(X,Y,re4,lon,lat);
%
%   mi = min(min(ren));
%   l =  isnan(ren);
%   ren(l) = mi-20;

%figure_w_normalized_uicontrolunits('pos',[150 500 1000 700])

exfig=figure_exists('Mcomp worldmap',0);
if exfig==0
    ui=figure_w_normalized_uicontrolunits( 'Name','Mcomp worldmap',...
        'NumberTitle','off');
    figure_w_normalized_uicontrolunits(ui)
end

clf
hold on; axis off
axesm('MapProjection','robinson')

%axesm('MapProjection','eqacylin','MapParallels',[],...
%   'MapLatLimit',[s4_south s3_north],'MapLonLimit',[s2_west s1_east])

ll = tmap < 0 & ren < 0;
ren(ll) = ren(ll)*0 + 20;
meshm(ren,tmapleg,size(tmap),tmap);

daspectm('m',30);
tightmap
view([0 90])
camlight; lighting phong
set(gca,'projection','perspective');

%   load worldlo
%   h = displaym(POline); set(h(1),'color',[0.9 0.9 0.9],'Linewidth',1.7)
%   h2 = displaym(PPpoint);
%   h = displaym(PPtext); trimcart(h);

% pl = plotm(ma(:,2),ma(:,1),'hw');
%  set(pl,'LineWidth',1.5,'MarkerSize',12,...
% 'MarkerFaceColor','w','MarkerEdgeColor','k')
%load coast.mat
%c = coast;
% plotm(c(:,1), c(:,2),'k','Linewidth',1);
% zdatam(handlem('allline'),10000) % keep line on surface
%zdatam(handlem('alltext'),10000) % keep line on surface

j = jet;
%j = j(64:-1:1,:);
%j = [ [ 0.9 0.9 0.9 ] ; j];
j = [ [ 0.9 0.9 0.9 ] ; j; [ 0.5 0.5 0.5] ];
%j = [ [ 0.9 0.9 0.9 ] ; j ];

caxis([ min(min(re4)) max(max(re4)) ]);
caxis([3.98 5.2])
colormap(j); brighten(0.1);

axis off; set(gcf,'color','k')

setm(gca,'ffacecolor','k')
setm(gca,'fedgecolor','w','flinewidth',1);

setm(gca,'mlabellocation',60)
setm(gca,'meridianlabel','on')
setm(gca,'plabellocation',15)
setm(gca,'parallellabel','on')
setm(gca,'Fontcolor','w','FontSize',7,'Fontweight','bold','Labelunits','degrees')
setm(gca,'Labelformat','none')
setm(gca,'Grid','on')
setm(gca,'GLineStyle','-')

h5 = colorbar;
set(h5,'position',[0.8 0.35 0.01 0.3],'TickDir','out','Ycolor','w','Xcolor','w',...
    'Fontweight','bold','FontSize',7);
set(gcf,'Inverthardcopy','off');


text:
function poly = edgeToPolyline(edge, N)
%EDGETOPOLYLINE Convert an edge to a polyline with a given number of segments
%
%   POLY = edgeToPolyline(EDGE, N)
%
%   Example
%     edge = [10 20 60 40];
%     poly = edgeToPolyline(edge, 10);
%     drawEdge(edge, 'lineWidth', 2);
%     hold on
%     drawPoint(poly);
%     axis equal;
%
%   See also
%     edges2d, drawEdge, drawPolyline
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2011-11-25,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

if N < 1
    error('number of segments must be greater than 1');
end

if length(edge) == 4
    % case of planar edges
    p1 = edge(1:2);
    p2 = edge(3:4);
    poly = [linspace(p1(1), p2(1), N+1)' linspace(p1(2), p2(2), N+1)'];

else
    % case of 3D edges
    p1 = edge(1:3);
    p2 = edge(4:6);
    poly = [...
        linspace(p1(1), p2(1), N+1)' ...
        linspace(p1(2), p2(2), N+1)' ...
        linspace(p1(3), p2(3), N+1)'];
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function [dpsi, deps] = nut2000a (date1, date2)

% nutation based on iau 2000a theory

% input

%  date1, date2 = tt julian date
%  (julian date = date1 + date2)

% output

%  dpsi = nutation in longitude in radians

%  deps = nutation in obliquity in radians

% reference

%  Nutation Series Evaluation in NOVAS 3.0
%  USNO Circular No. 181, December 15, 2009

% ported from NOVAS 3.0 Fortran subroutine

% Orbital Mechanics with MATLAB

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

global jplephem_inutate nals napl icpl cls

if (jplephem_inutate == 1)

    % read data files

    nals = csvread('nals.csv');

    napl = csvread('napl.csv');

    icpl = csvread('icpl.csv');

    cls = csvread('cls.csv');

    % transpose matrices

    nals = nals';

    napl = napl';

    icpl = icpl';

    cls = cls';

    % reset flag

    jplephem_inutate = 0;

end

% arc seconds to radians

das2r = 4.848136811095359935899141d-6;

% arc seconds in a full circle

turnas = 1296000.0d0;

% 2 * pi

d2pi = 6.283185307179586476925287d0;

% units of 0.1 microarcsecond to radians

u2r = das2r / 1.0d7;

% reference epoch (j2000)

dj0 = 2451545.0d0;

% days per julian century

djc = 36525.0d0;

% -------------------------
% luni-solar nutation model
% -------------------------

%  ---------------
%  planetary terms
%  ---------------

% interval between fundamental epoch j2000.0 and given date (jc)

t = ((date1 - dj0) + date2) / djc;

% -------------------
% luni-solar nutation
% -------------------

% fundamental (delaunay) arguments from simon et al. (1994)

% mean anomaly of the moon (radians)

el = mod (485868.249036d0 + t * (1717915923.2178d0 ...
    + t * (31.8792d0 + t * (0.051635d0 ...
    + t * (-0.00024470d0)))), turnas) * das2r;

% mean anomaly of the sun (radians)

elp = mod (1287104.79305d0 + t * (129596581.0481d0 ...
    + t * (-0.5532d0 + t * (0.000136d0 ...
    + t * (-0.00001149d0)))), turnas) * das2r;

% mean argument of the latitude of the moon (radians)

f = mod (335779.526232d0 + t * (1739527262.8478d0 ...
    + t * (-12.7512d0 + t * (-0.001037d0 ...
    + t * (0.00000417d0)))), turnas) * das2r;

% mean elongation of the moon from the sun (radians)

d = mod (1072260.70369d0 + t * (1602961601.2090d0 ...
    + t * (-6.3706d0 + t * (0.006593d0 ...
    + t * (-0.00003169d0)))), turnas) * das2r;

% mean longitude of the ascending node of the moon (radians)

om = mod (450160.398036d0 + t * (-6962890.5431d0 ...
    + t * (7.4722d0 + t * (0.007702d0 ...
    + t * (-0.00005939d0)))), turnas) * das2r;

% summation of luni-solar nutation series (in reverse order)

arg = mod ((nals(1, :)) * el  + (nals(2, :)) * elp + (nals(3, :)) * f ...
    + (nals(4, :)) * d  + (nals(5, :)) * om, d2pi);

sarg = sin(arg);

carg = cos(arg);

dp = sum((cls(1,:) + cls(2,:) * t) .* sarg + cls(3,:) .* carg);

de = sum((cls(4,:) + cls(5,:) * t) .* carg + cls(6,:) .* sarg);

% convert from 0.1 microarcsec units to radians

dpsils = dp * u2r;

depsls = de * u2r;

% ------------------
% planetary nutation
% ------------------

% mean anomaly of the moon (radians)

al = mod (2.35555598d0 + 8328.6914269554d0 * t, d2pi);

% mean anomaly of the sun (radians)

alsu = mod (6.24006013d0 + 628.301955d0 * t, d2pi);

% mean argument of the latitude of the moon (radians)

af = mod (1.627905234d0 + 8433.466158131d0 * t, d2pi);

% mean elongation of the moon from the sun (radians)

ad = mod (5.198466741d0 + 7771.3771468121d0 * t, d2pi);

% mean longitude of the ascending node of the moon (radians)

aom = mod (2.18243920d0 - 33.757045d0 * t, d2pi);

% general accumulated precession in longitude (radians)

apa = (0.02438175d0 + 0.00000538691d0 * t) * t;

% planetary longitudes, mercury through neptune (souchay et al. 1999)

alme = mod (4.402608842d0 + 2608.7903141574d0 * t, d2pi);

alve = mod (3.176146697d0 + 1021.3285546211d0 * t, d2pi);

alea = mod (1.753470314d0 +  628.3075849991d0 * t, d2pi);

alma = mod (6.203480913d0 +  334.0612426700d0 * t, d2pi);

alju = mod (0.599546497d0 +   52.9690962641d0 * t, d2pi);

alsa = mod (0.874016757d0 +   21.3299104960d0 * t, d2pi);

alur = mod (5.481293871d0 +    7.4781598567d0 * t, d2pi);

alne = mod (5.321159000d0 +    3.8127774000d0 * t, d2pi);

% summation of planetary nutation series (in reverse order)

arg = mod ((napl(1, :)) * al + (napl(2, :)) * alsu + (napl(3, :)) * af ...
    + (napl(4, :)) * ad + (napl(5, :)) * aom  + (napl(6, :)) * alme ...
    + (napl(7, :)) * alve + (napl(8, :)) * alea + (napl(9, :)) * alma ...
    + (napl(10, :)) * alju + (napl(11, :)) * alsa + (napl(12, :)) * alur ...
    + (napl(13, :)) * alne + (napl(14, :)) * apa, d2pi);
sarg = sin(arg);

carg = cos(arg);

dp = sum((icpl(1,:)) .* sarg + (icpl(2,:)) .* carg);

de = sum((icpl(3,:)) .* sarg + (icpl(4,:)) .* carg);

% convert from 0.1 microarcsec units to radians

dpsipl = dp * u2r;

depspl = de * u2r;

% add planetary and luni-solar components

dpsi = dpsipl + dpsils;

deps = depspl + depsls;

end

author: goGPS-Project


text:
classdef SimpAllInterpolationExplicit < MaterialInterpolation

    methods  (Access = public)

        function obj = SimpAllInterpolationExplicit(cParams)
            obj.init(cParams);
            obj.computeNstre();
            obj.computeSymbolicInterpolationFunctions();
        end

    end

    methods (Access = protected)

        function [mS,dmS] = computeMuSymbolicFunctionAndDerivative(obj)
            mS  = obj.computeSymMu();
            dmS = diff(mS);
        end

        function [kS,dkS] = computeKappaSymbolicFunctionAndDerivative(obj)
            kS  = obj.computeSymKappa();
            dkS = diff(kS);
        end

        function mu = computeSymMu(obj)
            m0   = obj.matProp.mu0;
            m1   = obj.matProp.mu1;
            k0   = obj.matProp.kappa0;
            k1   = obj.matProp.kappa1;
            eta0 = obj.computeEtaMu(m0,k0);
            eta1 = obj.computeEtaMu(m1,k1);
            c    = obj.computeCoeff(m0,m1,eta0,eta1);
            mu   = obj.computeRationalFunction(c);
        end

        function kappa = computeSymKappa(obj)
            m0    = obj.matProp.mu0;
            m1    = obj.matProp.mu1;
            k0    = obj.matProp.kappa0;
            k1    = obj.matProp.kappa1;
            eta0  = obj.computeEtaKappa(m0);
            eta1  = obj.computeEtaKappa(m1);
            c     = obj.computeCoeff(k0,k1,eta0,eta1);
            kappa = obj.computeRationalFunction(c);
        end

        function etaMu = computeEtaMu(obj,mu,kappa)
            N   = obj.ndim;
            num = -mu*(4*mu - kappa*N^2 - 2*mu*N^2 + 2*mu*N);
            den = 2*N*(kappa + 2*mu);
            etaMu = num/den;
        end

        function etaKappa = computeEtaKappa(obj,mu)
            N   = obj.ndim;
            num = 2*mu*(N-1);
            den = N;
            etaKappa = num/den;
        end

    end

    methods (Access = protected, Static)

        function c = computeCoeff(f0,f1,eta0,eta1)
            c.n01 =  -(f1 - f0)*(eta1 - eta0);
            c.n0  =  f0*(f1 + eta0);
            c.n1  =  f1*(f0 + eta1);
            c.d0  =  (f1 + eta0);
            c.d1  =  (f0 + eta1);
        end

        function f = computeRationalFunction(s)
            rho = sym('rho','positive');
            n01 = s.n01;
            n0  = s.n0;
            n1  = s.n1;
            d0  = s.d0;
            d1  = s.d1;
            num = n01*(1-rho)*(rho) + n0*(1-rho) + n1*rho;
            den = d0*(1-rho) + d1*rho;
            f   = num/den;
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function plot(obj, varargin)
% PLOT Plot of point cloud.
% ------------------------------------------------------------------------------
% DESCRIPTION/NOTES
% * Only active points are visualized.
% ------------------------------------------------------------------------------
% INPUT
% 1 ['Color', color]
%   Color of the points. The following choices are possible:
%
%     * Unicolor (format = 'ColorName'):
%       A single color is used for all points, e.g. 'yellow'. For possible
%       values, look for 'ColorSpec' in the documentation.
%
%     * Attribute (format = 'A.AttributeName'):
%       Coloring according to an attribute, e.g. 'A.roughness'.
%       Additionally, the following choices are possible:
%         * 'A.x'  : coloring according to x coordinate
%         * 'A.y'  : coloring according to y coordinate
%         * 'A.z'  : coloring according to z coordinate
%         * 'A.rgb': true color plot using the attributes r, g and b.
%
% 2 ['MarkerSize', markerSize]
%   Marker size of points.
%
% 3 ['MaxPoints', maxPoints]
%   Maximum number of points to plot. If the point cloud contains more than
%   maxPoi points, the points for plotting are randomly sampled among all
%   points.
%
% 4 ['ColormapName', colormapName]
%   Name of colormap defined as char (e.g. 'jet'). For possible values, look for
%   'colormap' in the documentation. Additionally, for visualization of classes,
%   the colormap 'classpal' can be used.
%   (This parameter is only considered for attribute dependent plots.)
%
% 5 ['CAxisLim', cAxisLim]
%   Limits of colorbar defined as vector with 2 elements.
%   (This parameter is only considered for attribute dependent plots.)
% ------------------------------------------------------------------------------
% EXAMPLES
% 1 Import a point cloud and plot points according to z coordinate (default).
%   pc = pointCloud('Lion.xyz');
%   pc.plot;
%
% 2 Unicolor plot.
%   pc = pointCloud('Lion.xyz');
%   pc.plot('Color', 'red');
%
% 3 Plot all points of a point cloud.
%   pc = pointCloud('Lion.xyz');
%   pc.plot('MaxPoints', Inf);
%
% 4 Plot only 1000 points of a point cloud with increased point size.
%   pc = pointCloud('Lion.xyz');
%   pc.plot('MaxPoints', 1000, 'MarkerSize', 10);
%
% 5 Plot points according to true color attributes r, g and b.
%   pc = pointCloud('Gieszkanne.xyz', 'Attributes', {'r' 'g' 'b'});
%   pc.plot('Color', 'A.rgb', 'MarkerSize', 5);
%
% 6 Plot points colored by an imported attribute.
%   pc = pointCloud('Lion.xyz', 'Attributes', {'nx' 'ny' 'nz' 'roughness'});
%   % Note: the imported attributes are saved as fields in the structure pc.A,
%   % e.g. the roughness is saved in pc.A.roughness.
%   pc.plot('Color', 'A.roughness');
% ------------------------------------------------------------------------------
% philipp.glira@gmail.com
% ------------------------------------------------------------------------------

% Start
[hFig, idxPC, firstCall] = plot_start;

if firstCall

    % Save colors to figure
    plot_saveColors(hFig)

    % Save default options to figure
    plot_saveOptions(hFig)

    % Create toolbar
    plot_createToolbar(hFig)

    % Set properties of figure
    plot_setFigureProperties(hFig)

end

% Save PC data to figure
plot_savePC(hFig, obj, idxPC); clear obj

% Create menu
plot_createMenu(hFig, firstCall, idxPC)

% Set parameters and plot
plot_setPrmAndPlot(hFig, idxPC, varargin{:})

% Temporary solution
if firstCall, set(gca, 'Visible', 'on'); end

end

author: pglira


text:
filename='BikeTriangle';
ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance'};
weights = [1 0.1];
constraint = {'volumeConstraint'};
optimizer = 'DualNestedInPrimal';
optimizerUnconstrained = 'SLERP';
incrementFactor = 1;
designVariable = 'LevelSet';
filterType = 'P1';
constraint_case = 'EQUALITY'; %'INEQUALITY';
line_search_initiator = 'STANDARD';
showBC = true;

nsteps = 50;
Vfrac_final = 0.2;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-3;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

monitoring_interval = 1;
optimalityInitial = 1e-3;
printing = false;
maxiter = 300;


text:
% APM Solver for simulation, estimation, and optimization with both
%  static (steady-state) and dynamic models. The dynamic modes can solve
%  index 2+ DAEs without numerical differentiation.
%
% y = apm_solve(app,imode)
%
% Function apm_solve uploads the model file (apm) and optionally
%   a data file (csv) with the same name to the web-server and performs
%   a forward-time stepping integration of ODE or DAE equations
%   with the following arguments:
%
%  Input:      app = model (apm) and data file (csv) name
%            imode = simulation mode {1..7}
%                               steady-state  dynamic  sequential
%                    simulate     1             4        7
%                    estimate     2             5        8 (under dev)
%                    optimize     3             6        9 (under dev)
%
% Output: y.names  = names of all variables
%         y.values = tables of values corresponding to y.names
%         y.nvar   = number of variables
%         y.x      = combined variables and values but variable
%                      names may be modified to make them valid MATLAB
%                      characters (e.g. replace '[' with '')
%
function y = apm_solve(app,imode)
    % check for number of arguments
    if nargin==0,
        disp('Please specify the APM model name such as:')
        disp('  y = apm_solve(myModel)')
        y = [];
        return
    end

    % server and application file names
    server = 'http://byu.apmonitor.com';
    app_model = [app '.apm'];
    app_data = [app '.csv'];

    % only one input argument (imode not specified)
    if nargin==1,
        if (exist(app_data,'file')),
            disp('Using default imode of 7 (sequential simulation)')
            imode = 7;
        else
            disp('Using default imode of 3 (steady-state optimization)')
            imode = 3;
        end
    end

    % application name
	app = lower(deblank(app));

    % clear previous application
    apm(server,app,'clear all');

    % check that model file exists (required)
    if (~exist(app_model,'file')),
        disp(['Error: file ' app_model ' does not exist']);
        y = [];
        return
    else
        % load model file
        apm_load(server,app,app_model);
    end

    % check if data file exists (optional)
    if (exist(app_data,'file')),
        % load data file
        csv_load(server,app,app_data);
    end

    % default options
    % use or don't use web viewer
    web = false;
    if web,
        apm_option(server,app,'nlc.web',2);
    else
        apm_option(server,app,'nlc.web',0);
    end
    % internal nodes in the collocation (between 2 and 6)
    apm_option(server,app,'nlc.nodes',3);
    % sensitivity analysis (default: 0 - off)
    apm_option(server,app,'nlc.sensitivity',0);
    % simulation mode (1=ss,  2=mpu, 3=rto)
    %                 (4=sim, 5=est, 6=nlc, 7=sqs)
    apm_option(server,app,'nlc.imode',imode);

    % attempt solution
    solver_output = apm(server,app,'solve');

    % check for successful solution
    status = apm_tag(server,app,'nlc.appstatus');

    if status==1,
        % open web viewer if selected
        if web,
            apm_web(server,app);
        end
        % retrieve solution and solution.csv
        y = apm_sol(server,app);
        return
    else
        apm_web_root(server,app);
        disp(solver_output);
        disp('Error: Did not converge to a solution');
        y = solver_output;
        return
    end

author: APMonitor

repo: arduino

save_path: github-repos/MATLAB/APMonitor-arduino


text:
function [MW] = Btuph2MW(Btuph)
% Convert power from British thermal units per hour to megawatts.
% Chad A. Greene 2012
MW = Btuph*2.930710702e-7;


text:
function mrk= mrk_defineClasses(mk, classDef, varargin)
%MRK_DEFINECLASSES - Convert BV Markers to class defining marker structure
%
%Synopsis:
% MRK_OUT= mrk_defineClasses(MRK_IN, CLASS_DEF, <OPT>)
%
%Arguments:
% MRK_IN: Marker structure as received by file_readBVmarkers.
%  MarkerFormat can be numeric or string.
% CLASS_DEF: Class array of size {2 x nClasses}. The first row
%  specifies the markers of each class, each cell being either
%  a cell array of strings (for MarkerFormat='string') or a vector
%  of integers (for MarkerFormat='numeric').
%  The second row specifies the class names, each cell begin a string.
%  (If the second row does not exist, generic class names are defined.)
% OPT: struct or property/value list of optional properties:
%  'RemoveVoidClasses': Void classes are removed from the list of classes,
%     default 0.
%  'KeepAllMarkers': Keep also for markers, which do not belong to any
%     of the specified classes, default 0.
%
%Returns:
% MRK_OUT: Marker structure with classes defined by labels
%     (fields 'y' and 'className')
%
%Example:
% file= 'Gabriel_01_07_24/selfpaced1sGabriel';
% [cnt,mk]= file_readBV(file);
% classDef= {[65 70], [74 192]; 'left','right'};
% mrk= mrk_defineClasses(mk, classDef);
%
% [cnt,mk]= file_readBV(file, 'MarkerFormat','string');
% classDef= {{'S 65','S 70'},{'S 74', 'S192'}; 'left','right'}
% mrk= mrk_defineClasses(mk, classDef);
% %% does the same

props= {'KeepAllMarkers'     0  'BOOL'
        'RemoveVoidClasses'  0  'BOOL'};
props_selectEvents= mrk_selectEvents;

if nargin==0,
  mrk= opt_catProps(props, props_selectEvents);
  return
end

opt= opt_proplistToStruct(varargin{:});
opt= opt_setDefaults(opt, props);
opt_checkProplist(opt, props);

misc_checkType(mk, 'STRUCT(time event)');
misc_checkType(mk.event, 'STRUCT(desc)', 'mk.event');
misc_checkType(classDef, 'CELL');

nClasses= size(classDef,2);
if iscell(mk.event.desc),
  mrk= struct('time', mk.time, 'event',struct('desc',{mk.event.desc}));
else
  mrk= struct('time', mk.time, 'event',struct('desc',mk.event.desc));
end

mrk.y= zeros(nClasses, numel(mrk.time));
for cc= 1:nClasses,
  if isnumeric(classDef{1,cc}),
    % vector as in {[10 11], [20:26];  'target', 'nontarget'}
    mrk.y(cc,:)= ismember(mk.event.desc, classDef{1,cc});
  elseif iscell(classDef{1,cc}),
    % cell of strings as in {{'S 10','S 11'}, {'S 20','S 21'};
    %                        'target',        'nontarget'}
    mrk.y(cc,:)= ismember(mk.event.desc, classDef{1,cc});
  else
    % single string as in {'S10', 'S20';  'target', 'nontarget'}
    mrk.y(cc,:)= ismember(mk.event.desc, classDef(1,cc));
  end
end

if size(classDef,1)>1,
  mrk.className= classDef(2,:);
else
  mrk.className= str_cprintf('class %d', 1:nClasses);
end

if ~opt.KeepAllMarkers,
  opt_selectEvents= opt_substruct(opt, props_selectEvents(:,1));
  mrk= mrk_selectEvents(mrk, 'valid', opt_selectEvents);
end
if opt.RemoveVoidClasses,
  mrk= mrk_removeVoidClasses(mrk);
end

author: bbci


text:
%
% y = nnormn(x,dim,p)
%
% NNORMN normalizes an array x by its p-vector norms along dimension <dim>.
%
% dim: dimension along which to calculate norm.  Default first nonsingleton
%   p: norm-type.  Default is 2.
%
% Equivalence: normc(x) == nnormn(x,1,2), normr(x) == nnormn(x,2,2)
%
% See also NORMC, NORMR, NNORM

% Created by Bill Winter December 2005
% Based on normc and normr
function x = nnormn(x,dim,p)
siz = size(x);
if nargin < 2, dim = find(siz > 1,1); end
if nargin < 3, p = 2; end
switch p
    case inf,   a = max(x,[],dim);              % max
    case -inf,  a = min(x,[],dim);              % min
    case 1,     a = sum(abs(x),dim);            % manhattan
    case 2,     a = sqrt(sum(abs(x).^2,dim));   % euclidean
    otherwise,  a = sum(abs(x).^p,dim).^(1/p);  % p-norm
end
a(a == 0) = 1;
N(1:length(siz)) = {':'};
N{dim} = ones(1,siz(dim));
x = x./a(N{:});


text:
function fileOffsets = dtiAppendPathwaysToPDB(fg, filename)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% dtiAppendDatabasePathways(toAppend, filename)
%
% toAppend: the pathway database we want to append to this file
% fileOffsets: an array of file offsets for each pathway.
%
% Author: DA
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

numstats = 0; % xxx
algoType = 0; % xxx

fid = fopen (filename, 'rb');
% figure out where num paths is, read it in
offset = fread (fid, 1, 'uint');
%fprintf ('Seeking to %d\n', offset);
fseek (fid, offset, -1);
[oldNumPaths, count] = fread (fid, 1, 'uint');
fclose (fid);
if count == 0
    %fprintf ('First time!');
    oldNumPaths = 0;
end
newNumPaths = oldNumPaths + length(fg.fibers);
fid = fopen (filename, 'r+b');
fseek (fid, offset, -1);
fwrite (fid, newNumPaths, 'uint');
fclose(fid);
fid = fopen (filename, 'ab');
fseek (fid, 0, 1);
fileOffsets = zeros(length(fg.fibers),1);

for p = 1:length(fg.fibers)
    fileOffsets(p) = ftell(fid);
    path_offset = 3 * 4  + numstats*8; % 3 int, numstats double
    fwrite(fid,path_offset,'int');

    nn = size(fg.fibers{p}, 2);

    fwrite(fid,nn,'int');
    fwrite(fid,algoType,'int'); % algo type
    fwrite(fid,0,'int'); % seedpointindex ???

    fwrite(fid,fg.fibers{p}(:), 'double');
    % Stats
    %for as = 1:numstats
     %   fwrite(fid,toAppend.pathways(p).path_stat_vector(as),'double');
    %end

    % Writing path nodes
    %fwrite(fid,[holdingchains(p).xpos*mmPerVox(1); holdingchains(p).ypos*mmPerVox(2); holdingchains(p).zpos*mmPerVox(3)],'double');
    %pos = asMatrixStruct(toAppend.pathways(p));
    % Change from 0 based to 1 based positions
    %pos = pos - repmat(toAppend.mm_scale(:),1,size(pos,2));
    %fwrite(fid,pos,'double');

    % Writing stats values per position
    %for as = 1:numstats
    %    if( toAppend.pathway_statistic_headers(as).is_computed_per_point )
    %        fwrite(fid,toAppend.pathways(p).point_stat_array(as,:),'double');
    %    end
    %end
end

%disp('Appended pathways to database.');
fclose(fid);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% ASTROTIK by Francesco Santilli
% R2BP (Restricted Two Bodies Problem)
% f2H converts true anomaly in hyperbolic anomaly for an hyperbolic orbit.
%
% Usage: H = f2H(f,e)
%
% where: f(k) = true anomaly [rad]
%        e = eccentricity [-] (e>1)
%        H(k) = hyperbolic anomaly [rad]

function H = f2H(f,e)

    if ~(nargin == 2)
        error('Wrong number of input arguments.')
    end

    check(f,1)
    check(e,0)

    if e <= 1
        error('e must be strictly major than 1.')
    end

    ee = sqrt((e-1)/(e+1));
    H = 2*atanh( ee*tan(f/2) );

    kk = f/(2*pi);
    k = round(kk) - ((kk-fix(kk)) == 0.5);
    H = H + k*(2*pi);

end