text:
classdef ConstantDragCoeffModel < AbstractDragCoefficientModel
    %ConstantDragCoeffModel Summary of this class goes here
    %   Detailed explanation goes here

    properties
        CdA(1,1) double
    end

    properties(Constant)
        enum = DragCoefficientModelEnum.Constant
    end

    methods
        function obj = ConstantDragCoeffModel(CdA)
            obj.CdA = CdA;
        end

        function CdA = getDragCoeff(obj, ut, rVect, vVect, bodyInfo, mass, altitude, pressure, density, vVectECEFMag, totalAoA, aoa, sideslip)
            arguments
                obj(1,1) ConstantDragCoeffModel
                ut(1,1) double
                rVect(3,1) double
                vVect(3,1) double
                bodyInfo(1,1) KSPTOT_BodyInfo
                mass(1,1) double
                altitude(1,1) double
                pressure(1,1) double
                density(1,1) double
                vVectECEFMag(1,1) double
                totalAoA(1,1) double = 0;
                aoa(1,1) double = 0;
                sideslip(1,1) double = 0;
            end

            CdA = obj.CdA;
        end

        function tf = usesTotalAoA(obj)
            tf = false;
        end

        function tf = usesAoaAndSideslip(obj)
            tf = false;
        end

        function useTf = openEditDialog(obj, lvdData)
            out = AppDesignerGUIOutput({false});
            lvd_EditConstDragPropertiesGUI_App(obj, lvdData, out);
            useTf = out.output{1};
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function stat_ds=cosmo_stat(ds, stat_name, output_stat_name)
% compute t-test or F-test (ANOVA) statistic
%
% stat_ds=cosmo_stats(ds, stat_name[, output_stat_name])
%
% Inputs:
%   ds                dataset struct with
%                       .samples PxQ, for P observations on Q features
%                       .sa.targets Px1 observation conditions (classes)
%                       .sa.chunks  Px1 observation chunks (e.g. subjects).
%   stat_name         One of [*]:
%                     't' : one-sample t-test against zero (nclasses==1),
%                           or paired t-test (nclasses==2)
%                     't2': two-sample t-test with equal variance,
%                           (nclasses==2) contrasting samples with unq(1)
%                           minus unq(2) where unq=unique(ds.sa.targets)
%                     'F' : one-way ANOVA or repeated measures ANOVA
%                           (nclasses>=2)
%                     [*] nclasses is the number of unique values in
%                         ds.sa.targets
%   output_stat_name  (optional) 'left', 'right', 'both', 'z', 'p', or
%                      the empty string '' (default).
%                     - 'left', 'right', and 'both' return a p-value with
%                        the specified tail.
%                     - 'p' returns a p-value, with tail='right' if
%                        stat_name='F' and tail='both' otherwise.
%                     - 'z' returns a z-score corresponding to the p-value
%                     - '' (empty) returns the t or F statistic.
%                     Missing values can be indicated by NaNs; if these are
%                     present and a p-value or z-score is returned, then
%                     these values are computed with a possible variable
%                     number of degrees of freedom across features.
%
% Returns:
%   stat_ds          dataset struct with fields:
%     .samples       1xQ statistic value, or (if output_stat_name is
%                    non-empty) z-score or p-value. See the Notes below
%                    for interpreting p-values.
%     .sa            struct with field:
%        .stats      One of 'Ftest(df1,df2)', 'Ttest(df)', 'Zscore()', or
%                    'Pval()', where df* are the degrees of freedom
%     .[f]a          identical to the input ds.[f]a, if present.
%
% Examples:
%     % one-sample t-test
%     % make a simple dataset
%     ds=struct();
%     ds.samples=reshape(mod(1:7:(12*3*7),13)',[],3)-3;
%     ds.sa.targets=ones(12,1);
%     ds.sa.chunks=(1:12)';
%     cosmo_disp(ds.samples);
%     %||   [ -2         4        -3
%     %||      5        -2         4
%     %||     -1         5        -2
%     %||      :         :         :
%     %||      9         2         8
%     %||      3         9         2
%     %||     -3         3         9 ]@12x3
%     %
%     % run one-sample t-test
%     s=cosmo_stat(ds,'t');
%     cosmo_disp(s);
%     %|| .samples
%     %||   [ 2.49      3.36      2.55 ]
%     %|| .sa
%     %||   .stats
%     %||     { 'Ttest(11)' }
%     %
%     % compute z-score of t-test
%     s=cosmo_stat(ds,'t','z');
%     cosmo_disp(s);
%     %|| .samples
%     %||   [ 2.17      2.73      2.21 ]
%     %|| .sa
%     %||   .stats
%     %||     { 'Zscore()' }
%     %
%     % compute (two-tailed) p-value of t-test
%     s=cosmo_stat(ds,'t','p');
%     cosmo_disp(s);
%     %|| .samples
%     %||   [ 0.03   0.00633    0.0268 ]
%     %|| .sa
%     %||   .stats
%     %||     { 'Pval()' }
%     %
%     % compute left-tailed p-value of t-test
%     s=cosmo_stat(ds,'t','left');
%     cosmo_disp(s);
%     %|| .samples
%     %||   [ 0.985     0.997     0.987 ]
%     %|| .sa
%     %||   .stats
%     %||     { 'Pval()' }
%
%     % one-way ANOVA
%     % each observation is independent and thus each chunk is unique;
%     % there are three conditions with four observations per condition
%     ds=struct();
%     ds.samples=reshape(mod(1:7:(12*3*7),13)',[],3)-3;
%     ds.sa.targets=repmat(1:3,1,4)';
%     ds.sa.chunks=(1:12)';
%     s=cosmo_stat(ds,'F');
%     cosmo_disp(s);
%     %|| .samples
%     %||   [ 0.472    0.0638      0.05 ]
%     %|| .sa
%     %||   .stats
%     %||     { 'Ftest(2,9)' }
%     % compute z-score
%     s=cosmo_stat(ds,'F','z'); % convert to z-score
%     cosmo_disp(s);
%     %|| .samples
%     %||   [ -0.354     -1.54     -1.66 ]
%     %|| .sa
%     %||   .stats
%     %||     { 'Zscore()' }
%
%
%     % two-sample t-test
%     % each observation is independent and thus each chunk is unique;
%     % there are two conditions with four observations per condition
%     ds=struct();
%     ds.samples=reshape(mod(1:7:(12*3*7),13)',[],3)-3;
%     ds.sa.targets=repmat(1:2,1,6)';
%     ds.sa.chunks=(1:12)';
%     s=cosmo_stat(ds,'t2','p'); % return p-value
%     cosmo_disp(s);
%     %|| .samples
%     %||   [ 0.0307  0.000242  7.07e-05 ]
%     %|| .sa
%     %||   .stats
%     %||     { 'Pval()' }
%     %
%     % for illustration, this test gives the same p-values as a
%     % repeated measures ANOVA
%     s=cosmo_stat(ds,'F','p');
%     cosmo_disp(s);
%     %|| .samples
%     %||   [ 0.0307  0.000242  7.07e-05 ]
%     %|| .sa
%     %||   .stats
%     %||     { 'Pval()' }
%
% Notes:
%  - If output_stat_name is not provided or empty, then this function runs
%    considerably faster than the builtin matlab functions
%    (ttest, ttest2, or anova1).
%  - When output_stat_name=='p' then the p-values returned are the same as
%    the builtin matlab functions anova1, ttest, and ttest2 with the
%    default tails.
%  - To run a one-sample t-tests against x (if x~=0), one has to
%    subtract x from ds.samples before using ds as input to this function
%  - The .sa.chunks and .sa.targets determine which test is performed:
%    * statname=='t': all chunks are unique => one-sample t-test
%                   : each chunk present twice => paired-sample t-test
%    * statname=='F': all chunks are unique => one-way ANOVA
%                   : each chunk present N times => repeated measures ANOVA
%    See cosmo_montecarlo_cluster_stat for examples on how .sa.targets and
%    .sa.chunks should be set for different statistics.
%  - Missing values can be indicated by NaNs, and if the output is a
%    p-value (or a z-score based on the p-value), then the output is
%    computed for different features using varying degrees of freedom.
%    For example, if the dataset has 10 samples and a one-sample t-test is
%    used, z-scores for samples with no NaN values is based on the
%    t-statistic with df=9, but those with two missing values (NaN values)
%    are based on the t-statistic with df=7. A use case is computing fMRI
%    group statistics where overlap across brains is not perfect at the
%    voxel-by-voxel level, in as imilar approach as AFNI's 3dttest++ with
%    the '-toz' option.
%  - This function computes feature-wise statistics that are not corrected
%    for multiple comparisons. To correct for multiple comparisons, see
%    cosmo_montecarlo_cluster_stat.
%
% See also: anova1, ttest, ttest2, cosmo_montecarlo_cluster_stat
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

    if nargin<3
        output_stat_name='';
    end

    [output_stat_name,tail]=get_stat_definition(stat_name,...
                                                output_stat_name);
    [samples,targets,chunks,type]=get_descriptors(ds);

    % run specified helper function
    if isfield(ds.sa,'contrast')
        contrast=ds.sa.contrast;
    else
        contrast=[];
    end

    % ensure stat_name matches the number of classes
    verify_class_count(stat_name,targets);

    % get ttest1, ttest2 or Ftest function handle
    stat_func=get_stat_func(stat_name);

    % apply statistic function handle
    [stat,df_struct,stat_label]=apply_stat_func(stat_func,...
                                    samples,targets,chunks,type,contrast);

    % convert to p-value or z-score, if necessary
    [stat,stat_label]=compute_output_stat(stat,df_struct,stat_label,...
                                                tail,output_stat_name);

    % store output
    stat_ds=struct();
    if isfield(ds,'a'), stat_ds.a=ds.a; end
    if isfield(ds,'fa'), stat_ds.fa=ds.fa; end
    stat_ds.samples=stat;
    stat_ds.sa.stats={stat_label};

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% helper functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [stat,df_struct,stat_label]=apply_stat_func(stat_func,...
                                                samples,targets,chunks,...
                                                type,contrast)
% applies the stat_func to samples, taking into account missing values
% (indicated by NaN) in samples
    samples_nan_msk=isnan(samples);
    has_missing=any(any(samples_nan_msk,1) & ~all(samples_nan_msk,1));
    if ~has_missing
        [stat,df_struct,stat_label]=apply_stat_func_no_nans(stat_func,...
                                                samples,targets,chunks,...
                                                type,contrast);
        return;
    end

    [nsamples,nfeatures]=size(samples);
    stat=zeros(1,nfeatures);

    zeroed_samples=samples;
    zeroed_samples(~samples_nan_msk)=0;

    nan_counts=sum(samples_nan_msk,1);
    unq_counts=unique(nan_counts);

    % build indicator matrix, so that a value in zeroed_samples(i,j)=v
    % - v=0 means that samples(i,j) is not NaN
    % - v=k, k>0 means that samples(i,j) is the k-th NaN value in that
    % column
    for c=1:numel(unq_counts)
        count=unq_counts(c);
        col_msk=nan_counts==count;

        msk=bsxfun(@and, samples_nan_msk, col_msk);
        zeroed_samples(msk)=repmat((1:count)',1,sum(col_msk));
    end

    % using the indicator matrix, split up the samples in blocks so that
    % in each block, the same rows in samples are NaN. Each block is
    % processed seperately in the 'for' loop below.

    [indices, zeroed_unq]=cosmo_index_unique(zeroed_samples');

    nclasses=max(targets);

    for k=1:numel(indices)
        zeroed=zeroed_unq(k,:)';
        cols=indices{k};

        nan_msk=zeroed>0;

        % keep samples where there are non-NaN values for the corresponding
        % chunk
        keep_msk=any(bsxfun(@eq,chunks(:),chunks(~nan_msk)'),2);
        keep_samples=samples(keep_msk,cols);

        % ensure chunks are numbered 1:max(chunks)
        [unused,unused,keep_chunks]=unique(chunks(keep_msk));

        % select targets
        keep_targets=targets(keep_msk);

        % compute t or F statistic
        [stat_k,df,stat_label_k]=stat_func(keep_samples,...
                                        keep_targets,...
                                        keep_chunks,...
                                        type,contrast);

        % if for some samples in a chunk the value was NaN and for others
        % it was not NaN, then the output is NaN
        has_mixed_nan_and_non_nan=any(nan_msk(keep_msk));

        % if a target value was completely non-present, then
        % the output is NaN
        keep_targets_msk=false(1,nclasses);
        keep_targets_msk(targets(keep_msk))=true;

        if has_mixed_nan_and_non_nan || any(~keep_targets_msk)
            stat_k(:)=NaN;
        end

        % store stat val
        stat(cols)=stat_k;

        if k==1
            % allocate space
            df_count=numel(df);
            df_matrix=zeros(df_count,nfeatures);
            stat_label=stat_label_k;
        else
            % stat label must always be the same
            assert(isequal(stat_label,stat_label_k));
        end

        ncols=numel(cols);

        df_matrix(:,cols)=repmat(df(:),1,ncols);
    end

    single_sample=samples(:,1);
    single_sample(:)=0;

    [unused,max_df]=stat_func(single_sample,targets,chunks,...
                                    type,contrast);

    df_struct=struct();
    df_struct.max_df=max_df(:);
    df_struct.feature_wise_df=df_matrix;

function [stat,df_struct,stat_label]=apply_stat_func_no_nans(stat_func,...
                                                samples,targets,chunks,...
                                                type,contrast)
    % much faster computation of statistics, if there are no NaNs
    [stat,max_df,stat_label]=stat_func(samples,targets,chunks,...
                                             type,contrast);

    nfeatures=size(samples,2);

    df_struct=struct();
    df_struct.max_df=max_df(:);
    df_struct.feature_wise_df=repmat(max_df(:),1,nfeatures);

function f=get_stat_func(stat_name)
    stat_name2func=struct();
    stat_name2func.t=@ttest1_wrapper;
    stat_name2func.t2=@ttest2_wrapper;
    stat_name2func.F=@ftest_wrapper;

    assert(isfield(stat_name2func,stat_name))

    f=stat_name2func.(stat_name);

function verify_class_count(stat_name,targets)
    stat_name2interval=struct();
    stat_name2interval.t=[1,2];
    stat_name2interval.t2=[2,2];
    stat_name2interval.F=[2,Inf];

    if ~isfield(stat_name2interval,stat_name)
        error('illegal statname ''%s'', supported are: %s',stat_name,...
                    cosmo_strjoin(fieldnames(stat_name2interval),', '));
    end

    valid_interval=stat_name2interval.(stat_name);

    count=max(targets);
    assert(count==numel(unique(targets)));

    min_count=valid_interval(1);
    if count<min_count
        error('statname ''%s'' requires at least %d unique targets',...
                        stat_name,min_count);
    end

    max_count=valid_interval(2);
    if count>max_count
        error('statname ''%s'' requires at most %d unique targets',...
                        stat_name,max_count);
    end

function [stat,stat_label]=compute_output_stat(stat,df_struct,stat_name,...
                                                 tail,output_stat_name)
% transform output is required

    if isempty(output_stat_name)
        % raw t or F value
        stat_fa_name=stat_name;

        % features with missing values (as indicated by df values that
        % are different from the maximum df possible) are set to NaN
        df_is_max_msk=bsxfun(@eq,df_struct.max_df,...
                                       df_struct.feature_wise_df);
        has_missing=~all(df_is_max_msk,1);
        stat(has_missing)=NaN;

        % set degrees of freedom
        df_str=arrayfun(@(x) sprintf('%d',x), df_struct.max_df,...
                                    'UniformOutput',false);
        stat_label=sprintf('%s(%s)',stat_fa_name,...
                                    cosmo_strjoin(df_str,','));
    else
        % transform to left-tailed p-value for each unique combination
        % of degrees of freedom
        stat=apply_cdf_wrapper_different_dfs(stat_name,stat,df_struct);

        switch output_stat_name
            case 'z'
                % transform to z-score
                stat=cosmo_norminv(stat);
                stat_fa_name='Zscore';
            case 'p'
                switch tail
                    case 'left'
                        % do nothing
                    case 'right'
                        % invert p-value
                        stat=1-stat;
                    case 'both'
                        % take whichever tail is more extreme
                        stat=(.5-abs(stat-.5))*2;
                    otherwise
                        assert(false,'this should not happen');
                end
                stat_fa_name='Pval';
            otherwise
                error('illegal output type %s', output_stat_name);
        end

        stat_label=sprintf('%s()',stat_fa_name);
    end

function cdf_stat=apply_cdf_wrapper_different_dfs(stat_name,stat,df_struct)
    % applies the cdf wrapper to each unique combination of degrees of
    % freedom stored in df_struct.feature_wise_df

    [idx, unique_dfs]=cosmo_index_unique(df_struct.feature_wise_df');

    % allocate space for output
    cdf_stat=zeros(size(stat))+NaN;

    % compute result
    for k=1:numel(idx)
        cols=idx{k};
        df_cell=num2cell(unique_dfs(k,:));

        col_stat=stat(cols);
        msk=~isnan(col_stat);

        cdf_stat(cols(msk))=cdf_wrapper(stat_name,col_stat(msk),df_cell{:});
    end

function y=cdf_wrapper(name, x, df1, df2)
    ensure_has_stats_functions();
    if ~(df1>0)
        y=zeros(size(x))+NaN;
        return;
    end

    switch name
        case 'Ttest'
            assert(nargin==3);
            y=tcdf(x, df1);

        case 'Ftest'
            assert(nargin==4);
            if ~(df2>0)
                y=zeros(size(x))+NaN;
            return;
            end

            y=fcdf(x, df1, df2);

        otherwise
            assert(false);
    end

function [stat,df,stat_label]=ttest1_wrapper(samples,targets,chunks,...
                                                           type,contrast)
    if ~isempty(contrast)
        error('contrast is not supported for t-stat');
    end

    nclasses=max(targets);
    if nclasses==2
        samples=pairwise_differences(samples,targets,chunks);
        nclasses=1;
    end

    if nclasses~=1
        error('t-stat: expected 1 or 2 classes, found %d',...
                    nclasses);
    end

    [stat,df]=quick_ttest(samples);
    stat_label='Ttest';

function [stat,df,stat_label]=ttest2_wrapper(samples,targets,chunks,...
                                                           type,contrast)
    if ~isempty(contrast)
        error('contrast is not supported for t-stat');
    end

    if strcmp(type,'within')
        error(['ttest2 stat: the values in chunks and targets suggest '...
                    'a within-subject design. If you want to '...
                    'run a paired-test, use ''t'' (not ''t2'') ',...
                    'as the second argument']);
    end

    m1=targets==1;
    m2=targets==2;

    [stat,df]=quick_ttest2(samples(m1,:),...
                          samples(m2,:));
    nclasses=max(targets);
    if nclasses~=2
        % missing targets, set all to NaN
        assert(nclasses<2)
        stat(:)=NaN;
    end

    stat_label='Ttest';

function [stat,df,stat_label]=ftest_wrapper(samples,targets,chunks,...
                                                            type,contrast)
    nclasses=max(targets);

    if nclasses<2
        error('F stat: expected >=2 classes, found %d',nclasses);
    end

    switch type
        case 'between'
            [stat,df]=quick_ftest_between(samples, targets, ...
                                        chunks,contrast);
        case 'within'
            [stat,df]=quick_ftest_within(samples, targets, ...
                                        chunks,contrast);

    end
    stat_label='Ftest';

function [f,df]=quick_ftest_between(samples,targets,chunks,contrast)
    % one-way ANOVA
    has_contrast=~isempty(contrast);
    contrast_sum=0;

    nclasses=max(targets);

    [ns,nf]=size(samples);
    mu=sum(samples,1)/ns; % grand mean

    b=zeros(nclasses,nf); % between-class sum of squares
    nsc=zeros(nclasses,1);
    wss=0; % within-class sum of squares

    for k=1:nclasses
        msk=k==targets;

        nsc(k)=sum(msk); % number of samples in this class
        sample=samples(msk,:);
        muc=sum(sample,1)/nsc(k); % class mean

        % between- and within-class sum of squares
        if has_contrast
            cmsk=contrast(msk);
            if ~all(cmsk(1)==cmsk)
                error('Contrast has differerent values in level %d',k);
            end
            contrast_sum=contrast_sum+cmsk(1);
            b(k,:)=sum(bsxfun(@times,contrast(msk),mu-muc),1);
        else
            b(k,:)=(mu-muc);
        end
        wss=wss+sum(bsxfun(@minus,muc,sample).^2,1);
    end

    if has_contrast
        if contrast_sum~=0
            error('contrast has sum %d, should be 0', contrast_sum);
        end
        bss=sum(b,1).^2/sum(contrast.^2);
        df1=1;
    else
        bss=sum(bsxfun(@times,nsc,b.^2),1);
        df1=nclasses-1;
    end

    df=[df1,ns-nclasses];

    mbss=bss/df(1);
    mwss=wss/df(2);

    f=zeros(1,nf)+NaN;
    msk=mbss>0;
    f(msk)=mbss(msk)./mwss(msk);

function [f,df]=quick_ftest_within(samples,targets,chunks,contrast)
    % repeated measures anova
    if ~isempty(contrast)
        error('contrast is not supported for within-subject design');
    end

    nchunks=max(chunks);
    nclasses=max(targets);

    nfeatures=size(samples,2);
    gm=mean(samples,1); % grand mean

    sst=zeros(1,nfeatures);
    ssw=zeros(1,nfeatures);
    for k=1:nclasses
        xk=samples(k==targets,:);
        n=size(xk,1);

        if n==0
            ssw(:)=NaN;
            break;
        end

        mu=sum(xk,1)/n;
        sst=sst+n*(gm-mu).^2;
        ssw=ssw+sum(bsxfun(@minus,mu,xk).^2);
    end

    sss=zeros(1,nfeatures);
    for k=1:nchunks
        xk=samples(k==chunks,:);
        n=size(xk,1);
        mu=mean(xk,1);
        sss=sss+n*(gm-mu).^2;
    end

    df1=(nclasses-1);
    mst=sst/df1;

    df2=df1*(nchunks-1);
    sse=ssw-sss;
    mse=sse/df2;

    msk=mse>0;
    f=zeros(1,nfeatures)+NaN;
    f(msk)=mst(msk)./mse(msk);
    df=[df1 df2];

function [t,df]=quick_ttest(x)
    % one-sample t-test against zero

    [ns,nf]=size(x);
    mu=sum(x,1)/ns; % grand mean

    df=ns-1;
    scaling=ns*df;

    % sum of squares
    ss=sum(bsxfun(@minus,x,mu).^2,1);

    t=zeros(1,nf)+NaN;
    msk=ss>0;
    t(msk)=mu(msk).*sqrt(scaling./ss(msk));

function [t,df]=quick_ttest2(x,y)
    % two-sample t-test with equal variance assumption

    [nx,nf]=size(x);
    ny=size(y,1);

    df=nx+ny-2;
    if nx==0 || ny==0
        t=zeros(1,nf)+NaN;
        return;
    end

    mux=sum(x,1)/nx; % mean of class x
    muy=sum(y,1)/ny; % "           " y

    scaling=(nx*ny)*df/(nx+ny);

    % sum of squares
    ss=sum([bsxfun(@minus,x,mux);bsxfun(@minus,y,muy)].^2,1);

    t=zeros(1,nf)+NaN;
    msk=ss>0;
    t(msk)=(mux(msk)-muy(msk)) .* sqrt(scaling./ss(msk));

function ensure_has_stats_functions()
    % - Octave has the required functionality in the octave-forge
    %   statistics toolbox and will raise an error if it is not installed.
    % - Matlab will raise an error message that the statistics toolbox is
    %   required
    persistent cached_has_stat_funcs;

    if isequal(cached_has_stat_funcs,[])
        if cosmo_wtf('is_matlab') ...
                    && isempty(which('tinv')) ...
                    && isempty(which('fcdf'))
            raise_exception_if_absent=true;
            cached_has_stat_funcs=cosmo_check_external('@stats', ...
                                raise_exception_if_absent);
        else
            cached_has_stat_funcs=false;
        end
    end

function [samples,targets,chunks,type]=get_descriptors(ds)
    cosmo_isfield(ds,{'samples','sa.chunks','sa.targets'},true);

    samples=ds.samples;

    % unique targets
    [unused,unusued,targets]=unique(ds.sa.targets);

    % unique chunks
    [unq_chunks,unusued,chunks]=unique(ds.sa.chunks);
    nc=max(chunks);

    if isequal(sort(ds.sa.chunks),unq_chunks)
        type='between';
    else
        combis=(targets-1)*nc+chunks;
        if isequal(sort(combis),(1:numel(combis))')
            type='within';
        else
            error(['Either all chunks must be unique, or each chunk '...
                        'must contain the same targets']);
        end
    end

function delta=pairwise_differences(samples,targets,chunks)
    % for one-sample t-test: compute differences between first and second
    % target
    n=size(samples,1);
    assert(numel(targets)==n);
    assert(numel(chunks)==n);

    assert(mod(n,2)==0);
    n2=n/2;

    idxs=zeros(n2,2);
    for k=1:n
        assert(idxs(chunks(k),targets(k))==0);
        idxs(chunks(k),targets(k))=k;
    end
    assert(isequal(size(idxs),[n2,2]));
    assert(all(idxs(:)>0));

    delta=samples(idxs(:,1),:)-samples(idxs(:,2),:);

function [output_stat_name,tail]=get_stat_definition(stat_name,...
                                                    output_stat_name)
    if any(cosmo_match({'left','right','both'},output_stat_name))
        tail=output_stat_name;
        output_stat_name='p';
    elseif strcmp(output_stat_name,'p')
        switch stat_name
            case 'F'
                tail='right'; % show anova1  behaviour w.r.t. p-values
            otherwise
                tail='both'; % show ttest[2] "                       "
        end
    else
        tail='both';
    end

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function [H,dH] = entub_vbmc(vp,grad_flags,jacobian_flag)
%ENTUB_VBMC Entropy upper bound for variational posterior

% Uses entropy upper bound of multivariate normal approximation

% Check if gradient computation is required
if nargout < 2                              % No 2nd output, no gradients
    grad_flags = false;
elseif nargin < 2 || isempty(grad_flags)    % By default compute all gradients
    grad_flags = true;
end
if isscalar(grad_flags); grad_flags = ones(1,4)*grad_flags; end

% By default assume variational parameters were transformed (before the call)
if nargin < 3 || isempty(jacobian_flag); jacobian_flag = true; end

D = vp.D;           % Number of dimensions
K = vp.K;           % Number of components
mu(:,:) = vp.mu;
sigma(1,:) = vp.sigma;
lambda(:,1) = vp.lambda(:);
w(1,:) = vp.w;

% Check which gradients are computed
if grad_flags(1); mu_grad = zeros(D,K); dS_mu = zeros(D,D,K); else, mu_grad = []; end
if grad_flags(2); sigma_grad = zeros(K,1); else, sigma_grad = []; end
if grad_flags(3); lambda_grad = zeros(D,1); else, lambda_grad = []; end
if grad_flags(4); w_grad = zeros(K,1); dS_w = zeros(D,D,K); else, w_grad = []; end

if K == 1
    % Entropy of single component, uses exact expression
    H = 0.5*D*(1 + log(2*pi)) + D*sum(log(sigma)) + sum(log(lambda));

    if grad_flags(2)
        sigma_grad(:) = D./sigma(:);
    end

    if grad_flags(3)
        % Should be dividing by LAMBDA, see below
        lambda_grad(:) = ones(D,1); % 1./lambda(:);
    end

    if grad_flags(4)
        w_grad = 0;
    end
else

    Mu = sum(bsxfun(@times,vp.w,vp.mu),2);
    Sigma = zeros(D,D);
    delta_mu = bsxfun(@minus,mu,Mu);
    for k = 1:K
        S_k = diag((lambda*sigma(k)).^2) + delta_mu(:,k)*delta_mu(:,k)';
        Sigma = Sigma + w(k)*S_k;
        if grad_flags(4); dS_w(:,:,k) = S_k; end
    end
    L = chol(Sigma);

    H = 0.5*D*(log(2*pi) + 1) + sum(log(diag(L)));

     if any(grad_flags)
         invK = L\(L'\eye(D));

         if grad_flags(1)
             for k = 1:K
                 mu_grad((1:D)+(k-1)*D) = 0.5*w(k).*(sum(bsxfun(@times,invK,delta_mu(:,k)'),2) + sum(bsxfun(@times,invK,delta_mu(:,k)),1)');
             end
         end

         if grad_flags(2)
             Q = sum(sum(invK.*diag(lambda.^2)));
             sigma_grad(:) = Q*(w.*sigma);
         end

         if grad_flags(3)
             lambda_grad(:) = diag(invK).*lambda.^2*sum(w.*(sigma.^2));
         end

         if grad_flags(4)
             for k = 1:K
                w_grad(k) = 0.5*sum(sum(invK.*dS_w(:,:,k)));
             end
         end
     end
end

if nargout > 1
    % Correct for standard log reparameterization of SIGMA
    if jacobian_flag && grad_flags(2)
        sigma_grad = bsxfun(@times,sigma_grad, sigma(:));
    end
    % Correct if NOT using standard log reparameterization of LAMBDA
    if ~jacobian_flag && grad_flags(3)
        lambda_grad = bsxfun(@rdivide,lambda_grad, lambda(:));
    end
    % Correct for standard softmax reparameterization of W
    if jacobian_flag && grad_flags(4)
        eta_sum = sum(exp(vp.eta));
        J_w = bsxfun(@times,-exp(vp.eta)',exp(vp.eta)/eta_sum^2) + diag(exp(vp.eta)/eta_sum);
        w_grad = J_w*w_grad;
    end
    dH = [mu_grad(:); sigma_grad(:); lambda_grad(:); w_grad(:)];
end

end

author: acerbilab

repo: vbmc

save_path: github-repos/MATLAB/acerbilab-vbmc


text:
%% This function will calculate the l2 norm error between the discovered system and actual test data
% Last Updated: 2019/04/22
% Coded By: K

function [Score]=Get_Score(dData_test,Data_test,u,Control,tspan,state0,Shuffle)
% Get the ODE simulation result
Noise=0;

% Using ODEs to get the simulation data
%[dData_Es,Data_Es]=Get_Sim_Data(@(t,z)Sindy_ODE_RHS(t,z,u),state0,u,tspan,Noise,Control,Shuffle);

% Get the score of the result
%[n,m]=size(Data_Es);
%Score=sum(norm(dData_test(1:n,:)-dData_Es))+sum(norm(Data_test(1:n,:)-Data_Es));

%% Using one step prediction to get the simulation data
dData_Es=Sindy_ODE_RHS(0,Data_test',u)';

% Get the score of the result
Score=sum(norm(dData_test-dData_Es));

author: dynamicslab

repo: SINDy-PI


text:
function hRD = plugin_rd
% Reference design definition

%   Copyright 2014-2015 The MathWorks, Inc.

% Call the common reference design definition function
hRD = AnalogDevicesAGCDemo.adrv9361z7035.common.plugin_rd('fmc lvds agc', 'Rx');
AnalogDevicesAGCDemo.adrv9361z7035.ccfmc_lvds_agc.rx.add_rx_io(hRD);

author: analogdevicesinc


text:
function psi_test ( )

%*****************************************************************************80
%
%% PSI_TEST tests R4_PSI and R8_PSI.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 September 2011
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../test_values' );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'PSI_TEST:\n' );
  fprintf ( 1, '  Test PSI_VALUES, R4_PSI, R8_PSI.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '             X         PSI(X)\n' );
  fprintf ( 1, '                    R4_PSI(X)         Diff\n' );
  fprintf ( 1, '                    R8_PSI(X)         Diff\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx1 ] = psi_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fx2 = r4_psi ( single ( x ) );
    fx3 = r8_psi ( x );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  %14.4f  %14.6g\n', x, fx1 );
    fprintf ( 1, '                  %14.6g  %14.6g\n', fx2, abs ( fx1 - fx2 ) );
    fprintf ( 1, '                  %14.6g  %14.6g\n', fx3, abs ( fx1 - fx3 ) );

  end

  rmpath ( '../test_values' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function tests = test_ft_preproc_bandpassfilter

% MEM 1gb
% WALLTIME 00:10:00
% DEPENDENCY ft_preproc_bandpassfilter

if nargout
  % assume that this is called by RUNTESTS
  tests = functiontests(localfunctions);
else
  % assume that this is called from the command line
  fn = localfunctions;
  for i=1:numel(fn)
    feval(fn{i});
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function testOptions(testCase)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nchan   = 8;
nsample = 1000;

dat = randn(nchan, nsample) + 1;
Fs  = 1000;
Fbp = [1 35];

instabilityfix  = [];
df              = [];
wintype         = [];
dev             = [];
plotfiltresp    = [];
usefftfilt      = [];

filttypes = {'brickwall', 'firws' 'fir' 'firls' 'but'};
filtdirs  = {'onepass' 'onepass-reverse' 'twopass' 'twopass-reverse' 'twopass-average' 'onepass-zerophase' 'onepass-reverse-zerophase' 'onepass-minphase'};
filtorders = {[] 1 2 4 8};

result = {};
opts   = {};
for i1 = 1:numel(filttypes)
  if isequal(filttypes{i1}, 'brickwall')
    opts{end+1}   =[filttypes{i1}];
    result{end+1} = ft_preproc_bandpassfilter(dat, Fs, Fbp, [], filttypes{i1}, [], instabilityfix, df, wintype, dev, plotfiltresp, usefftfilt);
    continue;
  elseif isequal(filttypes{i1}, 'but')
    sel2 = 1:5;
    sel3 = 2:4; % order 8 may become unstable, and thus requires an instabilityfix
  elseif isequal(filttypes{i1}, 'firws')
    sel2 = find(~contains(filtdirs, 'twopass')); % this leads to a filter that is not stable
    sel3 = 1; % the low number orders do not make sense at all for the finite impulse filters.
  else
    sel2 = 1:numel(filtdirs);
    sel3 = 1;
  end
  for i2 = sel2
    for i3 = sel3
      tmp = filtorders{i3};
      if isempty(tmp)
        tmp = 'defaultorder';
      elseif isnumeric(tmp)
        tmp = sprintf('order %d', tmp);
      end
      opts{end+1}   = sprintf('%s_%s_%s', filttypes{i1}, filtdirs{i2}, tmp);
      result{end+1} = ft_preproc_bandpassfilter(dat, Fs, Fbp, filtorders{i3}, filttypes{i1}, filtdirs{i2}, instabilityfix, df, wintype, dev, plotfiltresp, usefftfilt);
    end
  end
end

% result{end+1} = ft_preproc_bandpassfilter(dat, Fs, Fbp, 8, 'but'      , 'onepass'                  , 'split', [], [], [], [], []); % this one is instable
% result{end+1} = ft_preproc_bandpassfilter(dat, Fs, Fbp, 8, 'but'      , 'onepass-reverse'          , 'split', [], [], [], [], []);
% result{end+1} = ft_preproc_bandpassfilter(dat, Fs, Fbp, 8, 'but'      , 'twopass'                  , 'split', [], [], [], [], []);
% result{end+1} = ft_preproc_bandpassfilter(dat, Fs, Fbp, 8, 'but'      , 'twopass-reverse'          , 'split', [], [], [], [], []);
% result{end+1} = ft_preproc_bandpassfilter(dat, Fs, Fbp, 8, 'but'      , 'twopass-average'          , 'split', [], [], [], [], []);

% all iterations were done with (slightly) different options, hence the results should not be equal
for i=1:numel(result)
  for j=(i+1):numel(result)
    assert(~isequal(result{i}, result{j}), 'the results %d and %d should not be equal', i, j);
  end
end

for i=1:size(result{1},1)
  for k=1:numel(result)
    for m=(k+1):numel(result)
      b(k,m) = result{k}(i,:)/result{m}(i,:);
      b(m,k) = result{m}(i,:)/result{k}(i,:);
    end
    B(:,:,i) = b;
  end
end
n = numel(result);
figure;imagesc((mean(B,3)+mean(B,3)')/2);
set(gca, 'xtick', 1:n, 'ytick', 1:n, 'xticklabel', opts', 'yticklabel', opts', 'ticklabelinterpreter', 'none');
set(gcf, 'position', [230 47 993 750]);
colorbar
axis equal;axis tight

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%   GeomCG Tensor Completion. Copyright 2013 by
%   Michael Steinlechner
%   Questions and contact: michael.steinlechner@epfl.ch
%   BSD 2-clause license, see LICENSE.txt

disp('                ___   ___      |    geomCG: Tensor completion by Riemannian Optimization')
disp('   _  _  _     / __|/ __|      |    ')
disp('  / \/_\/ \|\/| (__ |(_/\      |    Daniel Kressner')
disp('  \_|\__\_/|  |\___|\___/      |    Michael Steinlechner')
disp('  \_|                          |    Bart Vandereycken')
disp('     ')
disp('geomCG is licensed under a BSD 2-clause license, see LICENSE.txt')
disp('For an explanation of the algorithm, we refer to')
disp('   ')
disp('      D. Kressner, M. Steinlechner, B. Vandereycken:')
disp('      Low-Rank Tensor Completion by Riemannian Optimization.')
disp('      MATHICSE Technical Report 20.2013, June 2013. Submitted to BIT Numerical Mathematics.')

disp('   ')

disp('Compiling mex files...')
mex calcFunction_mex.c
mex calcGradient_mex.c
mex calcInitial_mex.c
mex calcProjection_mex.c
mex getValsAtIndex_mex.c

addpath( cd )
addpath( [cd, filesep, 'examples'] )
disp('Make sure that the Tensor Toolbox is also properly installed and loaded into the current path')
disp('You can download the Tensor Toolbox from Sandia Labs. GeomCG was tested with version 2.5 only!')

disp('Finished. Try out the example code simpleGeomCG.m')

author: andrewssobral


text:
function x = subsref( x, S )

%   Disciplined convex/geometric programming information for SUBSREF:
%      The use of subscripts to extract elements or "slices" of any CVX
%      variable is identical to their use with numeric arrays. All
%      conventions are preserved, including the colon ':' and 'end'.

sx = x.size_;
if numel( S ) == 1 && isequal( S.type, '()' ) && all( strcmp( S.subs, ':' ) ),
    ns = numel(S.subs);
    sx( end + 1 : numel(ns) ) = 1;
    sx = [ sx(1:ns-1), prod(sx(ns:end)) ];
    x = cvx( sx, x.basis_ );
else
    try
        nx = prod( sx );
        ndxs = builtin( 'subsref', reshape( 1 : nx, sx ), S );
    catch errmsg
        throw( errmsg );
    end
    x = cvx( size( ndxs ), x.basis_( :, ndxs ) );
end

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

author: yu-jiang

repo: radpbook

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


text:
function imageHandle = prtPlotUtilImageEvaledClassifier(DS,linGrid,gridSize,cMap)
% Internal function,
% xxx Need Help xxx

nDims = size(linGrid,2);

switch nDims
    case 1
        imageHandle = imagesc(linGrid(:,1),linGrid(:,2),DS);
        set(gca,'YTickLabel',[])
        colormap(cMap)
    case 2
        xx = reshape(linGrid(:,1),gridSize);
        yy = reshape(linGrid(:,2),gridSize);
        imageHandle = imagesc(xx(1,:),yy(:,1),DS);
        colormap(cMap)
    case 3
        xx = reshape(linGrid(:,1),gridSize);
        yy = reshape(linGrid(:,2),gridSize);
        zz = reshape(linGrid(:,3),gridSize);
        imageHandle = slice(xx,yy,zz,DS,max(xx(:)),max(yy(:)),[min(zz(:)),mean(zz(:))]);
        view(3)
        colormap(cMap)
end
axis tight;
axis xy;

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
classdef PTKAirwaysLabelledByBronchus < PTKPlugin
    % PTKAirwaysLabelledByBronchus. Plugin to visualise the airway tree bronchus by bronchus
    %
    %     This is a plugin for the Pulmonary Toolkit. Plugins can be run using
    %     the gui, or through the interfaces provided by the Pulmonary Toolkit.
    %     See PTKPlugin.m for more information on how to run plugins.
    %
    %     Plugins should not be run directly from your code.
    %
    %     PTKAirwaysLabelledByBronchus calls the PTKAirways plugin to segment the
    %     airway tree. It then cretes an output labelled image where each
    %     bronchus is labelled with contrasting colour values, starting with 1
    %     for the trachea. The child segments of each bronchus are labelled
    %     consecutively starting from the parent bronchus plus 1. Eg. if a
    %     parent is labeled 3, then its children will be labelled 4,5,6 etc.
    %     This makes the heirarcy clear from the image.

    %     The output image is intended to be visualised using the Lines
    %     colourmap; for this reason the grey colours (multiples of 7) are
    %     skipped to avoid confusion when used as an overlay on top of a
    %     greyscale image.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    properties
        ButtonText = 'Airways <BR>by bronchus'
        ToolTip = 'Shows airways coloured by bronchus'
        Category = 'Airways'

        AllowResultsToBeCached = true
        AlwaysRunPlugin = false
        PluginType = 'ReplaceOverlay'
        HidePluginInDisplay = false
        FlattenPreviewImage = true
        PTKVersion = '1'
        ButtonWidth = 6
        ButtonHeight = 2
        GeneratePreview = true
        Visibility = 'Developer'
    end

    methods (Static)
        function results = RunPlugin(dataset, ~)
            results = dataset.GetTemplateImage(PTKContext.LungROI);
            airway_results = dataset.GetResult('PTKAirways');
            labeled_region = PTKAirwaysLabelledByBronchus.GetLabeledSegmentedImageFromAirwayTree(airway_results.AirwayTree, results);
            results.ChangeRawImage(labeled_region);
            results.ImageType = PTKImageType.Colormap;
        end
    end

    methods (Static, Access = private)

        function segmented_image = GetLabeledSegmentedImageFromAirwayTree(airway_tree, template)
            airway_tree.AddColourValues(1);
            segmented_image = zeros(template.ImageSize, 'uint8');

            segments_to_do = airway_tree;

            while ~isempty(segments_to_do)
                segment = segments_to_do(end);
                segments_to_do(end) = [];
                segmented_image(template.GlobalToLocalIndices(segment.GetAllAirwayPoints)) = segment.Colour;
                segments_to_do = [segments_to_do segment.Children];
            end
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function loadAlignment
% function loadAlignment
%
% loads bestrotvol.mat, converts it into a 4x4 homog transform matrix
% and incorporates it into mrSESSION.alignment.
%
% 9/17/98 rmk
% 4/29/99 BTB: Removed Inpts and Volpts fields, added tests for mrSESSION
%    field and coords.mat & coranal.mat files.
% 8/2001 djh: updated to mrLoadRet-3.0

global HOMEDIR mrSESSION

pathStr = fullfile(HOMEDIR,'bestrotvol.mat');
if ~exist(pathStr,'file')
    myErrorDlg('No bestrotvol file. Run mrAlign3.');
end
load(pathStr);

if isfield(mrSESSION, 'alignment')
    RESP = questdlg(['mrSESSION has an ''alignment'' field '...
            '(created when rotation was last saved).  '...
            'It was used to build any coords.mat and coranal.mat '...
            'files you may have in Volume, Gray, and Flat directories.'...
            'Permanently replace it with new alignment info?  '], ...
        'Question', 'Yes', 'No', 'Yes');
    if strcmp(RESP, 'Yes')
        rmfield(mrSESSION, 'alignment');
    else
        hmsgbox = msgbox(['Although you have updated the file bestrotvol.mat, ',...
                'you have NOT updated mrSESSION.alignment, nor deleted any old ', ...
                'coords.mat files from Volume, Gray, or Flat directories.'], 'Warning','warn');
        drawnow
        return
    end
end

% If we've made it this far, go ahead an replace mrSESSION.alignment field.
% 4x4 homog transform
mrSESSION.alignment = inplane2VolXform(rot,trans,scaleFac);
saveSession
disp('mrSESSION.alignment has been updated and saved.');

%%% Delete (after confirmation) old files that were built using previous alignment.
cleanAllFlats
cleanGray
cleanVolume

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [classifier] = saveclassifier(name,svs,alpha,xtrain,b,options,sec,t,lsiters,stats)
% {saveclassifier} generates structure with details of a kernel classifier.
%
%      classifier = saveclassifier(name,svs,alpha,xtrain,b,options)
%
%      classifier = saveclassifier(name,svs,alpha,xtrain,b,options,
%                                  sec,t,lsiters,stats)
%
%      name: string identifier of the classifier
%      svs: indices of the training points that are present on the
%           classifier function
%      alpha: coefficient vector (with |svs| elements)
%      xtrain: corresponding training examples
%      b: bias
%      options: the options data structure used to train the classifier,
%               that contains the regularization parameters, kernel
%	            function, kernel parameters etc...
%
%      [optional]
%      sec: training time
%      t: number of iteration performed while training
%      lsi: vector with the number of line searches for each iteration
%      stats: matrix of misc statistics
%
%      classifier: is a data structure with the following fields
%                  classifier.name: name
%                  classifier.svs: svs
%                  classifier.alpha: alpha
%                  classifier.xtrain: xtrain
%                  classifier.b: b
%                  classifier.options: options
%
%                  [if optional arguments are present]
%                  classifier.iters: t
%                  classifier.lsiters: lsiters
%                  classifier.traintime: sec
%                  classifier.stats: stats
%
% Author: Stefano Melacci (2009)
%         mela@dii.unisi.it
%         * based on the code of Vikas Sindhwani, vikas.sindhwani@gmail.com

classifier.name=name;
classifier.svs=svs;
classifier.alpha=alpha;
classifier.b=b;
classifier.xtrain=xtrain;
classifier.options=options;
if exist('t','var') && ~isempty(t), classifier.iters=t; end
if exist('lsiters','var') && ~isempty(lsiters), classifier.lsiters=lsiters; end
if exist('sec','var') && ~isempty(sec), classifier.traintime=sec; end
if exist('stats','var') && ~isempty(stats), classifier.stats=stats; end

author: viggin

repo: domain-adaptation-toolbox


text:
function [OutStruct] = matchObjBtwnTrials(inputImages,varargin)
	% Registers images to a set imaging session then matches objs between sessions to one another and outputs the alignment indicies for a single global cell across sessions. All images cropped to the minimum x,y dimension among all the input imaging session datasets.
	% Biafra Ahanonu
	% started: 2013.10.31
	% inputs
		% inputImages - cell array of [x y nFilters] matrices containing each set of filters, e.g. {imageSet1, imageSet2,...}
	% options
		% inputSignals - cell array of [nFilters frames] matrices containing each set of filter traces
	% outputs
		% OutStruct - structure containing
		% .globalIDs, [M N] matrix with M = number of global IDs and N = each session. Each m,n pair specifies the index of that global obj m in the data of session n. If .globalIDs(m,n)==0, means no match was found.
		% .trialIDs, a cell array that matches each column n in the matrix to a particular id, either automatic or input a cell array of strings specifying what each session is.
		% .coords, a [M C] matrix with M = number of global ID and C = 2 (1st column is x, 2nd column is y coordinates),

	[OutStruct] = ciapkg.classification.matchObjBtwnTrials(inputImages,'passArgs', varargin);
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function [xxx1,xxx2] =dymism(r1,r3,l,r4,r5,choice,rpm,loc)
%-----------------------------------------------------------
format bank;
% clf;
   t = 0:1:360;
   tg = t*pi/180.0;
    cx =[0,0,l,0];
   r2 = sqrt(l^2 - r3^2 + r1^2);
      jj = 0;
%-----------------------------------------------------------
        for j = 1:1:361
          jj = jj+1;
          x1 = r1 * cos(tg(j));
          x2 = r1 * sin(tg(j));
          [cx41, cx42, cx43] = x4fun(x1,x2,r3,l);
          [x41,x42] = qfun(cx41,cx42,cx43);
               x31 = x3fun(x1,x2,x41,r3,l);
               x32 = x3fun(x1,x2,x42,r3,l);
          a(jj,1) = jj;
          a(jj,2) = x1;
          a(jj,3) = x2;
          a(jj,4) = x31;
          a(jj,5) = x41;
          a(jj,6) = x32;
          a(jj,7) = x42;
          a(jj,8) = tang(x31-l,x41);
          a(jj,9) = tang(x32-l,x42);
%------------------------------------------------------------------------------------------------------
         if l < r4+r5 & r4*r5 > 0  % new check
          [cx61,cx62,cx63] = x6fun(x1,x2,x31,x41,r1,r4,r5);
               [x61,x62] = qfun(cx61,cx62,cx63);
               x51 = x5fun(x1,x2,x31,x41,x61,r1,r4,r5);
               x52 = x5fun(x1,x2,x31,x41,x62,r1,r4,r5);
          a(jj,10) = x51;
          a(jj,11) = x61;
          a(jj,12) = x52;
          a(jj,13) = x62;
          [ccx61,ccx62,ccx63] = x6fun(x1,x2,x32,x42,r1,r4,r5);
              [xx61,xx62] = qfun(ccx61,ccx62,ccx63);
                    xx51  = x5fun(x1,x2,x32,x42,xx61,r1,r4,r5);
                    xx52  = x5fun(x1,x2,x32,x42,xx62,r1,r4,r5);
           a(jj,14) = xx51;
           a(jj,15) = xx61;
           a(jj,16) = xx52;
           a(jj,17) = xx62;
       end; % new check
      end;
%-----------------------------------------------------------
if l < r4+r5 & r4*r5 > 0 % new check
  for i1 = 2:1:361
 [a(i1,10),a(i1,11),a(i1,12),a(i1,13)] = organize(a(i1-1,10),a(i1-1,11),a(i1,10),a(i1,11),a(i1,12),a(i1,13));
 [a(i1,14),a(i1,15),a(i1,16),a(i1,17)] = organize(a(i1-1,14),a(i1-1,15),a(i1,14),a(i1,15),a(i1,16),a(i1,17));
end;
end; % new check
%-----------------------------------------------------------
  if l < r4+r5 & r4*r5 > 0 & choice == 1
%   xxx = plot(a(:,2),a(:,3),a(:,4),a(:,5),a(:,10),a(:,11))
    [xxx1,xxx2] = cvelo(a(:,2),a(:,3),a(:,4),a(:,5),a(:,10),a(:,11),rpm,loc);
  end;
%--------------------------------------------------------------------------------
    if l < r4+r5 & r4*r5 > 0 & choice == 2
%      xxx = plot(a(:,2),a(:,3),a(:,4),a(:,5),a(:,12),a(:,13));
     [xxx1,xxx2] =  cvelo(a(:,2),a(:,3),a(:,4),a(:,5),a(:,12),a(:,13),rpm,loc)
    end;
%-------------------------------------------------------------------------------
    if l < r4+r5 & r4*r5 > 0 & choice == 3
%       xxx = plot(a(:,2),a(:,3),a(:,6),a(:,7),a(:,14),a(:,15));
     [xxx1,xxx2] =  cvelo(a(:,2),a(:,3),a(:,6),a(:,7),a(:,14),a(:,15),rpm,loc);
    end;
%----------------------------------------------------------------------------------
     if l < r4+r5 & r4*r5 > 0 & choice == 4
%      xxx = plot(a(:,2),a(:,3),a(:,6),a(:,7),a(:,16),a(:,17));
    [xxx1,xxx2] =   cvelo(a(:,2),a(:,3),a(:,6),a(:,7),a(:,16),a(:,17),rpm,loc);
     end;
%---------------------------------------------------------------------------
  save  link_age.dat  xxx1 -ascii
    disp('The answers are')
    disp('Given crank location')
    disp('Velocity')
    disp('velocity angle')
    disp('Accelleration')
    disp('Accelleration angle')
    disp('.....')
    disp('For three nodes');

author: Sable

repo: mcbench-benchmarks


text:
classdef MSDDetector < handle
    %MSDDETECTOR  Class implementing the MSD (Maximal Self-Dissimilarity) keypoint detector
    %
    % Maximal Self-Dissimilarity Interest Point Detector, as described in
    % [Tombari14].
    %
    % The algorithm implements a novel interest point detector stemming from
    % the intuition that image patches which are highly dissimilar over a
    % relatively large extent of their surroundings hold the property of being
    % repeatable and distinctive. This concept of "contextual
    % self-dissimilarity" reverses the key paradigm of recent successful
    % techniques such as the Local Self-Similarity descriptor and the
    % Non-Local Means filter, which build upon the presence of similar
    % (rather than dissimilar) patches. Moreover, it extends to contextual
    % information the local self-dissimilarity notion embedded in established
    % detectors of corner-like interest points, thereby achieving enhanced
    % repeatability, distinctiveness and localization accuracy.
    %
    % ## References
    % [Tombari14]:
    % > Federico Tombari and Luigi Di Stefano.
    % > "Interest Points via Maximal Self-Dissimilarities".
    % > In Asian Conference on Computer Vision, ACCV 2014.
    %
    % See also: cv.MSDDetector.MSDDetector, cv.FeatureDetector
    %

    properties (SetAccess = private)
        % Object ID
        id
    end

    methods
        function this = MSDDetector(varargin)
            %MSDDETECTOR  The full constructor
            %
            %     obj = cv.MSDDetector()
            %     obj = cv.MSDDetector('OptionName',optionValue, ...)
            %
            % ## Options
            % * __PatchRadius__ Patch radius. default 3
            % * __SearchAreaRadius__ Search Area radius. default 5
            % * __NMSRadius__ Non Maxima Suppression spatial radius. default 5
            % * __NMSScaleRadius__ Non Maxima Suppression scale radius.
            %   default 0
            % * __ThSaliency__ Saliency threshold. default 250.0
            % * __KNN__ Number of nearest neighbors. default 4
            % * __ScaleFactor__ Scale factor for building up the image
            %   pyramid. default 1.25
            % * __NScales__ Number of scales number of scales for building up
            %   the image pyramid (if set to -1, this number is automatically
            %   determined). default -1
            % * __ComputeOrientation__ Flag for associating a canoncial
            %   orientation to each keypoint. default false
            %
            % See also: cv.MSDDetector.detect
            %
            this.id = MSDDetector_(0, 'new', varargin{:});
        end

        function delete(this)
            %DELETE  Destructor
            %
            %     obj.delete()
            %
            % See also: cv.MSDDetector
            %
            if isempty(this.id), return; end
            MSDDetector_(this.id, 'delete');
        end

        function typename = typeid(this)
            %TYPEID  Name of the C++ type (RTTI)
            %
            %     typename = obj.typeid()
            %
            % ## Output
            % * __typename__ Name of C++ type
            %
            typename = MSDDetector_(this.id, 'typeid');
        end
    end

    %% Algorithm
    methods
        function clear(this)
            %CLEAR  Clears the algorithm state
            %
            %     obj.clear()
            %
            % See also: cv.MSDDetector.empty, cv.MSDDetector.load
            %
            MSDDetector_(this.id, 'clear');
        end

        function b = empty(this)
            %EMPTY  Checks if detector object is empty
            %
            %     b = obj.empty()
            %
            % ## Output
            % * __b__ Returns true if the detector object is empty (e.g in the
            %   very beginning or after unsuccessful read).
            %
            % See also: cv.MSDDetector.clear, cv.MSDDetector.load
            %
            b = MSDDetector_(this.id, 'empty');
        end

        function save(this, filename)
            %SAVE  Saves the algorithm parameters to a file
            %
            %     obj.save(filename)
            %
            % ## Input
            % * __filename__ Name of the file to save to.
            %
            % This method stores the algorithm parameters in the specified
            % XML or YAML file.
            %
            % See also: cv.MSDDetector.load
            %
            MSDDetector_(this.id, 'save', filename);
        end

        function load(this, fname_or_str, varargin)
            %LOAD  Loads algorithm from a file or a string
            %
            %     obj.load(fname)
            %     obj.load(str, 'FromString',true)
            %     obj.load(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __fname__ Name of the file to read.
            % * __str__ String containing the serialized model you want to
            %   load.
            %
            % ## Options
            % * __ObjName__ The optional name of the node to read (if empty,
            %   the first top-level node will be used). default empty
            % * __FromString__ Logical flag to indicate whether the input is a
            %   filename or a string containing the serialized model.
            %   default false
            %
            % This method reads algorithm parameters from the specified XML or
            % YAML file (either from disk or serialized string). The previous
            % algorithm state is discarded.
            %
            % See also: cv.MSDDetector.save
            %
            MSDDetector_(this.id, 'load', fname_or_str, varargin{:});
        end

        function name = getDefaultName(this)
            %GETDEFAULTNAME  Returns the algorithm string identifier
            %
            %     name = obj.getDefaultName()
            %
            % ## Output
            % * __name__ This string is used as top level XML/YML node tag
            %   when the object is saved to a file or string.
            %
            % See also: cv.MSDDetector.save, cv.MSDDetector.load
            %
            name = MSDDetector_(this.id, 'getDefaultName');
        end
    end

    %% Features2D: FeatureDetector
    methods
        function keypoints = detect(this, img, varargin)
            %DETECT  Detects keypoints in an image or image set
            %
            %     keypoints = obj.detect(img)
            %     keypoints = obj.detect(imgs)
            %     [...] = obj.detect(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __img__ Image (first variant), 8-bit grayscale or color image.
            % * __imgs__ Image set (second variant), cell array of images.
            %
            % ## Output
            % * __keypoints__ The detected keypoints. In the first variant, a
            %   1-by-N structure array. In the second variant of the method,
            %   `keypoints{i}` is a set of keypoints detected in `imgs{i}`.
            %
            % ## Options
            % * __Mask__ A mask specifying where to look for keypoints
            %   (optional). It must be a logical or 8-bit integer matrix with
            %   non-zero values in the region of interest. In the second
            %   variant, it is a cell-array of masks for each input image,
            %   `masks{i}` is a mask for `imgs{i}`. Not set by default.
            %
            % See also: cv.MSDDetector.MSDDetector
            %
            keypoints = MSDDetector_(this.id, 'detect', img, varargin{:});
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function f = cancelExponents(f)
%CANCELEXPONENTS  Cancel the exponents of a CHEBTECH.
%   F = CANCELEXPONENTS(F) returns the given F without change, since CHEBTECH
%   has no exponents.
%
% See also EXTRACTBOUNDARYROOTS.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function out = ALM_SADAL_smoothed(D, opts)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Test ALM in the paper ''Fast Alternating Linearization Methods for
%       Minimizing the Sum of Two Convex Functions'', Donald Goldfarb,
%       Shiqian Ma and Katya Scheinberg, Tech. Report, Columbia University,
%       2009 - 2010.
%
% Author: Shiqian Ma
% Date  : Apr. 20, 2010
% IEOR, Columbia University, Copyright (2010)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[m,n] = size(D);
mu = opts.mu; sigma = opts.sigma; rho = opts.rho;
Y = zeros(m,n); gradgY = zeros(m,n);
X = D - Y; Dnorm = norm(D,'fro');
Lambda = zeros(m,n); sv = opts.sv;

for itr = 1: opts.maxitr
    if choosvd(n, sv) == 1
        [U, gamma, V] = lansvd(mu*Lambda-Y+D,sv,'L');
    else
        [U, gamma, V] = svd(mu*Lambda-Y+D, 'econ');
    end
    gamma = diag(gamma);
    gamma_new = gamma-mu*gamma./max(gamma,mu+sigma);
    svp = length(find(gamma > mu));
    if svp < sv
        sv = min(svp + 1, n);
    else
        sv = min(svp + round(0.05*n), n);
    end

    Xp = X;
    X = U*diag(gamma_new)*V';
    Lambda = Lambda - (X+Y-D)/mu;
    muY = mu;
    B = Lambda - (X-D)/muY;
    Yp = Y;

    Y = muY*B - muY*min(rho, max(-rho,muY*B/(sigma+muY)));
    Lambda = Lambda - (X+Y-D)/muY;

    StopCrit = norm(D-X-Y,'fro')/Dnorm;
    relX = norm(X-opts.Xs,'fro')/norm(opts.Xs,'fro');
    relY = norm(Y-opts.Ys,'fro')/norm(opts.Ys,'fro');
%     fprintf('iter: %d, mu: %3.2e, rank(X):%d, relX: %3.2e, relY: %3.2e,crit:%3.2e\n', ...
%         itr, mu, length(find(gamma_new>1e-3)), relX, relY, norm(D-X-Y,'fro')/norm(D,'fro'));
    %fprintf('iter: %d, mu: %3.2e, crit:%3.2e\n', itr, mu, norm(D-X-Y,'fro')/norm(D,'fro'));

    mu = max(opts.muf, mu*opts.eta_mu);
    sigma = max(opts.sigmaf, sigma*opts.eta_sigma);

    if StopCrit < opts.epsilon
        out.X = X; out.Y = Y; out.iter = itr; out.relX = relX; out.relY = relY; out.StopCrit = StopCrit;
        return;
    end
end

out.X = X;
out.Y = Y;
out.iter = itr;
out.relX = relX;
out.relY = relY;
out.StopCrit = StopCrit;

repo: lrslibrary


text:
function [neighbors] = BuildNonCon2D(NGauss, tol)

% function [neighbors] = BuildNonCon2D(NGauss, tol)
% Purpose: find element to element connections through non-conforming interfaces
%          (** elements assumed straight sided **)

Globals2D;

% 1. Build Gauss nodes
[gz, gw] = JacobiGQ(0, 0, NGauss-1);

% 1.1 Find location of vertices of boundary faces
vx1 = VX(EToV(:,[1,2,3])); vx2 = VX(EToV(:,[2,3,1]));
vy1 = VY(EToV(:,[1,2,3])); vy2 = VY(EToV(:,[2,3,1]));

idB = find(EToE==((1:K)'*ones(1,Nfaces)));
x1 = vx1(idB)'; y1 = vy1(idB)';
x2 = vx2(idB)'; y2 = vy2(idB)';

% 1.2 Find those element-faces that are on boundary faces
[elmtsB,facesB] = find(EToE==((1:K)'*ones(1,Nfaces)));
Nbc = length(elmtsB);

sk = 1;
% 2.1 For each boundary face
for b1=1:Nbc
  % 2.2 Find element and face of this boundary face
  k1 = elmtsB(b1); f1 = facesB(b1);

  % 2.3 Find end coordinates of b1'th boundary face
  x11 = x1(b1);  y11 = y1(b1);  x12 = x2(b1);  y12 = y2(b1);

  % 2.4 Compute areas, lengths and face coordinates used in intersection
  % tests comparing b1'th boundary face with all boundary faces
  area1 = abs((x12-x11)*(y1-y11) - (y12-y11)*(x1-x11)); %scale
  area2 = abs((x12-x11)*(y2-y11) - (y12-y11)*(x2-x11));
  L   = (x12-x11)^2 + (y12-y11)^2 ;
  r21 = ((2*x1-x11-x12)*(x12-x11) + (2*y1-y11-y12)*(y12-y11))/L;
  r22 = ((2*x2-x11-x12)*(x12-x11) + (2*y2-y11-y12)*(y12-y11))/L;

  % 2.5 Find range of local face coordinate (bracketed between -1 and 1)
  r1 = max(-1,min(r21,r22)); r2 = min(1,max(r21,r22));

  % 2.6 Compute flag for overlap of b1 face with all other boundary faces
  flag = area1+area2+(r1<= -1 & r2<= -1)+(r1>=1 & r2>=1)+(r2-r1<tol);

  % 2.7 Find other faces with partial matches
  matches = setdiff(find(flag < tol),b1);
  Nmatches = length(matches(:));

  if(Nmatches>0)
    % 3.1 Find matches
    r1 = r1(matches); r2 = r2(matches);

    % 3.2 Find end points of boundary-boundary intersections
    xy11 = 0.5*[x11;y11]*(1-r1) +  0.5*[x12;y12]*(1+r1);
    xy12 = 0.5*[x11;y11]*(1-r2) +  0.5*[x12;y12]*(1+r2);

    % 3.3 For each face-face match
    for n=1:Nmatches

      % 3.4 Store which elements intersect
      k2 = elmtsB(matches(n)); f2 = facesB(matches(n));
      neighbors{sk}.elmtM = k1; neighbors{sk}.faceM = f1;
      neighbors{sk}.elmtP = k2; neighbors{sk}.faceP = f2;

      % 3.5 Build physical Gauss nodes on face fragment
      xg = 0.5*(1-gz)*xy11(1,n) + 0.5*(1+gz)*xy12(1,n);
      yg = 0.5*(1-gz)*xy11(2,n) + 0.5*(1+gz)*xy12(2,n);

      % 3.6 Find local coordinates of Gauss nodes
      [rg1,sg1] = FindLocalCoords2D(k1, xg, yg);
      [rg2,sg2] = FindLocalCoords2D(k2, xg, yg);

      % 3.7 Build interpolation matrices for volume nodes ->Gauss nodes
      gVM = InterpMatrix2D(rg1,sg1); neighbors{sk}.gVM  = gVM;
      gVP = InterpMatrix2D(rg2,sg2); neighbors{sk}.gVP  = gVP;

      % 3.8 Find face normal
      neighbors{sk}.nx = nx(1+(f1-1)*Nfp,k1);
      neighbors{sk}.ny = ny(1+(f1-1)*Nfp,k1);

      % 4.0 Build partial face data lift operator

      % 4.1 Compute weights for lifting
      partsJ = sqrt( (xy11(1,n)-xy12(1,n))^2 + (xy11(2,n)-xy12(2,n))^2 )/2;
      dgw = gw*partsJ/J(1,k1);

      % 4.2 Build matrix to lift Gauss data to volume data
      neighbors{sk}.lift = V*V'*(gVM')*diag(dgw);

      sk = sk+1;
    end
  end
end
return

repo: Aero-matlab


text:
function C = top_left_corner(A, r, c)
%TOP_LEFT_CORNER Select the top-left corner of a matrix
%
% C = top_left_corner(A, r, c)
%
% Inputs:
%  A  a rectangular matrix
%  r  the number of rows of the corner we want to select
%  c  the number of columns of the corner we want to select
% Outputs:
%  C  the r-by-c top left corner of A

end

author: odedstein

repo: sgi-introduction-course


text:
function hdr = read_yokogawa_header_new(filename)

% READ_YOKOGAWA_HEADER_NEW reads the header information from continuous,
% epoched or averaged MEG data that has been generated by the Yokogawa
% MEG system and software and allows that data to be used in combination
% with FieldTrip.
%
% Use as
%  [hdr] = read_yokogawa_header_new(filename)
%
% This is a wrapper function around the functions
% getYkgwHdrSystem
% getYkgwHdrChannel
% getYkgwHdrAcqCond
% getYkgwHdrCoregist
% getYkgwHdrDigitize
% getYkgwHdrSource
%
% See also READ_YOKOGAWA_DATA_NEW, READ_YOKOGAWA_EVENT

% **
% Copyright (C) 2005, Robert Oostenveld and 2010, Tilmann Sander-Thoemmes
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% FIXED
%  txt -> m
%  fopen iee-le

if ~ft_hastoolbox('yokogawa_meg_reader')
    ft_error('cannot determine whether Yokogawa toolbox is present');
end

handles = definehandles;

sys_info = getYkgwHdrSystem(filename);
id = sys_info.system_id;
ver = sys_info.version;
rev = sys_info.revision;
sys_name = sys_info.system_name;
model_name = sys_info.model_name;
clear('sys_info'); % remove structure as local variables are collected in the end

channel_info = getYkgwHdrChannel(filename);
channel_count = channel_info.channel_count;

acq_cond = getYkgwHdrAcqCond(filename);
acq_type = acq_cond.acq_type;

% these depend on the data type
sample_rate        = [];
sample_count       = [];
pretrigger_length  = [];
averaged_count     = [];
actual_epoch_count = [];

switch acq_type
  case handles.AcqTypeContinuousRaw
    sample_rate = acq_cond.sample_rate;
    sample_count = acq_cond.sample_count;
    if isempty(sample_rate) || isempty(sample_count)
      ft_error('invalid sample rate or sample count in ', filename);
      return;
    end
    pretrigger_length = 0;
    averaged_count = 1;

  case handles.AcqTypeEvokedAve
    sample_rate = acq_cond.sample_rate;
    sample_count = acq_cond.frame_length;
    pretrigger_length = acq_cond.pretrigger_length;
    averaged_count = acq_cond.average_count;
    if isempty(sample_rate) || isempty(sample_count) || isempty(pretrigger_length) || isempty(averaged_count)
      ft_error('invalid sample rate or sample count or pretrigger length or average count in ', filename);
      return;
    end
    if acq_cond.multi_trigger.enable
      ft_error('multi trigger mode not supported for ', filename);
      return;
    end
  case handles.AcqTypeEvokedRaw
    sample_rate = acq_cond.sample_rate;
    sample_count = acq_cond.frame_length;
    pretrigger_length = acq_cond.pretrigger_length;
    actual_epoch_count = acq_cond.average_count;
    if isempty(sample_rate) || isempty(sample_count) || isempty(pretrigger_length) || isempty(actual_epoch_count)
      ft_error('invalid sample rate or sample count or pretrigger length or epoch count in ', filename);
      return;
    end
    if acq_cond.multi_trigger.enable
      ft_error('multi trigger mode not supported for ', filename);
      return;
    end

  otherwise
    ft_error('unknown data type');
end
clear('acq_cond'); % remove structure as local variables are collected in the end

coregist = getYkgwHdrCoregist(filename);
digitize = getYkgwHdrDigitize(filename);
source = getYkgwHdrSource(filename);

% put all local variables into a structure, this is a bit unusual matlab programming style
tmp = whos;
orig = [];
for i=1:length(tmp)
  if isempty(strmatch(tmp(i).name, {'tmp', 'ans', 'handles'}))
    orig = setfield(orig, tmp(i).name, eval(tmp(i).name));
  end
end

% convert the original header information into something that FieldTrip understands
hdr = [];
hdr.orig         = orig;                % also store the original full header information
hdr.Fs           = orig.sample_rate;    % sampling frequency
hdr.nChans       = orig.channel_count;  % number of channels
hdr.nSamples     = [];                  % number of samples per trial
hdr.nSamplesPre  = [];                  % number of pre-trigger samples in each trial
hdr.nTrials      = [];                  % number of trials

switch orig.acq_type
  case handles.AcqTypeEvokedAve
    hdr.nSamples    = orig.sample_count;
    hdr.nSamplesPre = orig.pretrigger_length;
    hdr.nTrials     = 1;                % only the average, which can be considered as a single trial
  case handles.AcqTypeContinuousRaw
    hdr.nSamples    = orig.sample_count;
    hdr.nSamplesPre = 0;                % there is no fixed relation between triggers and data
    hdr.nTrials     = 1;                % the continuous data can be considered as a single very long trial
  case handles.AcqTypeEvokedRaw
    hdr.nSamples    = orig.sample_count;
    hdr.nSamplesPre = orig.pretrigger_length;
    hdr.nTrials     = orig.actual_epoch_count;
  otherwise
    ft_error('unknown acquisition type');
end

% construct a cell-array with labels of each channel
for i=1:hdr.nChans
% this should be consistent with the predefined list in ft_senslabel,
% with yokogawa2grad_new and with ft_channelselection
  if     hdr.orig.channel_info.channel(i).type == handles.NullChannel
    prefix = '';
  elseif hdr.orig.channel_info.channel(i).type == handles.MagnetoMeter
    prefix = 'M';
  elseif hdr.orig.channel_info.channel(i).type == handles.AxialGradioMeter
    prefix = 'AG';
  elseif hdr.orig.channel_info.channel(i).type == handles.PlannerGradioMeter
    prefix = 'PG';
  elseif hdr.orig.channel_info.channel(i).type == handles.RefferenceMagnetoMeter
    prefix = 'RM';
  elseif hdr.orig.channel_info.channel(i).type == handles.RefferenceAxialGradioMeter
    prefix = 'RAG';
  elseif hdr.orig.channel_info.channel(i).type == handles.RefferencePlannerGradioMeter
    prefix = 'RPG';
  elseif hdr.orig.channel_info.channel(i).type == handles.TriggerChannel
    prefix = 'TRIG';
  elseif hdr.orig.channel_info.channel(i).type == handles.EegChannel
    prefix = 'EEG';
  elseif hdr.orig.channel_info.channel(i).type == handles.EcgChannel
    prefix = 'ECG';
  elseif hdr.orig.channel_info.channel(i).type == handles.EtcChannel
    prefix = 'ETC';
  end
  hdr.label{i} = sprintf('%s%03d', prefix, i);

  % overwrite EEG-channel labels
  if hdr.orig.channel_info.channel(i).type == handles.EegChannel
    if ~isempty(hdr.orig.channel_info.channel(i).data.name)
      hdr.label{i} = hdr.orig.channel_info.channel(i).data.name;
    end
  end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% this defines some usefull constants
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function handles = definehandles
handles.output = [];
handles.sqd_load_flag = false;
handles.mri_load_flag = false;
handles.NullChannel         = 0;
handles.MagnetoMeter        = 1;
handles.AxialGradioMeter    = 2;
handles.PlannerGradioMeter  = 3;
handles.RefferenceChannelMark = hex2dec('0100');
handles.RefferenceMagnetoMeter       = bitor( handles.RefferenceChannelMark, handles.MagnetoMeter );
handles.RefferenceAxialGradioMeter   = bitor( handles.RefferenceChannelMark, handles.AxialGradioMeter );
handles.RefferencePlannerGradioMeter = bitor( handles.RefferenceChannelMark, handles.PlannerGradioMeter );
handles.TriggerChannel      = -1;
handles.EegChannel          = -2;
handles.EcgChannel          = -3;
handles.EtcChannel          = -4;
handles.NonMegChannelNameLength = 32;
handles.DefaultMagnetometerSize       = (4.0/1000.0);       % Square of 4.0mm in length
handles.DefaultAxialGradioMeterSize   = (15.5/1000.0);      % Circle of 15.5mm in diameter
handles.DefaultPlannerGradioMeterSize = (12.0/1000.0);      % Square of 12.0mm in length
handles.AcqTypeContinuousRaw = 1;
handles.AcqTypeEvokedAve     = 2;
handles.AcqTypeEvokedRaw     = 3;
handles.sqd = [];
handles.sqd.selected_start  = [];
handles.sqd.selected_end    = [];
handles.sqd.axialgradiometer_ch_no      = [];
handles.sqd.axialgradiometer_ch_info    = [];
handles.sqd.axialgradiometer_data       = [];
handles.sqd.plannergradiometer_ch_no    = [];
handles.sqd.plannergradiometer_ch_info  = [];
handles.sqd.plannergradiometer_data     = [];
handles.sqd.eegchannel_ch_no   = [];
handles.sqd.eegchannel_data    = [];
handles.sqd.nullchannel_ch_no   = [];
handles.sqd.nullchannel_data    = [];
handles.sqd.selected_time       = [];
handles.sqd.sample_rate         = [];
handles.sqd.sample_count        = [];
handles.sqd.pretrigger_length   = [];
handles.sqd.matching_info   = [];
handles.sqd.source_info     = [];
handles.sqd.mri_info        = [];
handles.mri                 = [];

author: spm


text:
function K = slkernel(varargin)
%SLKERNEL Computes the kernel for samples
%
% $ Syntax $
%   - K = slkernel(X0, kernel_type, ...)
%   - K = slkernel(X0, X, kernel_type, ...)
%
% $ Description $
%   - K = slkernel(X0, kernel_type, ...) Computes the Gram matrix for
%     the samples in matrix X0 using the kernel specified by kernel_type.
%     kernel_type can be name of a built-in kernel type, or the name,
%     handle of an user-specified function. The user-specified function
%     should take in two d x n matrix containing n pairs of vectors, and
%     output a 1 x n vector of their kernel values. The available builtin
%     kernels are given as follows:
%     \*
%     \t   Table 1.  Built-in Kernel Types          \\
%     \h     name     &    description              \\
%           'lin'     &  Linear Kernel:
%                        x1' * x2,
%                        with no parameters         \\
%           'gauss'   &  Gaussian RBF Kernel:
%                        exp(- ||x1 - x2||^2 / (2 * sigma^2)),
%                        with one parameter sigma   \\
%           'poly'    &  Polynomial Kernel:
%                        ((x1' * x2) + a)^k,
%                        with two parameters: k and a \\
%           'sigmoid' &  Sigmoidal Kernel:
%                        tanh(k * (x1' * x2) + a),
%                        with two parameters: k and a \\
%           'invquad' &  Inverse Quadratic Kernel:
%                        1 / sqrt(||x1 - x2||^2 + a),
%                        with one parameter: a        \\
%     \*
%
%   - K = slkernel(X0, X, kernel_type, p1, p2, ...) Computes the empirical
%     kernel mapping for samples X with respect to data set X0. Suppose
%     there are n0 samples in X0, n samples in X, then K will be an n0 * n
%     matrix. Each column in K corresponds to a column in X.
%
% $ History $
%   - Created by Dahua Lin on Jul 13rd, 2005
%   - Modified by Dahua Lin on May 2nd, 2005
%       - base on sltoolbox v4
%       - re-organize the code structure
%       - add the ability of user-specified kernel functions.
%

%% parse and verify input arguments

% for X0
if ~isnumeric(varargin{1})
    error('sltoolbox:invalidarg', ...
        'The first argument should be an numeric matrix');
end
X0 = varargin{1};
if ndims(X0) ~= 2
    error('sltoolbox:invaliddims', ...
        'The X0 should be a 2D matrix');
end

% for X
if isnumeric(varargin{2})
    X = varargin{2};
    if ndims(X) ~= 2
        error('sltoolbox:invaliddims', ...
            'X should be a 2D matrix');
    end
    ipkt = 3;       % argument index of kernel type
else
    X = X0;
    ipkt = 2;
end

% for kernel type
if nargin < ipkt || isempty(varargin{ipkt})
    error('sltoolbox:invalidarg', ...
        'kernel type is not specified');
end
kernel_type = varargin{ipkt};

% for extra parameters
if nargin == ipkt       % no extra parameters
    params = {};
else
    params = varargin(ipkt+1:end);
end

%% Determine the computation routine and delegate to it

% determine for the built-in ones

bik = false;
if ischar(kernel_type)
    switch kernel_type
        case 'lin'
            bik = true;   % bik means built-in kernel
            fh_kernel = @lin_kernel;
        case 'gauss'
            bik = true;
            fh_kernel = @gauss_kernel;
        case 'poly'
            bik = true;
            fh_kernel = @poly_kernel;
        case 'sigmoid'
            bik = true;
            fh_kernel = @sigmoid_kernel;
        case 'invquad'
            bik = true;
            fh_kernel = @invquad_kernel;
    end
end

% delegate
if bik
    K = fh_kernel(X0, X, params{:});
else
    K = slpweval(X0, X, kernel_type, params{:});
end

%% The built-in kernel functions

% Linear Kernel
function K = lin_kernel(X0, X)

K = X0' * X;

% Gaussian Radius Base Function Kernel
function K = gauss_kernel(X0, X, sigma)

D2 = slmetric_pw(X0, X, 'sqdist');
K = exp(- D2 / (2 * sigma * sigma));

% Polynomial Kernel
function K = poly_kernel(X0, X, k, a)

K = (X0' * X + a).^k;

% Sigmoidal Kernel
function K = sigmoid_kernel(X0, X, k, a)

K = tanh(k * (X0' * X) + a);

% Inverse Quadratic Kernel
function K = invquad_kernel(X0, X, a)

D2 = slmetric_pw(X0, X, 'sqdist');
K = 1 ./ sqrt(D2 + a);

author: lmthang


text:
function message = sysrequirements_for_testing(varargin)
% Checking whether the functions needed for testing are installed.
%     In order to run m-functions timing_arraylab_engines and
%     testing_memory_usage, you need to have installed in your system:
%         - MATLAB R2007a or later (you need the builtin function BSXFUN)
%         - GENOP (MATLAB Central file #10333)
%         - BSXFUN substitute (MATLAB Central file #23005)
%         - TIMEIT (MATLAB Central file #18798)

% Paolo de Leva
% University of Rome, Foro Italico, Rome, Italy
% 2009 Feb 22

disp ' '
disp(['You are running MATLAB version ' version])
errorflag = false;

  bsxfunflag = false;
  bsxmexflag = false;
   genopflag = false;
bsxtimesflag = false;
  timeitflag = false;

for i = 1:nargin
    switch varargin{i}
        case 'bsxfun',     bsxfunflag = true;
        case 'bsxmex',     bsxmexflag = true;
        case 'genop',       genopflag = true;
        case 'bsxtimes', bsxtimesflag = true;
        case 'timeit',     timeitflag = true;
    end
end

% Check to see if BSXFUN is a builtin function.
if bsxfunflag && ~exist('bsxfun', 'builtin')
    disp ' '
    w = (['WARNING: You need MATLAB 7.4 (R2007a) or a later version, because  '
          '         you need the builtin function BSXFUN (its replacements for'
          '         releases prior to R2007a are also tested here).           ']);
    disp (w)
    errorflag = true;
end

if genopflag && ~exist('genop', 'file')
    disp ' '
    w = (['WARNING: You need to download GENOP (MATLAB Central file #10333,'
          '         http://www.mathworks.com/matlabcentral/fileexchange)   '
          '         and move it into the current directory.                ']);
    disp (w)
    errorflag = true;
end

if bsxmexflag && (~exist('bsxfun', 'file') || ~exist('bsx_times', 'file'))
    disp ' '
    w = (['WARNING: You need to download Schwarz''s BSXFUN substitute (file #23005'
          '         http://www.mathworks.com/matlabcentral/fileexchange), install'
          '         it, add to the MATLAB search path the directory in which you '
          '         installed it (use PATHTOOL to set the path).                 '
          '         Do not forget to run MAKE_BSX_MEX to build BSX_TIMES.        ']);
      disp (w)
    errorflag = true;
end

if bsxtimesflag && ~exist('bsx_times', 'file')
    disp ' '
    w = (['WARNING: You need to download Schwarz''s BSXFUN substitute (file #23005'
          '         http://www.mathworks.com/matlabcentral/fileexchange), install'
          '         it, delete BSXFUN.M to avoid a conflict with the builtin     '
          '         BSXFUN, run MAKE_BSX_MEX to build BSX_TIMES, and move        '
          '         BSX_TIMES to the current directory.                          ']);
      disp (w)
    errorflag = true;
end

if timeitflag
    verdate = datenum(version('-date'), 'mmmm dd, yyyy');
    minimumdate = datenum(2005,6,21); % Correct version of functions TIC/TOC
    if verdate < minimumdate
        disp ' '
        w = (['WARNING: You need MATLAB 7.1 (R14SP3) or a later version, because '
              '         you need an accurate and precise version of the builtin  '
              '         functions TIC and TOC, supporting the undocumented syntax'
              '         T = TIC; ...; ELAPSED_TIME = TOC(T).                     ']);
        disp (w)
        errorflag = true;
    elseif ~exist('timeit', 'file')
        disp ' '
        w = (['WARNING: You need to download TIMEIT (MATLAB Central file #18798,   '
              '         http://www.mathworks.com/matlabcentral/fileexchange),      '
              '         install it, and add to the MATLAB search path the directory'
              '         in which you installed it (use PATHTOOL to set the path).  ']);
        disp (w)
        errorflag = true;
    end
end

if errorflag
    message = (['Before running this function you must install the above\n', ...
                'mentioned software']);
else
    message = false;
end


text:
% EXP    Exponential.
%    EXP(X) is the exponential of the elements of X, e to the X.
%    For complex Z=X+i*Y, EXP(Z) = EXP(X)*(COS(Y)+i*SIN(Y)).
%
%    See also EXPM1, LOG, LOG10, EXPM, EXPINT.
%
%    Reference page in Doc Center
%       doc exp
%
%    Other functions named exp
%
%       codistributed/exp    gpuArray/exp    sym/exp    ts/exp
%       fints/exp
%

author: jmaih


text:
function value = angle_rad_2d ( p1, p2, p3 )

%*****************************************************************************80
%
%% ANGLE_RAD_2D returns the angle swept out between two rays in 2D.
%
%  Discussion:
%
%    Except for the zero angle case, it should be true that
%
%      ANGLE_RAD_2D(P1,P2,P3) + ANGLE_RAD_2D(P3,P2,P1) = 2 * PI
%
%        P1
%        /
%       /
%      /
%     /
%    P2--------->P3
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real P1(2,1), P2(2,1), P3(2,1), define the rays
%    P1 - P2 and P3 - P2 which in turn define the angle.
%
%    Output, real VALUE, the angle swept out by the rays, measured
%    in radians.  0 <= VALUE < 2*PI.  If either ray has zero length,
%    then VALUE is set to 0.
%
  p(1,1) = ( p3(1,1) - p2(1,1) ) * ( p1(1,1) - p2(1,1) ) ...
         + ( p3(2,1) - p2(2,1) ) * ( p1(2,1) - p2(2,1) );

  p(2,1) = ( p3(1,1) - p2(1,1) ) * ( p1(2,1) - p2(2,1) ) ...
         - ( p3(2,1) - p2(2,1) ) * ( p1(1,1) - p2(1,1) );

  if ( p(1,1) == 0.0 & p(2,1) == 0.0 )
    value = 0.0;
    return
  end

  value = atan2 ( p(2,1), p(1,1) );

  if ( value < 0.0 )
    value = value + 2.0 * pi;
  end

  return
end


text:
function M = Q2M(Q)
% Generate a rotation matrix from a quaternion xi+yj+zk+w,
% where Q = [x y z], and w = 1-x^2-y^2-z^2.
% See: http://skal.planet-d.net/demo/matrixfaq.htm
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

%
% $Id$

Q = Q(1:3); % Assume rigid body
w = sqrt(1 - sum(Q.^2));
x = Q(1); y = Q(2); z = Q(3);
if w<1e-7,
    w = 1/sqrt(x*x+y*y+z*z);
    x = x*w;
    y = y*w;
    z = z*w;
    w = 0;
end;
xx = x*x; yy = y*y; zz = z*z; ww = w*w;
xy = x*y; xz = x*z; xw = x*w;
yz = y*z; yw = y*w; zw = z*w;
M = [...
(xx-yy-zz+ww)      2*(xy-zw)      2*(xz+yw) 0
    2*(xy+zw) (-xx+yy-zz+ww)      2*(yz-xw) 0
    2*(xz-yw)      2*(yz+xw) (-xx-yy+zz+ww) 0
           0              0              0  1];
return;


text:
% This example loads the 'nps plus' image set of patterns, and examines the
% relationships between these patterns and the standard Buckner Lab
% networks.

% Layer 1 of the river plot is the Wager lab patterns, and Layer 2 is the
% Buckner Lab networks

%% Load NPSplus

[npsplus_obj, layer1names, imgnames] = load_image_set('npsplus');

wh = 1:4;
npsplus_obj = get_wh_image(npsplus_obj, wh);
layer1names = layer1names(wh);

layer1names

%% Get rectangular similarity matrix

stats = image_similarity_plot(npsplus_obj, 'bucknerlab', 'noplot', 'cosine_similarity');
% Similarity matrix is stats.r

% Threshold to look at positive associations
sim_matrix = stats.r;
sim_matrix(sim_matrix < 0) = 0;

layer2names = stats.networknames;

[n2, n1] = size(sim_matrix);

%layer1colors = scn_standard_colors(n1);
%layer1colors = seaborn_colors(n1);

layer1colors = {[1 .4 .2] [.3 .3 1] [.7 .3 1] [.3 .7 .5]};

layer2colors = seaborn_colors(n2);

%% Reorder rows and columns for display, as needed

% Custom REORDERING of layer 1 for aesthetics
neworder = [1 3 2 4];
sim_matrix = sim_matrix(:, neworder);
layer1names = layer1names(neworder);
layer1colors = layer1colors(neworder);

%  REORDER rows here based on similarity to make plot look cleaner
[sim_matrix, layer2names, layer2colors] = riverplot_reorder_matrix(sim_matrix, layer2names, layer2colors);

% Custom REORDERING of layer 2 because the automatic one isn't super
neworder = [1 7 2 3 4 5 6];
sim_matrix = sim_matrix(neworder, :);
layer2names = layer2names(neworder);
layer2colors = layer2colors(neworder);

%% Create the plot

create_figure('riverplot');
set(gca, 'YDir', 'reverse');
layer1 = riverplot_draw_layer(2, n1, 'colors', layer1colors, 'y_loc', 1.5);
layer2 = riverplot_draw_layer(5, n2, 'colors', layer2colors);

ribbons = riverplot_ribbon_matrix(layer1, layer2, sim_matrix, 'colors', layer1colors, 'layer1fullcoverage', 'steepness', 0);

% Turn off lines
riverplot_toggle_lines(ribbons);

% Increase opacity
riverplot_set_ribbon_property(ribbons, 'FaceAlpha', .6);

set(gca, 'XLim', [0 8]);

% Add names
layer1 = riverplot_layer_names(layer1, layer1names);
layer2 = riverplot_layer_names(layer2, layer2names, 'right');

axis off

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [dat] = read_neuralynx_dma(filename, begsample, endsample, channel)

% READ_NEURALYNX_DMA reads specified samples and channels data from a Neuralynx DMA log file
%
% Use as
%    [hdr] = read_neuralynx_dma(filename)
%    [dat] = read_neuralynx_dma(filename, begsample, endsample)
%    [dat] = read_neuralynx_dma(filename, begsample, endsample, chanindx)
%
% The channel specification can be a vector with indices, or a single string with the value
%    'all', 'stx', 'pid', 'siz', 'tsh', 'tsl',
%    'cpu', 'ttl', 'x01',  ...,  'x10'
%
% This function returns the electrophysiological data in AD units
% and not in uV. You should look up the details of the headstage and
% the Neuralynx amplifier and scale the values accordingly.

% Copyright (C) 2005-2008, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The data is simply a stream of constant size records,  the size of the
% record is dependant on the number of Input Boards in the system.  8 in your
% case, there are 32 - 32 bit integers from each Input Board.
%
% // Digital Lynx Packets look like this for a 1 board system:
% // stx (or SOP)    0x00000800
% // Packet ID       0x00000001
% // siz             0x0000002A   (# A/D data ints + #extra wds = 32+10)
% // TimeStamp Hi    1 32 bit word
% // TimeStamp Low   1 32 bit word
% // cpu             1 32 bit word
% // ttl             1 32 bit word
% // 10 extras      10 32 bit words
% // A/D data       32 32 bit words per board
% // CRC             1 32 bit XOR of the entire packet including stx
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

needhdr = (nargin==1);
needdat = (nargin>=2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% since Cheetah 4.80 the DMS log files have a 16384 byte header
% determine whether that is the case in this file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fid = fopen(filename, 'rb', 'ieee-le');
buf = fread(fid, 10, 'uchar=>uchar');
if all(buf(1:4)==35)
  % the file starts with "####"
  hdroffset = 16384;
elseif all(buf>31 & buf<127)
  % the file does not start with "####" but the header seems to be ascii anyway
  hdroffset = 16384;
else
  hdroffset = 0;
end

if hdroffset
  % read the ascii header from the file, it does not contain much usefull information
  orig = neuralynx_getheader(filename);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% since Cheetah 5.0 the DMA log file does not always start with a complete
% packet but can start with an incomplete "junk" packet after the header.
% The junk due to that incomplete packet should then also be skipped.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if hdroffset
  % skip the standard ascii header
  fseek(fid, hdroffset, 'bof');
  % read a single block
  buf = fread(fid, 274, 'uint32=>uint32');
  % search for the first known values of a proper packet
  junk = find((buf(1:(end-1))==2048) & (buf(2:(end-0))==1), 1) - 1;
  if ~isempty(junk)
    % add the additional junk samples to the header offset
    hdroffset  = hdroffset + junk*4;
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% since modifying our hardware from 256 to 32 channel (19 March 2007), the
% number of channels in a DMA log file can vary. Start with the assumption
% that the file is from the 256 channel system with 8 boards that we have
% in Nijmegen and verify the CRC
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nboards    = 8;
nchanboard = 32;
blocksize  = nboards*nchanboard+18; % in 4 byte words
fseek(fid, hdroffset, 'bof');
% read a single block
buf = fread(fid, blocksize, 'uint32=>uint32');
% determine the number of boards by systematically checking the CRC value
while (nboards>0)
  if neuralynx_crc(buf(1:blocksize-1))==buf(blocksize)
    % the CRC is correct for this number of boards
    break
  else
    nboards    = nboards - 1;           % decrease by one
    blocksize  = nboards*nchanboard+18; % in 4 byte words
  end
end
fclose(fid);

if ~nboards
  ft_error('could not determine the number of boards and channels');
end

% deal with ascii and numeric input for the channel selection
if nargin>3 && ~ischar(channel)
  chanindx = channel;
  channel = 'selection';
elseif nargin>3 && ischar(channel)
  switch channel
    case 'stx'
      chanindx = 1;
    case 'pid'
      chanindx = 2;
    case 'siz'
      chanindx = 3;
    case 'tsh'
      chanindx = 4;
    case 'tsl'
      chanindx = 5;
    case 'cpu'
      chanindx = 6;
    case 'ttl'
      chanindx = 7;
    case 'x01'
      chanindx = 8;
    case 'x02'
      chanindx = 9;
    case 'x03'
      chanindx = 10;
    case 'x04'
      chanindx = 11;
    case 'x05'
      chanindx = 12;
    case 'x06'
      chanindx = 13;
    case 'x07'
      chanindx = 14;
    case 'x08'
      chanindx = 15;
    case 'x09'
      chanindx = 16;
    case 'x10'
      chanindx = 17;
    case 'crc'
      chanindx = blocksize;
    case 'all'
      chanindx = 1:blocksize;
    otherwise
      ft_error('unknown value in channel');
  end
elseif nargin<4
  channel = 'all';
  chanindx = 1:blocksize;
end

if needhdr
  % these have to be hardcoded, since the DMA logfile does not contain header information
  hdr              = [];
  if hdroffset
    % remember the ascii header from the file, it does not contain much usefull information
    hdr.orig       = orig;
  end
  hdr.Fs           = 32556;                   % sampling frequency
  hdr.nChans       = nboards*nchanboard;      % number of analog channels
  hdr.nSamples     = inf;                     % number of samples per trial
  hdr.nSamplesPre  = 0;                       % number of pre-trigger samples in each trial
  hdr.nTrials      = 1;                       % number of trials
  hdr.label        = {};                      % cell-array with labels of each channel
  for i=1:hdr.nChans
    chanlabel{i} = sprintf('csc%03d', i);
  end
  statuslabel = {
    'stx'
    'pid'
    'siz'
    'tsh'
    'tsl'
    'cpu'
    'ttl'
    'x01'
    'x02'
    'x03'
    'x04'
    'x05'
    'x06'
    'x07'
    'x08'
    'x09'
    'x10'
    };
  hdr.label       = cat(1, statuslabel(:), chanlabel(:), {'crc'});  % concatenate all channel labels
  hdr.nChans      = length(hdr.label);                              % this includes all channels
  % determine the length of the file, expressed in samples
  fid = fopen(filename, 'rb', 'ieee-le');
  fseek(fid, 0, 'eof');
  hdr.nSamples = floor((ftell(fid)-hdroffset)/(blocksize*4));

  % determine the timestamp of the first and of the last sample
  datoffset = 0*blocksize;
  fseek(fid, hdroffset + datoffset, 'bof');
  buf = fread(fid, blocksize, 'uint32=>uint32');
  beg_stx = buf(1);
  beg_tsh = buf(4);
  beg_tsl = buf(5);
  datoffset = (hdr.nSamples-1)*blocksize*4;
  fseek(fid, hdroffset + datoffset, 'bof');
  buf = fread(fid, blocksize, 'uint32=>uint32');
  end_stx = buf(1);
  end_tsh = buf(4);
  end_tsl = buf(5);
  fclose(fid);
  % check that the junk at the beginning was correctly detected
  if (beg_stx~=2048)
    ft_error('problem with STX at the begin of the file');
  end
  % check that the file is truely continuous, i.e. no gaps with junk in between
  if (end_stx~=2048)
    ft_error('problem with STX at the end of the file');
  end
  % combine the two uint32 words into a uint64
  hdr.FirstTimeStamp = timestamp_neuralynx(beg_tsl, beg_tsh);
  hdr.LastTimeStamp  = timestamp_neuralynx(end_tsl, end_tsh);
  hdr.TimeStampPerSample = double(hdr.LastTimeStamp-hdr.FirstTimeStamp)./(hdr.nSamples-1);  % this should be double, since it can be fractional

  % only return the header information
  dat = hdr;
  dat.orig.Offset = hdroffset;

elseif length(chanindx)>1
  % read the desired samples for all channels
  fid = fopen(filename, 'rb', 'ieee-le');
  datoffset = (begsample-1)*blocksize*4;
  fseek(fid, hdroffset + datoffset, 'bof');  % skip the header and samples that do not have to be read
  % read the data
  dat = fread(fid, [blocksize (endsample-begsample+1)], 'int32=>int32');
  fclose(fid);
  if size(dat,2)<(endsample-begsample+1)
    ft_error('could not read all samples');
  end
  if ~strcmp(channel, 'all')
    % select the subset of desired channels
    dat = dat(chanindx,:);
  end

elseif length(chanindx)==1
  % read the desired samples for only one channel
  fid = fopen(filename, 'rb', 'ieee-le');
  datoffset = (begsample-1)*blocksize*4;
  fseek(fid, hdroffset + datoffset, 'bof');  % skip the header and samples that do not have to be read
  fseek(fid, (chanindx-1)*4, 'cof');         % skip to the channel of interest
  % Note that the last block with 274*4 bytes can sometimes be incomplete, which is relevant when endsample=inf
  dat = fread(fid, [1 (endsample-begsample+1)], 'int32=>int32', (nboards*32+18-1)*4);
  if size(dat,2)<(endsample-begsample+1)
    ft_error('could not read all samples');
  end
  fclose(fid);
end


text:
function output = base64encode(input)
%BASE64ENCODE Encode a byte array using Base64 codec.
%
%    output = base64encode(input)
%
% The function takes a char, int8, or uint8 array INPUT and returns Base64
% encoded string OUTPUT. JAVA must be running to use this function. Note
% that encoding doesn't preserve input dimensions.
%
% See also base64decode

error(nargchk(1, 1, nargin));
error(javachk('jvm'));
if ischar(input), input = uint8(input); end

output = char(org.apache.commons.codec.binary.Base64.encodeBase64Chunked(input))';

end