text:
function [dat, roi_val, maskdat, beta_dat, beta_roi_val] = canlab_connectivity_preproc(dat, varargin)

% This function prepares data for connectivity analysis by removing nuisance
% variables and temporal filtering (high, low, or bandpass filter). This also
% can extract values from given masks and return averaged activity or pattern
% expression values.
%
% :Usage:
% ::
%
%    [preprocessed_dat, roi_val, maskdat] = canlab_connectivity_preproc(dat, varargin)
%    [preprocessed_dat, roi_val, maskdat, (optional) beta_dat, beta_roi_val] = canlab_connectivity_preproc(dat, varargin)
%
% :Features:
%
%   - can regress out nuisance variables with any additional nuisance matrix
%   - can remove signal from ventricle and white matter (default: use
%     canonical ventricle and white matter masks)
%   - can do temporal filtering, including high-pass, low-pass, or bandpass
%     filtering (it uses conn_filter.m from conn toolbox, which is using
%     fast fourier transform; see subfunction below)
%   - can extract data from given ROIs, and return averaged value or pattern
%     expression value (dot-product).
%   - Subsequent nuisance regression step is orthogonalized with respect to
%     pass-filter, consistent with principles in Lindquist 2019
%   - It is possible to extract ROI averages or linear patterns, save them,
%     and run this function afterwards. This will give the same result as
%     filtering each  voxel and then calculating the ROI average/pattern
%     response.
%   - (optional) can run additional GLM model with additional regressors if they are
%     specified. This runs the GLM *additionally*, therefore do not change
%     residualized image data or roi_vals. This GLM uses the same nuisance
%     matrix from above. output: beta_dat
%
% *Steps in order [with defaults]:*
%   1. Remove nuisance covariates (and linear trend if requested)
%   2. Remove ventricle and white matter - default: uses canonical images
%   3. High/low/bandpass filter on BOTH data and covs before regression (plot)
%   4. Residualization using regression (plot)
%   5. (optional) Run additional GLM using the same additional preprocessing
%   6. Windsorize based on distribution of full data matrix (plot)
%   7. Extract region-by-region average ROI or pattern expression data
%
% ..
%     Author and copyright information:
%
%     Copyright (C) 2014  Wani Woo
%
%     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:
%
%   **dat:**
%        fmri_data object with data
%
%   **dat.covariate:**
%        basic nuisance matrix
%
% :Optional Inputs:
%
%   **additional_nuisance**
%        When you have additional nuisance variables that you want
%        regress out from the data, you can use this option. This
%        option should be followed by a nuisance matrix (or values).
%        The matrix should have the same number of rows with the
%        number of images.
%
%   **vw**
%        When you want to regress out signals from ventricle and
%        white matter, you can use this option. Default is using canonical
%        white matter and ventricle images (but eroded versions).
%
%        You can also choose what to use to remove ventricle and
%        white matter signal between raw data or top 5 PCA
%        components (default). You can just put 'raw' if you want to
%        use raw signal than PCA compoenents.
%        also see: canlab_extract_ventricle_wm_timeseries.m
%        canlab_create_wm_ventricle_masks.m)
%        - *Example:* 'vw','raw'
%
%   **windsorize:**
%        Windsorizing entire data matrix to k x STD.
%        - *Example:* 'windsorize', 5 (windsorize to 5 STD)
%
%   **linear_trend:**
%        This option will include the linear trend to nuisance variables.
%
%   **hpf', 'lpf', or 'bpf:**
%        This option will do temporal filtering.
%        - 'hpf': high pass filter. This option should be followed by
%                 the lower bound of the frequency (e.g., .01 Hz [= 100
%                 sec]), then TR
%        - 'lpf': low pass filter. This option should be followed by
%                 the upper bound of the frequency (e.g., .25 Hz [= 4 sec]), then TR.
%        - 'bpf': bandpass filter. This should be followed by lower
%                 and upper bounds of the frequency (e.g., [.01 .25]), then TR.
%                 In all cases, after the frequency value, you need to provide TR in sec.
%        - *Example:* 'hpf', .01, TR
%                     'bpf', [.01 .25], TR
%
%   **regressors:**
%        - This option will run additional GLM with the regressors specified
%          by the experimenter. The number of rows should be the number of
%          images, and the number of columns are the number of regressors.
%        - with extract_roi, the function will extract roi values from beta
%          images as well.
%        - Example: 'regressors', X (you can build your own design matrix,
%                                   e.g., hrf convolution using onset and
%                                   duration, using onsets2fmridesign.m)
%
%   **extract_roi:**
%        This option will extract data from ROIs specified. This
%        option should be followed by one or more masks.
%        For one mask (potentially multiple ROIs, enter a char array with the mask name.
%        For multiple masks (1 or more), enter in a cell array of mask names.
%        You can specify methods with 'roi_methods' option.
%        - 'average_over' (default): calculate averaged value across the ROIs.
%        - 'pattern_expression': calculate dot-products between
%                                pattern mask and data
%        - 'unique_mask_values' (default): will divide a mask into
%                                          multiple regions that have different discrete values.
%        - 'contiguous_regions': will divide a mask into multiple
%                                contiguous regions.
%        - 'whole': will do average_over or pattern_expression across
%                   all the voxels within the mask.
%        - *Example:* 'extract_roi', mask, 'contiguous_regions'
%                     'extract_roi', mask, 'pattern_expression'
%
%   **no_preproc:**
%        If you want to skip the preprocessing part, and want to
%        extract ROI values only, you can use this option.
%
%   **no_plots:**
%        To suppress plots, enter 'no_plots'
%
% :Outputs:
%
%  **preprocessed_dat:**
%        fmri_data object after removing nuisance variables and
%        filtering temporal confounds.
%
%  **roi_val:**
%        returns values extracted from ROIs in cell arrays (if there are many different ROIs).
%        Each cell will have roi_val.dat, roi_val.mask_name, and roi_val.methods.
%
%  **maskdat:**
%        returns resampled roi (or pattern) masks
%
%  **beta_dat:**
%        returns beta images from the additional GLM (see 'regressors' option)
%
%
% :Examples:
% ::
%
%     roi_masks = which('weights_NSF_grouppred_cvpcr.img');
%     [preprocessed_dat, roi_val] = canlab_connectivity_preproc(dat, 'vw',
%               'bpf', [.008 .25], TR, 'extract_roi', roi_masks, 'pattern_expression');

% ..
%    PROGRAMMER'S NOTE
%    05/19/15 fixed a bug related to conn_filter
%    08/24/18 implemented simultaneous pass filter + regression
%    01/05/19: Tor and Marianne: Pass filter all covariates (even
%               spikes., Consistent with principles in Lindquist 2019
%               "sequential filtering" paper. See more notes below.
%    08/01/19: Wani reversed Tor and Marianne's revision above
%               after dicussion with Tor and others. Spikes should be removed
%               before "low-pass filtering" the nuisance covariates, and
%               added to the covariates later (see Ciric et al., 2018)
%               Lindquist et al., 2019 did not address the issue related to
%               the use of spike regressors combined with low pass
%               filtering. This issue should be further examined later.
%    08/01/19: Wani added an option for additional GLM
% ..

additional_R = [];
remove_vent_wm = false;
use_canonical_masks = true;
do_filter = false;
do_extract_roi = false;
do_preproc = true;
do_windsorize = false;
do_linear = false;
do_plots = true;
regressors = [];
do_additional_regress = false;

z = '===============================================================';
zz = '---------------------------------------------------------------';
fprintf('%s\n', z);

% parsing varargin
for i = 1:length(varargin)
    if ischar(varargin{i})
        switch varargin{i}
            case {'additional_nuisance', 'additional_R'}
                additional_R = varargin{i+1};

            case {'vw', 'ventricle_whitematter', 'remove_vw'}
                remove_vent_wm = true;

                % filter
            case {'hpf', 'lpf', 'bpf'}
                do_filter = true;
                TR = varargin{i+2};

                if any(strcmp(varargin, 'hpf')), filter(1) = varargin{i+1}; filter(2)= Inf; end
                if any(strcmp(varargin, 'lpf')), filter(1) = 0; filter(2) = varargin{i+1}; end
                if any(strcmp(varargin, 'bpf')), filter = varargin{i+1}; end

            case {'windsorize'}
                do_windsorize = true;
                k_std = varargin{i+1};

                % detrend
            case {'linear_trend'}
                do_linear = true;

                % extract roi values
            case {'extract_roi', 'roi', 'mask'}
                do_extract_roi = true;
                mask = varargin{i+1};

            case {'no_preproc'}
                do_preproc = false;

            case {'no_plots'}
                do_plots = false;

            case {'regressors'}
                regressors = varargin{i+1};
                do_additional_regress = true;
                regressors(:,all(regressors==1)) = []; % remove an intercept if any
                nreg = size(regressors,2);
        end
    end
end

%% Remove nuisance

if do_preproc

    tic;

    nobs = size(dat.dat, 2);
    if do_additional_regress
        if size(regressors,1) ~= nobs
            error('The number of rows of regressors is different from the number of images from dat. Please check.');
        end
    end

    fprintf('Canlab_connectivity_preproc: starting preprocessing ... \n');
    comb_R = [regressors dat.covariates additional_R]; % if all these are empty, comb_R should still be empty

    if do_additional_regress
        beta_dat = dat;
        beta_dat.dat = [];
    end

    dat.dat = double(dat.dat);  % enforce double to avoid data format problems
    dat.history{end+1} = 'Processed with canlab_connectivity_preproc';

    %% white matter, ventricle

    if do_plots
        create_figure('plot_qc',2,4);
        set(gcf, 'position', get(0,'screensize')); drawnow;
        plot_qc(dat,1,0,[]); drawnow;
    end

    if remove_vent_wm

        fprintf('::: removing ventricle & white matter signal ... \n');
        use_pca_comp_vw = true; % default

        if any(strcmp(varargin, 'raw')), use_pca_comp_vw = false; end

        if any(strcmp(varargin, 'datdir'))
            use_canonical_masks = false;
            subject_dir = varargin{find(strcmp(varargin, 'datdir'))+1};
        end

        if use_canonical_masks

            [value, components] = extract_gray_white_csf(dat, 'masks', ...
                {'gray_matter_mask.nii', 'canonical_white_matter_thrp5_ero1.nii', ...
                'canonical_ventricles_thrp5_ero1.nii'});
            wm_nuisance = value(:,2:3); % white, CSF
            wm_nuisance_comps = [components{2} components{3}]; % 5 components for white, CSF

        else

            if ~exist(fullfile(subject_dir, 'Structural/SPGR/white_matter.img'), 'file')...
                    && ~exist(fullfile(subject_dir, 'Structural/SPGR/ventricles.img'), 'file')
                wm_mask = filenames(fullfile(subject_dir, 'Structural/SPGR/wc2*.nii'), 'char', 'absolute');
                gm_mask = filenames(fullfile(subject_dir, 'Structural/SPGR/wc1*.nii'), 'char', 'absolute');
                canlab_create_wm_ventricle_masks(wm_mask, gm_mask, 'wm_thr', .9, 'vent_thr', .9);
            end

            [wm_nuisance, wm_nuisance_comps] = canlab_extract_ventricle_wm_timeseries(fullfile(subject_dir, 'Structural/SPGR'), dat, 'noplot');

        end

        if use_pca_comp_vw

            comb_R = [comb_R scale(wm_nuisance_comps)];

        elseif ~use_pca_comp_vw

            comb_R = [comb_R scale(wm_nuisance)];

        end

        dat.history{end+1} = 'Extracting white matter and ventricle signals';

    end

    %% linear trend
    if do_linear

        comb_R(:,end+1) = scale((1:nobs)');
        dat.history{end+1} = '(added linear detrending to nuisance covs)';

    end

    % (optional) global signal % not implemented

    if do_filter

            fprintf('::: temporal filtering brain data... \n');
            % conn_filter does mean-center by default.
            % We need mean back in

            % Pass-filter brain data
            % -------------------------------------------
            dat.dat = conn_filter(TR, filter, dat.dat', 'full')' ...
                + repmat(mean(dat.dat,2), 1, size(dat.dat,2));

            dat.history{end+1} = 'Done temporal filtering brain data';

            if do_plots
                plot_qc(dat,2,1,[]); drawnow;
            end
    end

    %% get residuals

    if ~isempty(comb_R)
        % Do regression/residualization

        if size(comb_R, 1) ~= nobs
            error('Nuisance covariates and data matrix have different nunmbers of cases.');
        end

        %% temporal filter (low and high): low pass, high pass, bandpass
        % Simultaneous temporal filtering and regression

        if do_filter

            fprintf('::: temporal filtering covariates ... \n');
            % conn_filter does mean-center by default.
            % We need mean back in

            comb_R = filter_cov_wo_spikes(comb_R, TR, filter);

            dat.history{end+1} = 'Done temporal filtering covariates';

        end

        % add intercept
        X = [comb_R ones(nobs,1)];

        if do_additional_regress

            % regression
            beta_dat.dat = data.dat * pinv(X)';
            beta_dat.dat = beta_dat.dat(:,1:nreg);
            beta_dat.covariates = X;
            X(:,1:nreg) = [];
        end

        % residualize
        dat.dat = dat.dat - dat.dat * pinv(X)' * X'; % this is more efficient than (dat.dat' - X*pinv(X)*dat.dat')'
        dat.history{end+1} = 'Regressed out nuisance covariates';

        if do_plots
            clim = plot_qc(dat,3,1,[]); drawnow;
        end

    else  % no regression

    end

    %% Windsorize entire data matrix to k STD
    if do_windsorize
        dat = windsorize(dat, k_std); % history added within the function
        if do_plots
            plot_qc(dat,4,0,clim); drawnow;
        end
    end

    t = toc;
    fprintf('Done: data preprocessing in %4.1d sec.\n\n', t);
    fprintf('%s\n\n', zz);
end

%% get ROI values

if do_extract_roi
    tic;
    fprintf('Canlab_preproc: extracting ROI signal ... \n');

    % default
    do_pattern = false;
    do_whole = false;
    unique_mask_values = true;
    contiguous_regions = false;

    if any(strcmp(varargin, 'pattern_expression')), do_pattern = true; unique_mask_values = false; end
    if any(strcmp(varargin, 'whole')), do_whole = true; unique_mask_values = false; end
    if any(strcmp(varargin, 'contiguous_regions')), contiguous_regions = true; unique_mask_values = false; end

    if ~iscell(mask)
        mask_temp = mask; clear mask;
        mask = cell(size(mask_temp,1),1); % preallocation
        for i = 1:size(mask_temp,1)
            mask{i} = mask_temp(i,:);
        end
    end

    [roi_val, maskdat] = extract_roi_sub(dat, mask, do_pattern, do_whole, unique_mask_values, contiguous_regions);

    if do_additional_regress
        beta_roi_val = extract_roi_sub(beta_dat, mask, do_pattern, do_whole, unique_mask_values, contiguous_regions);
    end

    % output: roi_val, maskdat
    % input: mask, dat, do_pattern, do_whole, unique_mask_values, contiguous_regions

    t2 = toc;
    fprintf('Done in %4.1d sec.\n', t2);
    fprintf('%s\n', zz);

    if do_preproc
        fprintf('Total time: 4.1%d sec.\n', t+t2);
        fprintf('%s\n', z);
    end
else
    roi_val = {};
end

end

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

function [y,fy]=conn_filter(rt, filter, x, option)

% from conn toolbox
% http://www.nitrc.org/projects/conn/

if nargin < 4, option='full'; end

fy=fft(x,[],1);
f=(0:size(x,1)-1);
f=min(f,size(x,1)-f);

switch(lower(option))
    case 'full'

        idx=find(f<filter(1)*(rt*size(x,1))|f>=filter(2)*(rt*size(x,1)));
        %idx=idx(idx>1);
        fy(idx,:)=0;
        y=real(ifft(fy,[],1))*2*size(x,1)*(min(.5,filter(2)*rt)-max(0,filter(1)*rt))/max(1,size(x,1)-numel(idx));

    case 'partial'

        idx=find(f>=filter(1)*(rt*size(x,1))&f<filter(2)*(rt*size(x,1)));
        %if ~any(idx==1), idx=[1,idx]; end
        y=fy(idx,:);

end

end

function obj = windsorize(obj, k)

% Windsorize entire data matrix to k x STD

m = sum(sum(obj.dat)) ./ numel(obj.dat); % save memory space

s = k .* std(obj.dat(:));

whbad = obj.dat < m - s;
obj.dat(whbad) = m - s;
whbad2 = obj.dat > m + s;
obj.dat(whbad2) = m + s;

nbad = sum(whbad(:)) + sum(whbad2(:));
percbad = 100 .* nbad ./ numel(obj.dat);

obj.history{end+1} = sprintf('Windsorized data matrix to %d STD; adjusted %3.0f values, %3.1f%% of values', k, nbad, percbad);
disp(obj.history{end});
end

function clim = plot_qc(dat, n, use7sd, clim)

% dat
% n: steps
% use7sd: use -6sd~+6sd as a range
% clim: specify color limit

titles = {'RAW:before Conn preproc', 'Pass filter', 'Residualize (Nuisance, VentWM, linear trend, if selected)', 'Windsorize'};

subplot(2, 4, n);
plot(mean(dat.dat));
set(gca, 'xlim', [0 size(dat.dat,2)]);
xlabel('timeseries');
ylabel('Global mean');
title(titles{n});

subplot(2, 4, n+4);

if use7sd
    m = nanmean(dat.dat(:));
    sd = nanstd(dat.dat(:));
    clim = [m-7*sd m+7*sd];
end

if isempty(clim)
    imagesc(dat.dat);
else
    imagesc(dat.dat, clim);
end

colormap(gray);
set(gca, 'xlim', [0 size(dat.dat,2)]);
xlabel('timeseries');
ylabel('voxels');
colorbar('southoutside');

end

function R = filter_cov_wo_spikes(R, TR, filter)

% Filter the covariates
% Eliminate bias from sequential filtering operations
% -------------------------------------------
% 1/15/2019: Tor and Marianne: Pass filter all covariates (even
% spikes), because filter transforms data in a way that matches
% what is done to spike regressors. So even if spike regressors
% "look funny", they adjust for the effects of "spikes" on the
% time series data appropriately, accounting for filtering.
%
% 8/1/2019: Wani reversed this to the previous version after
% dicussion with Tor and others. Spikes should be removed
% before "low-pass filtering" the nuisance covariates, and
% added to the covariates later (see Ciric et al., 2018)

% pass-filter covs (without outlier indices)
outlier_idx = sum(R==1)==1;
outliers = R(:,outlier_idx);
R(:,outlier_idx) = [];

R = conn_filter(TR, filter, R, 'full') ...
    + repmat(mean(R), size(R,1), 1);

% add outlier indices back in
R = [R outliers];

end

function [roi_val, maskdat] = extract_roi_sub(dat, mask, do_pattern, do_whole, unique_mask_values, contiguous_regions)

% subfunction for extract rois

% copy over apply_mask
maskdat = cell(numel(mask),1);
for i = 1:numel(mask)
    maskdat{i} = fmri_data(mask{i});
    isdiff = compare_space(dat, maskdat{i});

    if isdiff == 1 || isdiff == 2 % diff space, not just diff voxels
        maskdat{i} = resample_space(maskdat{i}, dat, 'nearest');

        if length(maskdat{i}.removed_voxels) == maskdat{i}.volInfo.nvox
            disp('Warning: resample_space returned illegal length for removed voxels. Fixing...');
            maskdat{i}.removed_voxels = maskdat{i}.removed_voxels(maskdat{i}.volInfo.wh_inmask);
        end
    end
end

roi_val = cell(numel(mask),1); % preallocation

for i = 1:numel(mask)
    if do_whole
        % single, whole-brain pattern
        if ~do_pattern

            roi_obj = apply_mask(dat, maskdat{i});
            roi_obj = remove_empty(roi_obj);
            roi_val{i}.dat = nanmean(roi_obj.dat)';
            roi_val{i}.mask_name = mask{i};
            roi_val{i}.methods = 'averaged_over using a whole mask';

        elseif do_pattern

            roi_val{i}.dat = apply_mask(dat, maskdat{i}, 'pattern_expression', 'ignore_missing');
            roi_val{i}.mask_name = mask{i};
            roi_val{i}.methods = 'pattern expression using a whole mask';

        end
    else
        % separate patterns for local regions

        if ~do_pattern
            if unique_mask_values

                roi_obj = extract_roi_averages(dat, maskdat{i});
                roi_val{i}.dat = cat(2, roi_obj(:).dat);

                % extract_roi_averages does not preserve missing cases: remove
                if any(dat.removed_images)
                    roi_val{i}.dat = roi_val{i}.dat(~dat.removed_images, :);
                end

                roi_val{i}.mask_name = mask{i};
                roi_val{i}.methods = 'averaged_over using unique_mask_values';

            elseif contiguous_regions

                roi_obj = extract_roi_averages(dat, maskdat{i}, 'contiguous_regions');
                roi_val{i}.dat = cat(2, roi_obj(:).dat);

                % extract_roi_averages does not preserve missing cases: remove
                if any(dat.removed_images)
                    roi_val{i}.dat = roi_val{i}.dat(~dat.removed_images, :);
                end

                roi_val{i}.mask_name = mask{i};
                roi_val{i}.methods = 'averaged_over using contiguous_regions';

            else
                warning('Unknown ROI extraction option. Not extracting ROI data...');
            end
        else
            % local patterns within each region

            roi_obj = extract_roi_averages(dat, maskdat{i}, 'pattern_expression', 'contiguous_regions');
            roi_val{i}.dat = cat(2, roi_obj(:).dat);

            % extract_roi_averages does not preserve missing cases: remove
            if any(dat.removed_images)
                roi_val{i}.dat = roi_val{i}.dat(~dat.removed_images, :);
            end

            roi_val{i}.mask_name = mask{i};
            roi_val{i}.methods = 'pattern expression using contiguous_regions';
        end
    end
end

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function bc = specific_bc(xbd,ybd)
%backwardstep_bc   Reference problem 5.2 boundary condition
%   bc = specific_bc(xbd,ybd);
%   input
%          xbd          x boundary coordinate vector
%          ybd          y boundary coordinate vector
%
%   specifies streamfunction associated with Backward step flow
%   IFISS function: DJS; 6 March 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
bc=0*xbd;
k=find(xbd==-1); bc(k)=2*ybd(k).*ybd(k).*(1-2*ybd(k)/3);
k=find(ybd==1); bc(k)=2/3;
return

author: johannesgerer

repo: jburkardt-m


text:
% STD_MAKETRIALINFO - create trial information structure using the
%                       .epoch structure of EEGLAB datasets
%
% Usage:
%   >> [STUDY,trialinfo] = std_maketrialinfo(STUDY, ALLEEG);
%   >> [~,trialinfo]     = std_maketrialinfo([], EEG);
%
% Inputs:
%   STUDY      - EEGLAB STUDY set
%   ALLEEG     - vector of the EEG datasets included in the STUDY structure
%
% Inputs:
%   STUDY      - EEGLAB STUDY set updated. The fields which is created or
%                updated is STUDY.datasetinfo.trialinfo
%   trialinfo  - Trial information structure (cell array)
%
% Authors: Arnaud Delorme, SCCN/INC/UCSD, April 2010

% Copyright (C) Arnaud Delorme arno@ucsd.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function [STUDY,alltrialinfo] = std_maketrialinfo(STUDY, ALLEEG)

%% test if .epoch field exist in ALLEEG structure
epochfield = cellfun(@isempty, { ALLEEG.epoch });
if any(epochfield)
    fprintf('Warning: some datasets are continuous and trial information cannot be created. Please use [ALLEEG.trials] to see which ones.\n');
    return;
end

for iEEG = 1:length(ALLEEG)
    % fill in empty field values and fill in with values in the same epoch
    ALLEEG(iEEG) = eeg_uniformepochinfo(ALLEEG(iEEG));
end

%% check if conversion of event is necessary
ff = {};
flagConvert = true;
for index = 1:length(ALLEEG), ff = union(ff, fieldnames(ALLEEG(index).event)); end
for iField = 1:length(ff)

    fieldChar = zeros(1,length(ALLEEG))*NaN;
    for index = 1:length(ALLEEG)
        if isfield(ALLEEG(index).event, ff{iField})
            if ischar(ALLEEG(index).event(1).(ff{iField}))
                 fieldChar(index) = 1;
            else fieldChar(index) = 0;
            end
        end
    end
    if ~all(fieldChar(~isnan(fieldChar)) == 1) && ~all(fieldChar(~isnan(fieldChar)) == 0)
        % need conversion to char
        for index = 1:length(ALLEEG)
            if fieldChar(index) == 0
                if flagConvert, disp('Warning: converting some event fields to strings - this may be slow'); flagConvert = false; end
                for iEvent = 1:length(ALLEEG(index).event)
                    ALLEEG(index).event(iEvent).(ff{iField}) = num2str(ALLEEG(index).event(iEvent).(ff{iField}));
                end
            end
        end
    end
end

%% Make trial info
for index = 1:length(ALLEEG)
    tmpevent = ALLEEG(index).event;
    eventlat = abs(eeg_point2lat( [ tmpevent.latency ], [ tmpevent.epoch ], ALLEEG(index).srate, [ALLEEG(index).xmin ALLEEG(index).xmax]));
    events   = ALLEEG(index).event;
    ff = fieldnames(events);
    ff = setdiff_bc(ff, { 'latency' 'urevent' 'epoch' });
    trialinfo = [];

    % process time locking event fields
    % ---------------------------------
    indtle    = find(eventlat == 0);
    epochs    = [ events(indtle).epoch ];
    extractepoch = true;

    % Double checking and changing threshold
    if length(epochs) < ALLEEG(index).trials
        indtle    = find(eventlat < 0.02);
        epochs    = [ events(indtle).epoch ];
    end

    if length(epochs) ~= ALLEEG(index).trials
        % special case where there are not the same number of time-locking
        % event as there are data epochs
        if length(unique(epochs)) ~= ALLEEG(index).trials
            extractepoch = false;
            disp('std_maketrialinfo: not the same number of time-locking events as trials, trial info ignored');
        else
            % pick one event per epoch
            [tmp tmpind] = unique_bc(epochs(end:-1:1)); % reversing the array ensures the first event gets picked
            tmpind = length(epochs)+1-tmpind;
            indtle = indtle(tmpind);
            if length(indtle) ~= ALLEEG(index).trials
                extractepoch = false;
                disp('std_maketrialinfo: not the same number of time-locking events as trials, trial info ignored');
            end
        end
    end
    if extractepoch
        commands = {};
        for f = 1:length(ff)
            eval( [ 'eventvals = {events(indtle).' ff{f} '};' ]);
            %if isnumeric(eventvals{1})
                %eventvals = cellfun(@num2str, eventvals, 'uniformoutput', false);
            %    eventvals = [ eventvals{:} ];
            %end
            commands = { commands{:} ff{f} eventvals };
        end
        trialinfo = struct(commands{:});
        STUDY.datasetinfo(index).trialinfo = trialinfo;
    end
    alltrialinfo{index} = trialinfo;

%    % same as above but 10 times slower
%     for e = 1:length(ALLEEG(index).event)
%         if eventlat(e) < 0.0005 % time locking event only
%             epoch = events(e).epoch;
%             for f = 1:length(ff)
%                 fieldval  = getfield(events, {e}, ff{f});
%                 trialinfo = setfield(trialinfo, {epoch}, ff{f}, fieldval);
%             end
%         end
%     end
end

% check HED field with unicode characters
% ---------------------------------------
if isfield(STUDY.datasetinfo, 'trialinfo')
    if isfield(STUDY.datasetinfo(1).trialinfo(1), 'HED')
        changeUnicode = false;
        for iDat = 1:length(STUDY.datasetinfo)
            for iTrial = 1:length(STUDY.datasetinfo(iDat).trialinfo)
                STUDY.datasetinfo(iDat).trialinfo(iTrial).HED(STUDY.datasetinfo(iDat).trialinfo(iTrial).HED > 255) = 32;
                changeUnicode = true;
            end
        end
        if changeUnicode
            disp('Unicode character detected in HED event field, removing them to improve stability');
        end
    end
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [tfr,t,f] = tfrmhs(x,t,N,g,h,trace);
%TFRMHS	Margenau-Hill-Spectrogram time-frequency distribution.
%	[TFR,T,F]=TFRMHS(X,T,N,G,H,TRACE) computes the Margenau-Hill-Spectrogram
%	distribution of a discrete-time signal X, or the cross
%	Margenau-Hill-Spectrogram representation between two signals.
%
%	X     : Signal if auto-MHS, or [X1,X2] if cross-MHS.
%	T     : time instant(s)          (default : 1:length(X)).
%	N     : number of frequency bins (default : length(X)).
%	G,H   : analysis windows, normalized so that the representation
%               preserves the signal energy.
%	                (default : Hamming(N/10) and Hamming(N/4)).
%	TRACE : if nonzero, the progression of the algorithm is shown
%                                        (default : 0).
%	TFR   : time-frequency representation. When called without
%               output arguments, TFRMHS runs TFRQVIEW.
%	F     : vector of normalized frequencies.
%
%	Example:
%	 sig=fmlin(128,0.1,0.4); g=tftb_window(21,'Kaiser');
%	 h=tftb_window(63,'Kaiser'); tfrmhs(sig,1:128,64,g,h,1);
%
%	See also all the time-frequency representations listed in
%	 the file CONTENTS (TFR*)

%	F. Auger, May-August 1994, July 1995.
%	Copyright (c) 1996 by CNRS (France).
%
%  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 2 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, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

if (nargin == 0),
 error('At least 1 parameter required');
end;
[xrow,xcol] = size(x);
if (xcol==0)|(xcol>2),
 error('X must have one or two columns');
end

if (nargin <= 2),
 N=xrow;
elseif (N<0),
 error('N must be greater than zero');
elseif (2^nextpow2(N)~=N),
 fprintf('For a faster computation, N should be a power of two\n');
end;

hlength=floor(N/4); hlength=hlength+1-rem(hlength,2);
glength=floor(N/10);glength=glength+1-rem(glength,2);

if (nargin == 1),
 t=1:xrow; g = tftb_window(glength); h = tftb_window(hlength); trace = 0;
elseif (nargin == 2)|(nargin == 3),
 g = tftb_window(glength); h = tftb_window(hlength); trace = 0;
elseif (nargin == 4),
 h = tftb_window(hlength); trace = 0;
elseif (nargin == 5),
 trace = 0;
end;

[trow,tcol] = size(t);
if (trow~=1),
 error('T must only have one row');
end;

[grow,gcol]=size(g); Lg=(grow-1)/2;
if (gcol~=1)|(rem(grow,2)==0),
 error('G must be a smoothing window with odd length');
end;

[hrow,hcol]=size(h); Lh=(hrow-1)/2; h=h/h(Lh+1);
if (hcol~=1)|(rem(hrow,2)==0),
  error('H must be a smoothing window with odd length');
end;

Lgh=min(Lg,Lh); points=-Lgh:Lgh;
Kgh=sum(h(Lh+1+points).*conj(g(Lg+1+points))); h=h/Kgh;

tfr= zeros (N,tcol); tfr2= zeros(N,tcol);
if trace, disp('Pseudo Margenau-Hill distribution'); end;
for icol=1:tcol,
 ti= t(icol);
 if trace, disprog(icol,tcol,10); end;
 tau=-min([round(N/2)-1,Lg,ti-1]):min([round(N/2)-1,Lg,xrow-ti]);
 indices= rem(N+tau,N)+1;
 tfr(indices,icol)=x(ti+tau,1).*conj(g(Lg+1+tau));
 tau=-min([round(N/2)-1,Lh,ti-1]):min([round(N/2)-1,Lh,xrow-ti]);
 indices= rem(N+tau,N)+1;
 tfr2(indices,icol)=x(ti+tau,xcol).*conj(h(Lh+1+tau));
end;
if trace, fprintf('\n'); end;
tfr=real(fft(tfr).*conj(fft(tfr2)));

if (nargout==0),
 tfrqview(tfr,x,t,'tfrmhs',g,h);
elseif (nargout==3),
 if rem(N,2)==0,
  f=[0:N/2-1 -N/2:-1]'/N;
 else
  f=[0:(N-1)/2 -(N-1)/2:-1]'/N;
 end;
end;

author: HeLiangHIT


text:
%% ZeroVolShape
% Subclass of <Shape.html |Shape|> representing a shape without volume.

%%% Description
% |ZeroVolShape| is the superclass of shapes without volume such as a plane,
% point, line.  |ZeroVolShape| is designed to help visualization of shapes
% without volume, because such shapes are hard to visualize.

%%% See Also
% <Plane.html |Plane|>, <Rectangle.html |Rectangle|>, <Line.html |Line|>,
% <Point.html |Point|>, <maxwell_run.html |maxwell_run|>

classdef ZeroVolShape < Shape

% 	methods (Abstract)
% 		draw2d(this, axes_handle, normal_axis, intercept)
% 		draw3d(this, axes_handle)
% 	end

	methods
		function this = ZeroVolShape(lprim_cell, lsf, dl_max)
			function level = lsf_zv(x, y, z, force_draw)
				if nargin < 4  % no force_draw
					force_draw = false;
				end

				level = lsf(x, y, z);
				level_copy = level(:);
				if force_draw && all(level_copy <= 0)
					% If the level set function is negative, nothing is drawn,
					% so shift the function upward.
					max_level = max(level_copy);
					max_level = max(level_copy(level_copy<max_level));  % second largest
% 					max_level = max(level_copy(level_copy<max_level));  % third largest
					level = level - max_level;
				end
			end

			if nargin < 3  % no dl_max
				super_args = {lprim_cell, @lsf_zv};
			else
				super_args = {lprim_cell, @lsf_zv, dl_max};
			end
			this = this@Shape(super_args{:});
		end
	end
end

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


text:
function [NC,NE] = snap_points_to_close_edges(C,E,epsilon)
  % SNAP_POINTS_TO_CLOSE_EDGES snap points to edges closer than a given epsilon
  % breaking that edge into two edges and updating given the edge list.
  %
  % [NC,NE] = snap_points_to_close_edges(C,E)
  % [NC,NE] = snap_points_to_close_edges(C,E,epsilon)
  %
  % Inputs:
  %   C  #C by dim list of point positions
  %   E  #E by 2 list of edges
  %   epsilon  minium distance allowed after collapses are complete, default is to
  %     use fraction of maximum edge length
  % Outputs:
  %   NC  #NC by dim list of new point positions
  %   NE  #NE by 2 list of new edges
  %
  % Example:
  %   %% Break edges
  %   % point position list
  %   C = [1,0;-1,0;0,1;eps,eps;0,-1;1,1;-1,-1];
  %   % edge list
  %   E = [1 2; 3 4; 4 5;6 7];
  %   % plot original
  %   subplot(1,2,1);
  %   plot([C(E(:,1),1) C(E(:,2),1)]',[C(E(:,1),2) C(E(:,2),2)]', ...
  %     '-','LineWidth',1);
  %   % break edges at close points
  %   [NC,NE] = snap_points_to_close_edges(C,E,0.2);
  %   % plot result
  %   subplot(1,2,2);
  %   plot([NC(NE(:,1),1) NC(NE(:,2),1)]',[NC(NE(:,1),2) NC(NE(:,2),2)]', ...
  %     '-','LineWidth',1);

  if ~exist('epsilon','var') || isempty(epsilon)
    if size(E,2) == 2
      EE = E;
    else
      EE = edges(F);
    end
    % maximum edge length
    maxD = max(sqrt(sum((C(EE(:,1),:) - C(EE(:,2),:)).^2,2)));
    epsilon = maxD/100;
  end

  % avoid sqrts
  sqr_eps = epsilon.^2;

  % make room for outputs
  NC = C;
  NE = E;

  % set previous count to phony value to enter while loop
  prev_count = size(NE,1)+1;

  % Continue edge break iterations until we're no longer breaking anything
  while prev_count ~= size(NE,1)
    prev_count = size(NE,1);
    % compute projection of each point to each line segment
    [T,sqrD] = project_to_lines(NC,NC(NE(:,1),:),NC(NE(:,2),:));
    % each vertex seen by each edge
    NCNE = repmat(NC,[1 1 size(NE,1)]);
    % edge start positions
    S = NC(NE(:,1),:);
    % edge destination positions
    D = NC(NE(:,2),:);
    % distance of each point to each edge start
    sqrDS = ...
      squeeze(sum((NCNE - permute(repmat(S,[1 1 size(NC,1)]),[3 2 1])).^2,2));
    % distance of each point to each edge dest
    sqrDD = ...
      squeeze(sum((NCNE - permute(repmat(D,[1 1 size(NC,1)]),[3 2 1])).^2,2));
    % replace distances to edges when point is closest to start or dest endpoints
    % respectively
    sqrD(T<0) = sqrDS(T<0);
    sqrD(T>1) = sqrDD(T>1);
    % inf-out self-distances
    sqrD( ...
      sub2ind(size(sqrD),[NE(:,1);NE(:,2)],[1:size(NE,1) 1:size(NE,1)]')) = inf;
    % for each edge find the closest point
    [minD,break_q] = min(sqrD);
    % mask telling whether closest point for each edge is close enough
    close = minD<sqr_eps;
    % don't let edges consider breaking at far points
    break_q(~close) = -1;
    % each edge claims a closest point, find edges which claim their closest
    % point "first"
    [~,first] = unique(break_q,'first');
    % default is to not break an edge
    break_e = false(1,size(NE,1));
    % only break first edge to claim closest point
    break_e(first) = true;
    % only break edges at close points
    break_e = break_e & close;
    % new edges are: old edges, first parts of broken edges, second parts of
    % broken edges
    NE = [ ...
      NE(~break_e,:) ; ...
      NE(break_e,1) break_q(break_e)'; ...
      break_q(break_e)' NE(break_e,2)];
  end
end

author: alecjacobson

repo: gptoolbox


text:
%% Blob Detection
% This program demonstrates how to use BLOB to detect and filter region.
%
% Sources:
%
% * <https://github.com/opencv/opencv/blob/3.1.0/samples/cpp/detect_blob.cpp>
%

%% Image
img = cv.imread(fullfile(mexopencv.root(),'test','detect_blob.png'), 'Color',true);
img = rot90(img);

%% Detector Parameters
% this is default parameters for |cv.SimpleBlobDetector|
params = struct();
params.ThresholdStep = 10;
params.MinThreshold = 10;
params.MaxThreshold = 220;
params.MinRepeatability = 2;
params.MinDistBetweenBlobs = 10;
params.FilterByColor = false;
params.BlobColor = 0;
params.FilterByArea = false;
params.MinArea = 25;
params.MaxArea = 5000;
params.FilterByCircularity = false;
params.MinCircularity = 0.9;
params.MaxCircularity = 1e37;
params.FilterByInertia = false;
params.MinInertiaRatio = 0.1;
params.MaxInertiaRatio = 1e37;
params.FilterByConvexity = false;
params.MinConvexity = 0.95;
params.MaxConvexity = 1e37;

%%
% we are going to detect blobs with 6 different param configurations
p = repmat(params, 6, 1);

%%
% 1st: we want all
p(1).FilterByArea = true;
p(1).MinArea = 1;
p(1).MaxArea = size(img,1) * size(img,2);

%%
% 2nd: we want area between 500 and 2900 pixels
p(2).FilterByArea = true;
p(2).MinArea = 500;
p(2).MaxArea = 2900;

%%
% 3rd: we want only circular object
p(3).FilterByCircularity = true;

%%
% 4th: we want ratio inertia
p(4).FilterByInertia = true;
p(4).MinInertiaRatio = 0;
p(4).MaxInertiaRatio = 0.2;

%%
% 5th: we want convexity
p(5).FilterByConvexity = true;
p(5).MinConvexity = 0;
p(5).MaxConvexity = 0.9;

%%
% 6th: we want blob with gravity center color equal to 0 (dark blobs)
p(6).FilterByColor = true;
p(6).BlobColor = 0;

%%
% helper function to convert params struct to pairs of name/value options
getopts = @(s) reshape([fieldnames(s), struct2cell(s)]', 1, []);

%% Detect
% blob detectors loop
for i=1:numel(p)
    % create detector using specified options
    opts = getopts(p(i));
    sbd = cv.SimpleBlobDetector(opts{:});

    % detect keypoint
    kpts = sbd.detect(img);

    % draw results
    out = cv.drawKeypoints(img, kpts);
    for k=1:numel(kpts)
        out = cv.circle(out, kpts(k).pt, round(kpts(k).size), ...
            'Color',randi([0 255],[1 3]));
    end

    % show output
    figure('Name',sprintf('BLOB%d',i))
    imshow(out)

    % create a title to describe current parameters
    str = {};
    if p(i).FilterByArea
        str{end+1} = sprintf('Area range [%g to %g]', ...
            p(i).MinArea, p(i).MaxArea);
    end
    if p(i).FilterByCircularity
        str{end+1} = sprintf('Circularity range [%g to %g]', ...
            p(i).MinCircularity, p(i).MaxCircularity);
    end
    if p(i).FilterByColor
        str{end+1} = sprintf('Blob color %g', p(i).BlobColor);
    end
    if p(i).FilterByConvexity
        str{end+1} = sprintf('Convexity range [%g to %g]', ...
            p(i).MinConvexity, p(i).MaxConvexity);
    end
    if p(i).FilterByInertia
        str{end+1} = sprintf('Inertia ratio range [%g to %g]', ...
            p(i).MinInertiaRatio, p(i).MaxInertiaRatio);
    end
    title(str)
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
clc;
clearvars;
rng default;
A = magic(6);

[n1, n2] = size(A);

sz = n1 * n2;

factors = [0.05 0.1 0.2 0.3];
A
for f=factors
    m = round(f * sz);
    fprintf('Trying for factor=%f, m=%d\n', f, m);
    Omega = spx.stats.rand_subset(sz, m);
    [U, S, V] = svd(A);
    y = spx.fast.partial_svd_compose(U, diag(S), V, Omega);
    spx.io.print.vector(y);
    A2 = spx.sparse.join_data_indices(y, Omega, n1, n2);
    full(A2)
end
A

author: indigits

repo: sparse-plex


text:
function [ xstar, seed ] = rk2_tv_step ( x, t, h, q, fv, gv, seed )

%*****************************************************************************80
%
%% RK2_TV_STEP takes one step of a stochastic Runge Kutta scheme.
%
%  Discussion:
%
%    The Runge-Kutta scheme is second-order, and suitable for time-varying
%    systems.
%
%    d/dx X(t,xsi) = F ( X(t,xsi), t ) + G ( X(t,xsi), t ) * w(t,xsi)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 June 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Jeremy Kasdin,
%    Runge-Kutta algorithm for the numerical integration of
%    stochastic differential equations,
%    Journal of Guidance, Control, and Dynamics,
%    Volume 18, Number 1, January-February 1995, pages 114-120.
%
%    Jeremy Kasdin,
%    Discrete Simulation of Colored Noise and Stochastic Processes
%    and 1/f^a Power Law Noise Generation,
%    Proceedings of the IEEE,
%    Volume 83, Number 5, 1995, pages 802-827.
%
%  Parameters:
%
%    Input, real X, the value at the current time.
%
%    Input, real T, the current time.
%
%    Input, real H, the time step.
%
%    Input, real Q, the spectral density of the input white noise.
%
%    Input, external real FV, the name of the deterministic
%    right hand side function.
%
%    Input, external real GV, the name of the stochastic
%    right hand side function.
%
%    Input/output, integer SEED, a seed for the random
%    number generator.
%
%    Output, real XSTAR, the value at time T+H.
%
  a21 = 1.0;
  a31 = 0.5;
  a32 = 0.5;

  q1 = 2.0;
  q2 = 2.0;

  t1 = t;
  x1 = x;
  [ n1, seed ] = r8_normal_01 ( seed );
  w1 = n1 * sqrt ( q1 * q / h );
  k1 = h * fv ( t1, x1 ) + h * gv ( t1, x1 ) * w1;

  t2 = t1 + a21 * h;
  x2 = x1 + a21 * k1;
  [ n2, seed ] = r8_normal_01 ( seed );
  w2 = n2 * sqrt ( q2 * q / h );
  k2 = h * fv ( t2, x2 ) + h * gv ( t2, x2 ) * w2;

  xstar = x1 + a31 * k1 + a32 * k2;

  return
end


text:
%% PLANEWAVEMAXWELL1 plane wave solutions to Maxwell equations in a cube.
%
% Test Maxwell function can solve problems with complex solution and
% indefinite case.
%
% See also  planewaveMaxwell, planewaveMaxwell2
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

close all;

%% Generate an initial mesh
[node,elem,HB] = cubemesh([-1,1,-1,1,-1,1],1);
bdFlag = setboundary3(node,elem,'Neumann');

%% Parameters
maxIt = 4;
N = zeros(maxIt,1);
h = zeros(maxIt,1);
energyErr = zeros(maxIt,1);
L2Err = zeros(maxIt,1);
energyErrImag = zeros(maxIt,1);
L2ErrImag = zeros(maxIt,1);

%% Get the data of the pde
global d P omega
d = [1, pi/2, pi/2];  % in sphereical coordinate
P = [1, 0, 0];
omega = 1;
r = d(1); theta = d(2); phi = d(3);
d = [r*sin(phi)*cos(theta), r*sin(phi)*sin(theta), r*cos(phi)];
% pde = planewavedataC; % plane wave with complex coefficients
pde = planewavedata; % plane wave with real coefficients and complex solution

%% Finite Element Method
for k = 1:maxIt
    % refine grid
    [node,elem,bdFlag] = uniformrefine3(node,elem,bdFlag);
    % solve the equation
%     [u,edge,A,M] = Maxwell(node,elem,HB,pde,bdFlag);
    [u,edge,eqn] =  Maxwell1(node,elem,bdFlag,pde);
    % compute the error
    uI = edgeinterpolate1(pde.exactu,node,edge);
    L2Err(k) = sqrt(abs(real(u-uI)'*eqn.M*real(u-uI)));
    energyErr(k) = sqrt(abs(real(u-uI)'*eqn.A*real(u-uI)) + L2Err(k)^2);
    L2ErrImag(k) = sqrt(abs(imag(u-uI)'*eqn.M*imag(u-uI)));
    energyErrImag(k) = sqrt(abs(imag(u-uI)'*eqn.A*imag(u-uI)) + L2ErrImag(k)^2);
%     energyErr(k) = getHcurlerror3ND1(node,elem,pde.curlu,u);
%     L2Err(k) = getL2error3ND1(node,elem,pde.exactu,u);
    N(k) = length(u);
    h(k) = 1./(size(node,1)^(1/3)-1);
end

%% Plot convergence rates
figure(1); clf;
showrateh2(h,energyErr,1,'r-+','|| Re(u-u_h)||_A',...
           h,L2Err,1,'b-+','|| Re(u-u_h)||');
figure(2); clf;
showrateh2(h,energyErrImag,1,'r-+','|| Img(u-u_h)||_A',...
           h,L2ErrImag,1,'b-+','|| Img(u-u_h)||');

repo: ifem


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% THIS IS A GITLAB CI SCRIPT!! DO NOT RUN IT MANUALLY !!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% store call path
call_path = pwd;
disp('Path after calling: ')
disp(call_path)
% deactivate warnings to make logs more readable
warning ('off','all');
%% run script
try
    %% Load project
    disp('******************************************');
    disp('Load simulink project');
    disp('******************************************');
    drawnow;
    proj=simulinkproject(ProjectName);
    % make sure scripts are run from project root
    cd(projectRoot);
    %% Run Testcases
    disp('Path before starting the testcases: ')
    disp(pwd)
    % R2017b does not support output detail parameter
    if(strcmp(version('-release'), '2018b'))
        results = runtests(TestCase, 'OutputDetail', 4);
    else
        results = runtests(TestCase);
    end
    %% Format results
    % create table and remove details column (no information in string
    % format available)
    results_tab = table(results);
    results_tab.Details = [];
    % the partial removal of HTML tags leads to the string being not
    % interpreted as formatted
    tmp = erase(evalc('disp(results_tab)'), '<strong>');
    % remove the rest of the HTML tags to make the output readable in the
    % command line
    tmp2 = erase(tmp, '</strong>');
    disp(tmp2);

    %% Handle exit
    % this if clause has to be evaluated first and exit code has to be stored since variables are
    % removed during close(proj)
    if(any([results.Failed]))
        % close project actively to prevent matlab asking to save changes
        close(proj);
        % close everything which is there forcefully, otherwise unsaved changes
        % dialogs can block CI
        bdclose('all');
        Simulink.data.dictionary.closeAll('-discard');
        exit_code = 1;
    else
        % close project actively to prevent matlab asking to save changes
        close(proj);
        % close everything which is there forcefully, otherwise unsaved changes
        % dialogs can block CI
        bdclose('all');
        Simulink.data.dictionary.closeAll('-discard');
        exit_code = 0;
    end
 catch e
    disp('******************************************');
    disp('!ERROR: Test failed');
    disp('******************************************');
    disp('******************************************');
    disp('Exception: ');
    disp(getReport(e))
    disp('******************************************');
    % close project actively to prevent matlab asking to save changes
    close(proj);
    % close everything which is there forcefully, otherwise unsaved changes
    % dialogs can block CI
    bdclose('all');
    Simulink.data.dictionary.closeAll('-discard');
    exit_code = 1;
end

quit(exit_code)

author: TUMFTM


text:
function v = s2v(s,a)

% s2v - structure array to vector
%
%   v = s2v(s,a);
%
%   v(i) = getfield(s(i), a);
%
%   Copyright (c) 2010 Gabriel Peyre

v = zeros(length(s),1);
for i=1:length(s)
    v(i) = getfield(s(i), a);
end

author: gpeyre

repo: matlab-toolboxes


text:
function [propValue, srcObj] = get(obj, propName)
% Accessor for reading BstPanel attributes.

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2008

% Check if the field propName exists in the object
if (ismember(propName, fields(obj)))
    % If it exists, try to return its value
    if (length(obj) == 1)
        propValue = obj.(propName);
        srcObj    = obj;
    else
        propValue = {obj.(propName)};
        srcObj    = obj;
    end

else
    %iPanel = find(cellfun(@(f)isfield(f, propName), obj.sControls), 1);
    if isfield(obj.sControls, propName)
        propValue = obj.sControls.(propName);
        srcObj    = obj;
    else
        propValue = [];
        srcObj    = [];
    end
end

end

author: brainstorm-tools


text:
function [vRate,vTime] = plot_LogOmori(mCatalog,mMain)
    % function [vRate,vBin] = plot_LogOmori(mCatalog,fTMain)
    % -------------------------------------------------------
    % Create daily rate plot of aftershock sequence in loglog space
    %
    % Outgoing:
    % vRate : Events per day
    % vTime : Time bin

    vAfdate = datenum(floor(mCatalog(:,3)),mCatalog(:,4),mCatalog(:,5),mCatalog(:,8),mCatalog(:,9),mCatalog(:,10));
    fTmain = datenum(floor(mMain(:,3)),mMain(:,4),mMain(:,5),mMain(:,8),mMain(:,9),mMain(:,10));
    vTime = vAfdate-fTmain;

    [vRate,vTime] = hist(vTime,0:1:ceil(max(vTime)));

    figure
    loglog(vTime,vRate,'s','Linestyle','none','Markersize',8,'Color',[0 0 0]);
    xlabel('Time [years]','FontSize',14,'Fontweight','bold');
    ylabel('Number of events [per day]','FontSize',14,'Fontweight','bold');

author: CelsoReyes


text:
function write_rgbe_line(buffer_line, fid)
%
%       write_rgbe_line(buffer_line, file)
%
%       This function write an image line using RLE and RGBE encoding
%
%        Input:
%           -buffer_line: a color line of an image
%           -file: the outputting file
%
%     Copyright (C) 2014  Francesco Banterle
%
%     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/>.
%

n = length(buffer_line);

cur_pointer = 1;

while(cur_pointer <= n )
    run_length = 0;
    run_length_old = 0;

    run_start = cur_pointer;

    %we need to find a long run; length>3
    while((run_length<4) && (run_start <= n ))
        run_start = run_start + run_length;
        run_length_old = run_length;

        i_end = min([run_start + 127 - 1, n]);

        run_length = 1;

        %finding a run
        for i=(run_start + 1):i_end
            if(buffer_line(run_start) == buffer_line(i))
                run_length = run_length + 1;
            else
                break
            end
        end
    end

    %do we have a short run <4 before a long one?
    if((run_length_old > 1) && (run_length_old == (run_start - cur_pointer)))
        length_to_write = run_length_old + 128;
        value_to_write = buffer_line(cur_pointer);
        fwrite(fid, length_to_write, 'uint8');
        fwrite(fid, value_to_write, 'uint8');

        cur_pointer = run_start;
    end

    %writing non-runs
    while(cur_pointer < run_start)
        non_run_length = run_start - cur_pointer;

        if(non_run_length > 128)
            non_run_length = 128;
        end

        fwrite(fid, non_run_length, 'uint8');
        fwrite(fid, buffer_line(cur_pointer : (cur_pointer + non_run_length - 1)), 'uint8');

        cur_pointer = cur_pointer + non_run_length;
    end

    %writing the found long run
    if(run_length > 3)
        length_to_write = run_length + 128;
        value_to_write = buffer_line(run_start);
        fwrite(fid, length_to_write, 'uint8');
        fwrite(fid, value_to_write, 'uint8');

        cur_pointer = cur_pointer + run_length;
    end

end

end

author: banterle


text:
% QUATERNION CONJUGATE
function [qCon] = qConjugate(q)
% This function defines the conjugate of the given quaternion.
% Associated block:
% "Quaternion Conjugate"
qCon = zeros(4,1);
qCon(2) = -q(2);
qCon(3) = -q(3);
qCon(4) = -q(4);
end

repo: OpenMAS


text:
function [err,data] = afni_ni_checkfordata()
   global comm;

   data = [];
   err = 1;
   if (comm.sockcon < 0),
      comm.sockcon = pnet('tcpsocket',sprintf('%d', comm.socknum));
      if (comm.sockcon < 0),
         fprintf(2,'Could not open socket %d\n', comm.socknum);
         return;
      end
      if (comm.dbg ),
         fprintf(1,  'have new socket connection %d\n', comm.sockcon);
      end
   else
      if (comm.dbg > 1) fprintf(1,'have socket %d\n', comm.sockcon); end
   end
   %have socket, try connection
   if (comm.con < 0),
      comm.con = pnet(comm.sockcon,'tcplisten','noblock');  % Don't wait forever
      err = 0;
      if (comm.dbg > 1) fprintf(1,'No connection yet \n'); end
      return;  %nothing yet
   else
      if (comm.dbg > 1) fprintf(1,'have connection handle %d\n', comm.con); end
   end
   status = pnet(comm.con,'status');
   if (status == 0),
      % bad status, perhaps disconnected or first time, reopen socket
      fprintf(2,'Stream closed, attempt to reopen on second pass.\n');
      comm.con = -1;
      err = 0;
      return;
   else
      if (comm.dbg > 1) fprintf(2,'status %d\n', status); end
   end

   %have connection, get data
   pnet(comm.con,'setreadtimeout',0.25);
      % Don't wait forever, I had 3 seconds but that was too slow
      % 0.25 worked just fine on localhost 127.0.0.1

   %check for stuff
   b = pnet(comm.con,'read');
   while (length(b)) data=[data b]; b = pnet(comm.con,'read');  end

   err = 0;
   return

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


function res = fact1(x)
if x > 0
   res = x * fact1(x-1);
else
   res = 1;
end


text:
function newFaces = minConvexHull(points, varargin)
%MINCONVEXHULL Return the unique minimal convex hull of a set of 3D points.
%
%   FACES = minConvexHull(PTS)
%   NODES is a set of 3D points  (as a Nx3 array). The function computes
%   the convex hull, and merge contiguous coplanar faces. The result is a
%   set of polygonal faces, such that there are no coplanar faces.
%   FACES is a cell array, each cell containing the vector of indices of
%   nodes given in NODES for the corresponding face.
%
%   FACES = minConvexHull(PTS, PRECISION)
%   Adjust the threshold for deciding if two faces are coplanar or
%   parallel. Default value is 1e-14.
%
%   Example
%     % extract square faces from a cube
%     [n, e, f] = createCube;
%     f2 = minConvexHull(n);
%     drawMesh(n, f2);
%
%     % Subdivides and smooths a mesh rpresenting a cube
%     [n, e, f] = createCube;
%     [n2, f2] = subdivideMesh(n, triangulateFaces(f), 4);
%     [n3, f3] = smoothMesh(n2, f2);
%     figure; drawMesh(n3, f3);
%     axis equal; view(3);
%     % merge coplanar faces, making apparent the faces of the original cube
%     f4 = minConvexHull(n3);
%     figure; drawMesh(n3, f4);
%     axis equal; view(3);
%
%
%   See also
%   meshes3d, mergeCoplanarFaces, drawMesh, convhull, convhulln
%

% ------
% Author: David Legland
% E-mail: david.legland@inra.fr
% Created: 2006-07-05
% Copyright 2006-2022 INRA - CEPIA Nantes - MIAJ (Jouy-en-Josas)

% set up precision
acc = 1e-14;
if ~isempty(varargin)
    acc = varargin{1};
end

% triangulated convex hull. It is not uniquely defined.
faces = convhulln(points);

% compute centroid of the nodes
pointsCentroid = centroid(points);

% number of base triangular faces
N = size(faces, 1);

% compute normals of given faces
normals = planeNormal(createPlane(...
    points(faces(:,1),:), points(faces(:,2),:), points(faces(:,3),:)));

% initialize empty faces
newFaces = {};

% Processing flag for each triangle
% 1 : triangle to process, 0 : already processed
% in the beginning, every triangle face need to be processed
flag = ones(N, 1);

% iterate on each triangular face of the convex hull
for iFace = 1:N

    % check if face was already performed
    if ~flag(iFace)
        continue;
    end

    % indices of faces with same normal
    ind = find(abs(vectorNorm3d(cross(repmat(normals(iFace, :), [N 1]), normals)))<acc);
    ind = ind(ind~=iFace);

    % keep only coplanar faces (test coplanarity of points in both face)
    ind2 = iFace;
    for j = 1:length(ind)
        if isCoplanar(points([faces(iFace,:) faces(ind(j),:)], :), acc)
            ind2 = [ind2 ind(j)]; %#ok<AGROW>
        end
    end

    % compute order of the vertices in current face
    faceVertices = unique(faces(ind2, :));
    [tmp, I]  = angleSort3d(points(faceVertices, :)); %#ok<ASGLU>

    % create the new face, ensuring it is a row vector
    face = faceVertices(I);
    face = face(:)';

    % ensure face has normal pointing outwards
    outerNormal = meshFaceCentroids(points, face) - pointsCentroid;
    if dot(meshFaceNormals(points, face), outerNormal, 2) < 0
        face = face([1 end:-1:2]);
    end

    % add a new face to the list
    newFaces = [newFaces {face}]; %#ok<AGROW>

    % mark processed faces
    flag(ind2) = 0;
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function kernel = calibrate4D(obj, AtA, nCha, coil)
% retrieve calibration kernel from 4D calibration region
%
% (c) Thomas Kuestner
% ---------------------------------------------------------------------

sampling = ones([obj.kernelSize,nCha],obj.measPara.precision);

dummyK = zeros([obj.kernelSize,nCha],obj.measPara.precision);
dummyK((end+1)/2,(end+1)/2,(end+1)/2,coil) = 1;
idxY = find(dummyK);
sampling(idxY) = 0;
idxA = find(sampling);

Aty = AtA(:,idxY); Aty = Aty(idxA); % correlation values to target point, take complete neighbourhood and not just aquired ones
AtA = AtA(idxA,:); AtA =  AtA(:,idxA); % kick out the searched point

kernel = sampling*0;

lambda = norm(AtA,'fro')/size(AtA,1)*obj.calibTyk;

cnd = cond(AtA + eye(size(AtA))*lambda);
if(cnd == inf)
    rawkernel = pinv(AtA + eye(size(AtA))*lambda)*Aty; % grappa weighting values
else
    rawkernel = (AtA + eye(size(AtA))*lambda)\Aty; % grappa weighting values
end
kernel(idxA) = rawkernel;

end

author: thomaskuestner


text:
%% Grain Boundary Tutorial
% A quick guide to grain boundary analysis
%
%% Grain boundaries generation
%
% To work with grain boundaries we need some ebsd data and have to detect
% grains within the data set.

% load some example data
mtexdata twins

% detect grains
[grains,ebsd.grainId,ebsd.mis2mean] = calcGrains(ebsd('indexed'))

% smooth them
grains = grains.smooth

% visualize the grains
plot(grains,grains.meanOrientation)

%%
% Now we can extract from the grains its boundary and save it to a seperate
% variable

gB = grains.boundary

%%
% The output tells us that we have 3219 Magnesium to Magnesium boundary
% segments and 606 boundary segements where the grains are cutted by the
% scanning boundary. To restrict the grain boundaries to a specific phase
% transistion you shall do

gB_MgMg = gB('Magnesium','Magnesium')

%% Properties of grain boundaries
%
% A variable of type grain boundary contains the following properties
%
% * misorientation
% * direction
% * segLength
%
% These can be used to colorize the grain boundaries. By the following
% command we plot the grain boundaries colorized by the misorientation
% angle

plot(gB_MgMg,gB_MgMg.misorientation.angle./degree,'linewidth',2)
mtexColorbar

%%

hold on
plot(gB('notIndexed'),'lineColor','blue','linewith',5)
hold off

%%

grains.innerBoundary

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function [MATreordered,MATindices,MATcost] = reorderMAT(MAT,H,cost)
%REORDERMAT         Reorder matrix for visualization
%
%   [MATreordered,MATindices,MATcost] = reorderMAT(MAT,H,cost);
%
%   This function reorders the connectivity matrix in order to place more
%   edges closer to the diagonal. This often helps in displaying community
%   structure, clusters, etc.
%
%   Inputs:     MAT,            connection matrix
%               H,              number of reordering attempts
%               cost,           'line' or 'circ', for shape of lattice
%                               (linear or ring lattice)
%
%               MATreordered    reordered connection matrix
%               MATindices      reordered indices
%               MATcost         cost of reordered matrix
%
%
%   Olaf Sporns, Indiana University

N = length(MAT);
diagMAT = diag(diag(MAT));
MAT = MAT-diagMAT;

% generate cost function
if strcmp(cost,'line');
    profil = fliplr(normpdf(1:N,0,N/2));
end;
if strcmp(cost,'circ');
    profil = fliplr(normpdf(1:N,N/2,N/4));
end;
COST = toeplitz(profil,profil);

% initialize lowCOST
lowMATcost = sum(sum(COST.*MAT));

% keep track of starting configuration
MATstart = MAT;
starta = 1:N;

% reorder
for h=1:H
    a = 1:N;
    % choose two positions at random and flip them
    r = randperm(N);
    a(r(1)) = r(2);
    a(r(2)) = r(1);
    MATcostnew = sum(sum(MAT(a,a).*COST));
    if (MATcostnew < lowMATcost)
        MAT = MAT(a,a);
        r2 = starta(r(2));
        r1 = starta(r(1));
        starta(r(1)) = r2;
        starta(r(2)) = r1;
        lowMATcost = MATcostnew;
    end;
end;	% h

MATreordered = MATstart(starta,starta) + diagMAT(starta,starta);
MATindices = starta;
MATcost = lowMATcost;


text:
function obj = meas(varargin)

% Constructor for meas class object.
% You must always pass one argument if you want to create a new object.

if nargin==0 % Used when objects are loaded from disk
  obj = init_fields;
  obj = class(obj, 'meas');
  return;
end

firstArg = varargin{1};
if isa(firstArg, 'meas') %  used when objects are passed as arguments
  obj = firstArg;
  return;
end

% We must always construct the fields in the same order,
% whether the object is new or loaded from disk.
% Hence we call init_fields to do this.
obj = init_fields;

% attach class name tag, so we can call member functions to
% do any initial setup
obj = class(obj, 'meas');

% Now the real initialization begins
obj.value = varargin{1};
obj.error = varargin{2};

return;
%%%%%%%%%

function obj = init_fields()
% Initialize all fields to dummy values
obj.value = [];
obj.error = [];

author: Sable

repo: mcbench-benchmarks


text:
function [h, p, ci, stats] = ttest2(D, varname, wh_keep1, wh_keep2, varargin)
% Two sample ttest for two samples of one subject-level variable
%
% :Usage:
% ::
%
%    ttest2(D, varname, wh_keep1, wh_keep2, [optional inputs])
%
% ..
%     Author and copyright information:
%
%     Copyright (C) 2013 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:
%
%   **D:**
%        a canlab_dataset object
%
%   **varname:**
%        the name of a valid variable to get from dataset
%
%   **wh_keep1:**
%        subjects forming first sample
%
%   **wh_keep2:**
%        subjects forming second sample
%
% :Optional Inputs:
%
%   **noverbose:**
%         will suppress print out of results and bargraph
%
%   **varargin:**
%         other variables passed directly to MATLAB's ttest2
%
% :Outputs:
%
%   same as MATLAB's ttest2 output
%

if any(wh_keep1 & wh_keep2), warning('YOUR SAMPLES ARE OVERLAPPING!!'); end

verbose=1;
if any(strcmp('noverbose', varargin))
    verbose=0;
    varargin(find(strcmp('noverbose', varargin))) = [];
end

x1 = get_var(D, varname, wh_keep1);
x2 = get_var(D, varname, wh_keep2);

[h, p, ci, stats] = ttest2(x1, x2, varargin{:});

if verbose, ttest2_printout(x1,x2, 1); end

end

author: canlab


text:
function output = callsnopt(model)

% Build trees etc. This common format is used for fmincon, snopt and ipopt
% interfaces
model = yalmip2nonlinearsolver(model);

if ~model.derivative_available
    disp('Derivate-free call to snopt not yet implemented')
    error('Derivate-free call to snopt not yet implemented')
end

% Figure out which variables are artificially introduced to normalize
% arguments in callback operators to simplify chain rules etc. We can do
% our function evaluations and gradient computations in our lifted world,
% but only expose the model in the original variables to the nonlinear
% solver.
% model = compressLifted(model);

if model.options.savedebug
    save snoptdebug model
end

showprogress('Calling SNOPT',model.options.showprogress);

x0 = model.x0;
solveopt = 1;
xlow = model.lb;
xupp = model.ub;
xmul   = zeros(length(xupp),1);
xstate = zeros(length(xupp),1);

Fupp = [ inf;
         repmat(0,length(model.bnonlinineq),1);
         repmat(0,nnz(model.K.q),1);
         repmat(0,length(model.bnonlineq),1);
         repmat(0,length(model.b),1);
         repmat(0,length(model.beq),1)];

Flow = [ -inf;
         repmat(-inf,length(model.bnonlinineq),1);
         repmat(-inf,nnz(model.K.q),1);
         repmat(0,length(model.bnonlineq),1);
         repmat(-inf,length(model.b),1);
         repmat(0,length(model.beq),1)];

Fmul   = zeros(length(Fupp),1);
Fstate = zeros(length(Fupp),1);
ObjAdd = 0;
ObjRow = 1;
A = [];
iAfun = [];
jAvar = [];
% Sparsity pattern of jacobian
if ~isempty(Fupp)
    G = jacobiansparsityfromnonlinear(model);
else
    G = [];
end
% Add a row for objective. No sparsity declared
G = [ones(1,size(G,2));G];
[iGfun,jGvar] = find(G);
model.sparsityElements = find(G);

usrf = 'snopt_callback';
snopt_callback([],model);

global latest_xevaled
global latest_x_xevaled
latest_xevaled = [];
latest_x_xevaled = [];

solvertime = tic;
if strcmpi(model.solver.version,'cmex')
    % Some old interface? Keep for safety
    if model.options.verbose == 0
        snscreen('off')
    else
        snscreen('on');
    end
    snseti('Minimize',1)
    [xout,F,xmul,Fmul,inform, xstate, Fstate, ns, ninf, sinf, mincw, miniw, minrw] = snoptcmex( solveopt, x0, xlow, xupp, xmul, xstate, Flow, Fupp, Fmul, Fstate,ObjAdd, ObjRow, A, iAfun(:), jAvar(:),iGfun(:), jGvar(:), usrf );
else
    Astruct.A = A;
    Astruct.row = iAfun;
    Astruct.col = jAvar;
    if model.options.verbose == 0
       model.options.snopt.screen = 'off';
    end
    [xout,F,inform,xmul,Fmul] = snopt(x0, xlow, xupp, xmul, xstate,Flow, Fupp, Fmul, Fstate, usrf,ObjAdd, ObjRow,A,G,model.options.snopt);
end

solvertime = toc(solvertime);

lambda = Fmul(2:end);

if ~isempty(xout) && ~isempty(model.lift);
    x = zeros(length(model.linearindicies),1);
    x(model.lift.linearIndex) = xout;
    x(model.lift.liftedIndex) = model.lift.T*xout + model.lift.d;
    x = RecoverNonlinearSolverSolution(model,x);
else
    x = RecoverNonlinearSolverSolution(model,xout);
end

problem = 0;

% Internal format for duals
D_struc = [];

% Check, currently not exhaustive...
switch inform
    case {1}
        problem = 0;
    case {11,12,13,14,40,43,91} % 1 is sent when I test
        problem = 1;
    case 21
        problem = 2;
    case {31,32}
        problem = 3;
    case {2,3,33,41,42,43,44}
        problem = 4;
    case {71,72,73,74}
        problem = 16;
    otherwise
        problem = 11;
end

% Save all data sent to solver?
if model.options.savesolverinput
    solverinput.model = model;
else
    solverinput = [];
end

% Save all data from the solver?
if model.options.savesolveroutput
    solveroutput.x = x;
    solveroutput.F = F;
    solveroutput.xmul = xmul;
    solveroutput.inform=inform;
    solveroutput.Fstate=Fstate;
    solveroutput.ns = ns;
    solveroutput.ninf = ninf;
    solveroutput.sinf = sinf;
    solveroutput.mincw = mincw;
    solveroutput.miniw = miniw;
    solveroutput.miniw = miniw;

else
    solveroutput = [];
end

% Standard interface
output = createOutputStructure(x,D_struc,[],problem,'SNOPT',solverinput,solveroutput,solvertime);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function stickModify(handle, values, connect)

% STICKMODIFY Helper code for visualisation of a stick man.
% FORMAT
% DESC updates a stick man representation in a 3-D plot.
% ARG handle : a vector of handles to the structure to be updated.
% ARG vals : the x,y,z channels to update the skeleton with.
% ARG connect : the connectivity of the skeleton.
%
% SEEALSO : stickVisualise
%
% COPYRIGHT : Neil D. Lawrence, 2005, 2006

% MOCAP

vals = reshape(values, size(values, 2)/3, 3);

indices = find(connect);
[I, J] = ind2sub(size(connect), indices);
set(handle(1), 'Xdata', vals(:, 1), 'Ydata', vals(:, 2), 'Zdata', ...
                 vals(:, 3));
%set(handle(1), 'visible', 'off')

for i = 1:length(indices)
  set(handle(i+1), 'Xdata', [vals(I(i), 1) vals(J(i), 1)], ...
            'Ydata', [vals(I(i), 2) vals(J(i), 2)], ...
            'Zdata', [vals(I(i), 3) vals(J(i), 3)]);
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [ n_data, x, fx ] = bessel_j0_values ( n_data )

%*****************************************************************************80
%
%% BESSEL_J0_VALUES returns some values of the J0 Bessel function.
%
%  Discussion:
%
%    In Mathematica, the function can be evaluated by:
%
%      BesselJ[0,x]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real X, the argument of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 21;

  fx_vec = [ ...
     -0.1775967713143383E+00, ...
     -0.3971498098638474E+00, ...
     -0.2600519549019334E+00, ...
      0.2238907791412357E+00, ...
      0.7651976865579666E+00, ...
      0.1000000000000000E+01, ...
      0.7651976865579666E+00, ...
      0.2238907791412357E+00, ...
     -0.2600519549019334E+00, ...
     -0.3971498098638474E+00, ...
     -0.1775967713143383E+00, ...
      0.1506452572509969E+00, ...
      0.3000792705195556E+00, ...
      0.1716508071375539E+00, ...
     -0.9033361118287613E-01, ...
     -0.2459357644513483E+00, ...
     -0.1711903004071961E+00, ...
      0.4768931079683354E-01, ...
      0.2069261023770678E+00, ...
      0.1710734761104587E+00, ...
     -0.1422447282678077E-01 ];

  x_vec = [ ...
     -5.0E+00, ...
     -4.0E+00, ...
     -3.0E+00, ...
     -2.0E+00, ...
     -1.0E+00, ...
      0.0E+00, ...
      1.0E+00, ...
      2.0E+00, ...
      3.0E+00, ...
      4.0E+00, ...
      5.0E+00, ...
      6.0E+00, ...
      7.0E+00, ...
      8.0E+00, ...
      9.0E+00, ...
     10.0E+00, ...
     11.0E+00, ...
     12.0E+00, ...
     13.0E+00, ...
     14.0E+00, ...
     15.0E+00 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    x = 0.0;
    fx = 0.0;
  else
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end

author: johannesgerer


text:
function [ B_cone, B_rod ] = BleachingParameters( A_cone, A_rod )
%
%       [ B_cone, B_rod ] = BleachingParameters( A_cone, A_rod )
%
%       This function computes sigmoid response
%
%       input:
%           -A_cone: adaptation for cones in cd/m^2
%           -A_rod: adaptation for rods in cd/m^2
%
%       output:
%           -B_cone: bleaching parameter for cones
%           -B_rod: bleaching parameter for rods
%
%     Copyright (C) 2011-14  Francesco Banterle
%
%     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/>.
%

B_cone = 2 * 1e6 / (2*1e6 + A_cone);
B_rod = 0.04 / (0.04 + A_rod);

end


text:
function N = R1(rad)
% Suppose right handed coordinate frames A and B have the same x axis,
% and rotate the A frame along the x axis by rad obtains the B frame, then
% a point in B, p_B, can be transformed from its coordinates in A frame by
% p_B = R1(rad) * p_A.
N = [1 0 0;0 cos(rad) sin(rad);0 -sin(rad) cos(rad)];

repo: ekfmonoslam


text:
function h=window(N,type, varargin)

% window
%
% h=window(N,type, varargin)
%
% Generate window of length N (in samples)
%
%
% Possible types are :
% 'Hamming', 'Hanning', 'Nuttall',  'Papoulis', 'Harris',
% 'Rect',    'Triang',  'Bartlett', 'BartHann', 'Blackman'
% 'Gauss',   'Parzen',  'Kaiser',   'Dolph',    'Hanna'.
% 'Nutbess', 'spline'
%
% For the gaussian window, an optionnal parameter K
% sets the value at both extremities. The default value is 0.005
%
% For the Kaiser-Bessel window, an optionnal parameter
% sets the scale. The default value is 3*pi.
%
% For the Spline windows, h=window(N,'spline',nfreq,p)
% yields a spline weighting function of order p and frequency
% bandwidth proportional to nfreq.
%
% Example : w=window(256,'Gauss',0.005);

if nargin < 1
	help(mfilename);
	return
end

if  N<=0
	error('N should be strictly positive.');
end

if nargin < 2
	type= 'Hamming';
end

type=upper(type);

switch(type)
	case {'RECTANG','RECT'}
	 	h=ones(N,1);
	case 'HAMMING'
		h=0.54 - 0.46*cos(2.0*pi*(1:N)'/(N+1));
	case {'HANNING', 'HANN'}
		h=0.50 - 0.50*cos(2.0*pi*(1:N)'/(N+1));
	case 'KAISER'
		if (nargin==3), beta=varargin{1}; else beta=3.0*pi; end;
		ind=(-(N-1)/2:(N-1)/2)' *2/N; beta=3.0*pi;
		h=bessel(0,j*beta*sqrt(1.0-ind.^2))/real(bessel(0,j*beta));
	case 'NUTTALL',
		ind=(-(N-1)/2:(N-1)/2)' *2.0*pi/N;
		h=+0.3635819 ...
		+0.4891775*cos(    ind) ...
		+0.1363995*cos(2.0*ind) ...
		+0.0106411*cos(3.0*ind) ;
	case 'BLACKMAN'
		ind=(-(N-1)/2:(N-1)/2)' *2.0*pi/N;
		h= +0.42 + 0.50*cos(ind) + 0.08*cos(2.0*ind) ;
	case 'HARRIS'
		ind=(1:N)' *2.0*pi/(N+1);
		h=+0.35875 ...
		-0.48829 *cos(    ind) ...
		+0.14128 *cos(2.0*ind) ...
		-0.01168 *cos(3.0*ind);
	case {'BARTLETT','TRIANG'}
		h=2.0*min((1:N),(N:-1:1))'/(N+1);
	case 'BARTHANN'
		h=  0.38 * (1.0-cos(2.0*pi*(1:N)/(N+1))') ...
		+ 0.48 * min((1:N),(N:-1:1))'/(N+1);
	case 'PAPOULIS'
		ind=(1:N)'*pi/(N+1); h=sin(ind);
	case 'GAUSS'
		if (nargin==3), K=varargin{1}; else K=0.005; end;
		h= exp(log(K) * linspace(-1,1,N)'.^2 );
	case 'PARZEN'
		ind=abs(-(N-1)/2:(N-1)/2)'*2/N; temp=2*(1.0-ind).^3;
		h= min(temp-(1-2.0*ind).^3,temp);
	case 'HANNA'
		if (nargin==3), L=varargin{1}; else L=1; end;
		ind=(0:N-1)';h=sin((2*ind+1)*pi/(2*N)).^(2*L);
	case {'DOLPH','DOLF'}
		if (rem(N,2)==0), oddN=1; N=2*N+1; end;
		if (nargin==3), A=10^(param/20); else A=1.0e-3; end;
		K=N-1; Z0=cosh(acosh(1.0/A)/K); x0=acos(1/Z0)/pi; x=(0:K)/N;
		indices1=find((x<x0)|(x>1-x0));
		indices2=find((x>=x0)&(x<=1-x0));
		h(indices1)= cosh(K*acosh(Z0*cos(pi*x(indices1))));
		h(indices2)= cos(K*acos(Z0*cos(pi*x(indices2))));
		h=fftshift(real(ifft(A*real(h))));h=h'/h(K/2+1);
 		if oddN, h=h(2:2:K); end;
	case 'NUTBESS',
		if (nargin==3), beta=varargin{1}; nu=0.5;
		elseif (nargin==4), beta=varargin{1}; nu=varargin{2};
		else beta=3*pi; nu=0.5;
		end;
		ind=(-(N-1)/2:(N-1)/2)' *2/N;
		h=sqrt(1-ind.^2).^nu .* ...
		real(bessel(nu,j*beta*sqrt(1.0-ind.^2)))/real(bessel(nu,j*beta));
	case 'SPLINE'
		if (nargin < 3),
			error('Three or four parameters required for spline windows');
		elseif (nargin==3)
			nfreq=param; p=pi*N*nfreq/10.0;
			else nfreq=varargin{1}; p=varargin{2};
		end
		ind=(-(N-1)/2:(N-1)/2)';
		h=sinc((0.5*nfreq/p)*ind) .^ p;
	otherwise
		error('unknown window type');
end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function pass = test_feval(pref)
% Test feval

if ( nargin == 0 )
    pref = chebfunpref;
end
tol = 1e3*pref.cheb3Prefs.chebfun3eps;

seedRNG(42);

f = chebfun3(@(x,y,z) x, [-1 2 -pi/2 pi -3 1]);
pass(1) = (abs(f(0,0,0)) < tol*vscale(f));

pass(2) = (abs(f(pi/6,pi/12,-1)-pi/6) < tol*vscale(f));

f = chebfun3(@(x,y,z) y, [-1 2 -pi/2 pi -3 1]);
pass(3) = (abs(f(0,0,0)) < tol);

pass(4) = (abs(f(pi/6,pi/12,-1)-pi/12) < tol*vscale(f));

f = chebfun3(@(x,y,z) z, [-1 2 -pi/2 pi -3 1]);
pass(5) = (abs(f(0,0,0)) < tol);

pass(6) = (abs(f(pi/6,pi/12,-1)+1) < tol*vscale(f));

% some harder tests.
f = @(x,y,z) cos(x) + sin(x.*y) + sin(z.*x);
g = chebfun3(f);
pts = 2*rand(3,1) - 1;
pass(7) = (abs(f(pts(1),pts(2),pts(3)) - g(pts(1),pts(2),pts(3)))<tol*vscale(g));

% Are we evaluating on arrays correctly?
r = rand(10,1);
s = rand(10,1);
t = rand(10,1);
[rr, ss, tt]=meshgrid(r,s,t);
pass(8) = (norm((f(r,s,t) - g(r,s,t))) < tol*vscale(g));

pass(9) = (max(max(max(abs(f(rr,ss,tt) - g(rr,ss,tt))))) < tol*vscale(g));

% Does this work off [-1,1]^2
g = chebfun3(f,[-pi/6 pi/2 -pi/12 sqrt(3) -3 1]); % strange domain.
r = 0.126986816293506; s = 0.632359246225410; t = 0.351283361405006;
% three fixed random number in domain.
pass(10) = (abs(f(r,s,t) - g(r,s,t))<tol*vscale(g));

% Are we evaluating on arrays correctly
pass(11) = (norm((f(r,s,t) - g(r,s,t)))<tol*vscale(g));

pass(12) = (max(max(max(abs(f(rr,ss,tt) - g(rr,ss,tt)))))< tol*vscale(g));

%% vector inputs
ff = @(x,y,z) sin(pi*(x+y+z));
f = chebfun3(ff);
xx = linspace(-1, 1, 100)';
yy = linspace(-1, 1, 100)';
zz = linspace(-1, 1, 100)';
F = f(xx,yy,zz);
FF = ff(xx,yy,zz);
pass(13) = norm(F - FF) < 100*tol;

% row vector inputs
ff = @(x,y,z) sin(pi*(x+y+z));
f = chebfun3(ff);
xx = linspace(-1, 1, 100);
yy = linspace(-1, 1, 100);
zz = linspace(-1, 1, 100);
F = f(xx,yy,zz);
FF = ff(xx,yy,zz);
pass(14) = norm(F - FF) < 100*tol;

% 'trig' flag + vector inputs
ff = @(x,y,z) sin(pi*(x+y+z));
f = chebfun3(ff, 'trig');
xx = linspace(-1, 1, 100)';
yy = linspace(-1, 1, 100)';
zz = linspace(-1, 1, 100)';
F = f(xx,yy,zz);
FF = ff(xx,yy,zz);
pass(15) = norm(F - FF) < 100*tol;

% random vector inputs
ff = @(x,y,z) sin(pi*(x+y+z));
f = chebfun3(ff);
xx = rand(100, 1);
yy = rand(100, 1);
zz = rand(100, 1);
F = f(xx,yy,zz);
FF = ff(xx,yy,zz);
pass(16) = norm(F - FF) < 100*tol;

% 'trig' flag + random vector inputs
ff = @(x,y,z) sin(pi*(x+y+z));
f = chebfun3(ff, 'trig');
xx = rand(100, 1);
yy = rand(100, 1);
zz = rand(100, 1);
F = f(xx,yy,zz);
FF = ff(xx,yy,zz);
pass(17) = norm(F - FF) < 100*tol;

%% Matrix input with meshgrid
% x and y are matrix inputs (meshgrid) and z is a matrix containing copies
% of a single scalar
ff = @(x,y,z) sin(pi*(x+y+z));
f = chebfun3(ff);
xx = linspace(-1, 1, 100)';
[xx, yy] = meshgrid(xx);
zz = xx(1,1)*ones(size(xx));
F = f(xx,yy,zz);
FF = ff(xx,yy,zz);
pass(18) = norm(F - FF) < 100*tol;

% y and z are matrix inputs (meshgrid) and x is just a matrix of just one
% scalar
ff = @(x,y,z) sin(pi*(x+y+z));
f = chebfun3(ff);
zz = linspace(-1, 1, 100)';
[yy, zz] = meshgrid(zz);
xx = zz(1,1)*ones(size(xx));
F = f(xx,yy,zz);
FF = ff(xx,yy,zz);
pass(19) = norm(F - FF) < 100*tol;

% x and z are matrix inputs (meshgrid) and y is a matrix of just one scalar
ff = @(x,y,z) sin(pi*(x+y+z));
f = chebfun3(ff);
zz = linspace(-1, 1, 100)';
[xx, zz] = meshgrid(zz);
yy = zz(1,1)*ones(size(xx));
F = f(xx,yy,zz);
FF = ff(xx,yy,zz);
pass(20) = norm(F - FF) < 100*tol;

%% Matrix input with ndgrid
% x and y are matrix inputs (ndgrid) and z is a matrix containing copies
% of a single scalar
ff = @(x,y,z) sin(pi*(x+y+z));
f = chebfun3(ff);
xx = linspace(-1, 1, 100)';
[xx, yy] = ndgrid(xx);
zz = xx(1,1)*ones(size(xx));
F = f(xx,yy,zz);
FF = ff(xx,yy,zz);
pass(21) = norm(F - FF) < 100*tol;

% y and z are matrix inputs (ndgrid) and x is just a matrix of just one
% scalar
ff = @(x,y,z) sin(pi*(x+y+z));
f = chebfun3(ff);
zz = linspace(-1, 1, 100)';
[yy, zz] = ndgrid(zz);
xx = zz(1,1)*ones(size(xx));
F = f(xx,yy,zz);
FF = ff(xx,yy,zz);
pass(22) = norm(F - FF) < 100*tol;

% x and z are matrix inputs (ndgrid) and y is a matrix of just one scalar
ff = @(x,y,z) sin(pi*(x+y+z));
f = chebfun3(ff);
zz = linspace(-1, 1, 100)';
[xx, zz] = ndgrid(zz);
yy = zz(1,1)*ones(size(xx));
F = f(xx,yy,zz);
FF = ff(xx,yy,zz);
pass(23) = norm(F - FF) < 100*tol;

%% Tensor inputs
% Tensor input generated by meshgrid
ff = @(x,y,z) sin(pi*(x+y+z));
f = chebfun3(ff);
[xx, yy, zz] = meshgrid(linspace(-1, 1, 100));
F = f(xx,yy,zz);
FF = ff(xx,yy,zz);
pass(24) = norm(F(:) - FF(:)) < 100*tol;

% 'trig' flag + meshgrid
ff = @(x,y,z) sin(pi*(x+y+z));
f = chebfun3(ff, 'trig');
[xx, yy, zz] = meshgrid(linspace(-1, 1, 100));
F = f(xx,yy,zz);
FF = ff(xx,yy,zz);
pass(25) = norm(F(:) - FF(:)) < 100*tol;

% % Tensor input generated by ndgrid
ff = @(x,y,z) sin(pi*(x+y+z));
f = chebfun3(ff);
[xx, yy, zz] = ndgrid(linspace(-1, 1, 100));
F = f(xx,yy,zz);
FF = ff(xx,yy,zz);
pass(26) = norm(F(:) - FF(:)) < 100*tol;

% 'trig' flag + ndgrid
ff = @(x,y,z) sin(pi*(x+y+z));
f = chebfun3(ff, 'trig');
[xx, yy, zz] = ndgrid(linspace(-1, 1, 100));
F = f(xx,yy,zz);
FF = ff(xx,yy,zz);
pass(27) = norm(F(:) - FF(:)) < 100*tol;

% random tensor inputs
ff = @(x,y,z) sin(pi*(x+y+z));
f = chebfun3(ff);
xx = rand(10, 20, 30);
yy = rand(10, 20, 30);
zz = rand(10, 20, 30);
F = f(xx,yy,zz);
FF = ff(xx,yy,zz);
pass(28) = norm(F(:) - FF(:)) < 100*tol;

% 'trig' flag + random tensor inputs
ff = @(x,y,z) sin(pi*(x+y+z));
f = chebfun3(ff, 'trig');
xx = rand(10, 20, 30);
yy = rand(10, 20, 30);
zz = rand(10, 20, 30);
F = f(xx,yy,zz);
FF = ff(xx,yy,zz);
pass(29) = norm(F(:) - FF(:)) < 100*tol;

%% Cross sections
% Fixed x
ff = @(x,y,z) sin(x+y+z) + x + y; % A function with different variable-ranks
dom = [-1 1 -4 -2 6 8];
f = chebfun3(ff, dom);
f2D = f(0.5, :, :);
fChebfun2 = chebfun2(@(y,z) ff(0.5, y, z), dom(3:6));
pass(30) = norm(fChebfun2 - f2D) < 100*tol;

% APA's bug report for a function that has x-rank 1:
ff = @(x, y, z) cos(y + z).*sin(z).*exp(x);
f = chebfun3(ff, dom);
f2D = f(0.5, :, :);
fChebfun2 = chebfun2(@(y,z) ff(0.5, y, z), dom(3:6));
pass(31) = norm(fChebfun2 - f2D) < 100*tol;

% APA's bug report for a function that has y-rank 1:
ff = @(x, y, z) cos(x + z).*sin(z).*exp(y);
f = chebfun3(ff, dom);
f2D = f(0.5, :, :);
fChebfun2 = chebfun2(@(y,z) ff(0.5, y, z), dom(3:6));
pass(32) = norm(fChebfun2 - f2D) < 100*tol;

% Similar to APA's bug report for a function that has z-rank 1:
ff = @(x, y, z) cos(x + y).*sin(y).*exp(z);
f = chebfun3(ff, dom);
f2D = f(0.5, :, :);
fChebfun2 = chebfun2(@(y,z) ff(0.5, y, z), dom(3:6));
pass(33) = norm(fChebfun2 - f2D) < 100*tol;

% Fixed y
ff = @(x,y,z) sin(x+y+z) + x + y; % A function with different variable-ranks
f = chebfun3(ff, dom);
f2D = f(:, -3, :);
fChebfun2 = chebfun2(@(x,z) ff(x, -3, z), [dom(1:2) dom(5:6)]);
pass(34) = norm(fChebfun2 - f2D) < 100*tol;

% A function that has x-rank 1:
ff = @(x, y, z) cos(y + z).*sin(z).*exp(x);
f = chebfun3(ff, dom);
f2D = f(:, -3, :);
fChebfun2 = chebfun2(@(x, z) ff(x, -3, z), [dom(1:2) dom(5:6)]);
pass(35) = norm(fChebfun2 - f2D) < 100*tol;

% A function that has y-rank 1:
ff = @(x, y, z) cos(x + z).*sin(z).*exp(y);
f = chebfun3(ff, dom);
f2D = f(:, -3, :);
fChebfun2 = chebfun2(@(x, z) ff(x, -3, z), [dom(1:2) dom(5:6)]);
pass(36) = norm(fChebfun2 - f2D) < 100*tol;

% A function that has z-rank 1:
ff = @(x, y, z) cos(x + y).*sin(y).*exp(z);
f = chebfun3(ff, dom);
f2D = f(:, -3, :);
fChebfun2 = chebfun2(@(x, z) ff(x, -3, z), [dom(1:2) dom(5:6)]);
pass(37) = norm(fChebfun2 - f2D) < 100*tol;

% Fixed z
%ff = @(x,y,z) sin(x+y+z);
ff = @(x,y,z) sin(x+y+z) + z + y; % A function with different variable-ranks
f = chebfun3(ff, dom);
f2D = f(:, :, 7);
fChebfun2 = chebfun2(@(x,y) ff(x, y, 7), dom(1:4));
pass(38) = norm(fChebfun2 - f2D) < 100*tol;

% A function that has x-rank 1:
ff = @(x, y, z) cos(y + z).*sin(z).*exp(x);
f = chebfun3(ff, dom);
f2D = f(:, :, 7);
fChebfun2 = chebfun2(@(x, y) ff(x,y, 7), dom(1:4));
pass(39) = norm(fChebfun2 - f2D) < 100*tol;

% A function that has y-rank 1:
ff = @(x, y, z) cos(x + z).*sin(z).*exp(y);
f = chebfun3(ff, dom);
f2D = f(:, :, 7);
fChebfun2 = chebfun2(@(x, y) ff(x, y, 7), dom(1:4));
pass(40) = norm(fChebfun2 - f2D) < 100*tol;

% A function that has z-rank 1:
ff = @(x, y, z) cos(x + y).*sin(y).*exp(z);
f = chebfun3(ff, dom);
f2D = f(:, :, 7);
fChebfun2 = chebfun2(@(x, y) ff(x, y, 7), dom(1:4));
pass(41) = norm(fChebfun2 - f2D) < 100*tol;

% Fixed x and y
ff = @(x,y,z) sin(x+y+z);
f = chebfun3(ff);
f1 = f(0.5, 0.5, :);
fChebfun = chebfun(@(z) sin(1+z));
pass(42) = norm(fChebfun - f1) < 100*tol;

% Fixed x and z
ff = @(x,y,z) sin(x+y+z);
f = chebfun3(ff);
f1 = f(0.5, :, 0.5);
fChebfun = chebfun(@(z) sin(1+z));
pass(43) = norm(fChebfun - f1) < 100*tol;

% Fixed y and z
ff = @(x,y,z) sin(x+y+z);
f = chebfun3(ff);
f1 = f(:, 0.5, 0.5);
fChebfun = chebfun(@(z) sin(1+z));
pass(44) = norm(fChebfun - f1) < 100*tol;

% No fixed variables
ff = @(x,y,z) sin(x+y+z);
f = chebfun3(ff);
fNew = f(:, :, :);
pass(45) = norm(fNew - f) < 100*tol;

%% Evaluate at parametric 1D chebfuns
f = chebfun3(@(x,y,z) x+y.*z);
curve = chebfun(@(t) [cos(t) sin(t) t/(8*pi)], [0, 8*pi]);
f1D = f(curve(:, 1), curve(:, 2), curve(:, 3));
fExact = curve(:, 1) + curve(:, 2).*curve(:, 3);
pass(46) = norm(fExact - f1D) < 100*tol;

%% Test from #1898:
ff = @(x, y, z) sin(x + 2*y + z);
f = chebfun3(ff);
x = linspace(-1, 1, 2).';
y = linspace(-1, 1, 4).';
z = 0;
[xx, yy, zz] = meshgrid(x, y, z);
vals = f(xx, yy, zz);
pass(47) = all(size (xx) == size(vals));
pass(48) = norm(ff(xx,yy,zz) - vals) < tol;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% pop_kaiserbeta() - Estimate Kaiser window beta
%
% Usage:
%   >> [beta, dev] = pop_kaiserbeta; % pop-up window mode
%   >> beta = pop_kaiserbeta(dev);
%
% Inputs:
%   dev       - scalar maximum passband deviation/ripple
%
% Output:
%   beta      - scalar Kaiser window beta
%   dev       - scalar maximum passband deviation/ripple
%
% References:
%   [1] Proakis, J. G., & Manolakis, D. G. (1996). Digital Signal
%       Processing: Principles, Algorithms, and Applications (3rd ed.).
%       Englewood Cliffs, NJ: Prentice-Hall
%
% Author: Andreas Widmann, University of Leipzig, 2005
%
% See also:
%   pop_firws, firws, pop_firwsord, windows

%123456789012345678901234567890123456789012345678901234567890123456789012

% Copyright (C) 2005 Andreas Widmann, University of Leipzig, widmann@uni-leipzig.de
%
% 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 2 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, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function [beta, dev] = pop_kaiserbeta(dev)

    beta = [];

    if nargin < 1 || isempty(dev)
        drawnow;
        uigeom = {[1 1]};
        uilist = {{'style' 'text' 'string' 'Max passband deviation/ripple:'} ...
                  {'style' 'edit' 'string' ''}};
        result = inputgui(uigeom, uilist, 'pophelp(''pop_kaiserbeta'')', 'Estimate Kaiser window beta -- pop_kaiserbeta()');
        if length(result) == 0, return, end
        if ~isempty(result{1})
            dev = str2num(result{1});
        else
            error('Not enough input arguments.');
        end
    end

    devdb = -20 * log10(dev);
    if devdb > 50
        beta = 0.1102 * (devdb - 8.7);
    elseif devdb >= 21
        beta = 0.5842 * (devdb - 21)^0.4 + 0.07886 * (devdb - 21);
    else
        beta = 0;
    end

end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [EB] = PB2EB(PB)
% Convert computery things from petabytes to exabytes.
% Chad A. Greene 2012
EB = PB*2^-10 ;


text:
function pdf = truncated_normal_a_pdf ( x, mu, sigma, a )

%*****************************************************************************80
%
%% TRUNCATED_NORMAL_A_PDF evaluates the lower truncated Normal PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 August 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the argument of the PDF.
%
%    Input, real MU, S, the mean and standard deviation of the
%    parent Normal distribution.
%
%    Input, real A, the lower truncation limit.
%
%    Output, real PDF, the value of the PDF.
%
  alpha = ( a - mu ) / sigma;
  xi = ( x - mu ) / sigma;

  alpha_cdf = normal_01_cdf ( alpha );
  xi_pdf = normal_01_pdf ( xi );

  pdf = xi_pdf / ( 1.0 - alpha_cdf ) / sigma;

  return
end

repo: jburkardt-m


text:
function [ferns,hsPr] = fernsClfReplaceFeature( data, hs, ferns, features_remove, varargin )
% Train random fern classifier.
%
% See "Fast Keypoint Recognition in Ten Lines of Code" by Mustafa Ozuysal,
% Pascal Fua and Vincent Lepetit, CVPR07.
%
% Dimensions:
%  M - number ferns
%  S - fern depth
%  F - number features
%  N - number input vectors
%  H - number classes
%
% USAGE
%  [ferns,hsPr] = fernsClfTrain( data, hs, [varargin] )
%
% INPUTS
%  data     - [NxF] N length F feature vectors
%  hs       - [Nx1] target output labels in [1,H]
%  varargin - additional params (struct or name/value pairs)
%   .S        - [10] fern depth (ferns are exponential in S)
%   .M        - [50] number of ferns to train
%   .thrr     - [0 1] range for randomly generated thresholds
%   .bayes    - [1] if true combine probs using bayes assumption
%
% OUTPUTS
%  ferns    - learned fern model w the following fields
%   .fids     - [MxS] feature ids for each fern for each depth
%   .thrs     - [MxS] threshold corresponding to each fid
%   .pFern    - [2^SxHxM] learned log probs at fern leaves
%   .bayes    - if true combine probs using bayes assumption
%   .inds     - [NxM] cached indices for original training data
%   .H        - number classes
%  hsPr     - [Nx1] predicted output labels
%
% EXAMPLE
%  N=5000; H=5; d=2; [xs0,hs0,xs1,hs1]=demoGenData(N,N,H,d,1,1);
%  fernPrm=struct('S',4,'M',50,'thrr',[-1 1],'bayes',1);
%  tic, [ferns,hsPr0]=fernsClfTrain(xs0,hs0,fernPrm); toc
%  tic, hsPr1 = fernsClfApply( xs1, ferns ); toc
%  e0=mean(hsPr0~=hs0); e1=mean(hsPr1~=hs1);
%  fprintf('errors trn=%f tst=%f\n',e0,e1); figure(1);
%  subplot(2,2,1); visualizeData(xs0,2,hs0);
%  subplot(2,2,2); visualizeData(xs0,2,hsPr0);
%  subplot(2,2,3); visualizeData(xs1,2,hs1);
%  subplot(2,2,4); visualizeData(xs1,2,hsPr1);
%
% See also fernsClfApply, fernsInds
%
% Piotr's Image&Video Toolbox      Version 2.50
% Copyright 2010 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Lesser GPL [see external/lgpl.txt]

% get additional parameters and check dimensions
dfs={'features_replace',[],'doreplace',[],'fids_replace',[]};
[features_replace,doreplace,fids_replace]=getPrmDflt(varargin,dfs,1);
[N,F]=size(data); assert(length(hs)==N);
M=size(ferns.inds,2); S=log2(size(ferns.pFern,1));
H=max(hs); assert(all(hs>0)); assert(S<=20);
assert(size(ferns.inds,1)==N);

% indices we want to remove
if isempty(doreplace),
  doreplace = ismember(ferns.fids,features_remove);
else
  assert(numel(doreplace)==M*S);
end
nreplace = nnz(doreplace);

% features to replace
if isempty(fids_replace),
  % allowed features to replace with
  if isempty(features_replace),
    features_replace = uint32(setdiff(1:F,features_remove));
  else
    features_replace = uint32(features_replace);
  end
  fids_replace = features_replace(floor(rand(nreplace,1)*numel(features_replace)+1));
else
  assert(numel(fids_replace)==nreplace);
end

% replace features in fids
ferns.fids(doreplace) = fids_replace;

% update inds to reflect these replacements
ferns.inds = fernsInds_update(data,ferns.fids,ferns.thrs,ferns.inds,doreplace);

% get counts for each leaf for each class for each fern
% KB: use hist to avoid inner loop
pFern = nan(2^S,H,M);
edges = 1:2^S;
for m = 1:M,
  for h = 1:H,
    pFern(:,h,m) = histc(ferns.inds(hs==h,m),edges);
  end
end
pFern = pFern + ferns.bayes;

% KB: store the unnormalized counts
ferns.counts = pFern;

% convert fern leaf class counts into probabilities
if( ferns.bayes<=0 )
  pFern = bsxfun(@rdivide,pFern,sum(pFern,2));
else
  pFern = bsxfun(@rdivide,pFern,sum(pFern,1));
  pFern=log(pFern);
end

% store pFern and compute output values
ferns.pFern=pFern; clear pFern;
if(nargout==2), hsPr=fernsClfApply([],ferns,ferns.inds); end

end

author: kristinbranson

repo: JAABA


text:
function data = spm_chi2_plot(action,arg1,arg2,arg3,arg4)
% Display a plot showing convergence of an optimization routine.
% FORMAT spm_chi2_plot('Init',title,xlabel,ylabel)
% Initialises the plot in the 'Interactive' window.
%
% FORMAT spm_chi2_plot('Set',value)
% Updates the plot.
%
% FORMAT spm_chi2_plot('Clear')
% Clears the 'Interactive' window.
%
%-----------------------------------------------------------------------
% @(#)spm_chi2_plot.m	1.2 John Ashburner 97/09/11

global pb_pointer pb_name ax
%-----------------------------------------------------------------------
if (nargin == 0)
	spm_chi2_plot('Init');
else
	% initialize
	%---------------------------------------------------------------
	if (strcmp(lower(action),'init'))
		if (nargin<4)
			arg3 = 'Iteration #';
			if (nargin<3)
				arg2 = 'Chi-squared';
				if (nargin<2)
					arg1 = 'Optimizing';
				end
			end
		end
		fg = spm_figure('FindWin','Interactive');
		if ~isempty(fg)
			pb_pointer = get(fg,'Pointer');
			pb_name    = get(fg,'Name');
			spm_chi2_plot('Clear');
			set(fg,'Pointer','watch');
			set(fg,'Name',pb_name);
			ax = axes('Position', [0.15 0.1 0.8 0.75],...
				'Box', 'on','Parent',fg);
			lab = get(ax,'Xlabel');
			set(lab,'string',arg3,'FontSize',10);
			lab = get(ax,'Ylabel');
			set(lab,'string',arg2,'FontSize',10);
			lab = get(ax,'Title');
			set(lab,'string',arg1);
			line('Xdata',[], 'Ydata',[],...
				'LineWidth',2,'Tag','Chi2Plot','Parent',ax);
			drawnow;
		end

	% reset
	%---------------------------------------------------------------
	elseif (strcmp(lower(action),'set'))
		if (nargin<2)
			arg1 = 0;
		end
		F = spm_figure('FindWin','Interactive');
		br = findobj(F,'Tag','Chi2Plot');
		if (~isempty(br))
			xd = get(br,'Xdata');
			yd = [get(br,'Ydata') arg1];
			xd = [xd (length(xd)+1)];
			set(br,'Ydata',yd,'Xdata',xd);
			drawnow;
		end

	% clear
	%---------------------------------------------------------------
	elseif (strcmp(lower(action),'clear'))
		fg = spm_figure('FindWin','Interactive');
		spm_figure('Clear',fg);
		set(fg,'Pointer',pb_pointer);
		set(fg,'Name',pb_name);
		drawnow;
	elseif (strcmp(lower(action),'data'))
		F = spm_figure('FindWin','Interactive');
		br = findobj(F,'Tag','Chi2Plot');
		if (~isempty(br))
			data = get(br,'Ydata');
		else
			data = [];
		end
	end
end