text:
function newnums=roundtowardvec(X,roundvec,type)
%function newnums=roundtowardvec(X,[roundvec],[type])
%
% This function rounds number(s) toward given values. If more than one
% number is given to round, it will return the matrix with each rounded
% value, otherwise it will return the single rounded value. It will ignore
% NaNs and return them back with NaNs.
%
% Inputs: X: the number(s) that you want rounded
%
%         roundvec:(opt) the values to round X to. If none given, it will
%           default to -inf:1:inf (and use the built in functions).
%
%         type:(opt) specifies which kind of rounding you want
%           the function to use.
%
%           Choices are: 'round' - round to nearest value
%                        'floor' - round toward -Inf
%                        'ceil'  - round toward Inf
%                        'fix'   - round toward 0
%                        'away'  - round away from 0 (ceil if positive and floor if negative)
%                     (see help files for more clarity)
%
%           If no type is given, the function will default to rounding to
%           the nearest value.
%
% Outputs: newnums: rounded values, in same shape as X input matrix

% For some reason, loops seem to be faster than vectorizing this code

if nargin==0
    help roundtowardvec; %if nothing given, tell what to give
    return
elseif isempty(X)
    newnums=[]; %if given empty, return empty without going through whole script
    return
end
if ~exist('type','var') || isempty(type)
    type='round';  %%round to nearest value if not specified
end
if ~exist('roundvec','var') || isempty(roundvec)
    if strcmpi(type,'round')
        newnums=round(X);
        %to nearest integer
    elseif strcmpi(type,'away')
        newnums=ceil(abs(X)).*sign(X);
        %nearest integer away from 0
    elseif strcmpi(type,'fix')
        newnums=fix(X);
        %nearest integer toward 0
    elseif strcmpi(type,'floor')
        newnums=floor(X);
        %nearest integer toward -inf
    elseif strcmpi(type,'ceil')
        newnums=ceil(X);
        %nearest integer toward inf
    else
        error(sprintf('Round type not recognized. Options are:\n''round'' - round to nearest value\n''floor'' - round toward -Inf\n''ceil''  - round toward Inf\n''fix''   - round toward 0\n''away''  - round away from 0')) %#ok<SPERR>
    end
else
    %%make newnums size of X
    newnums=X;
    if strcmpi(type,'round') %to nearest value
        roundvec=reshape(unique(roundvec),1,[]);
        for i=numel(X):-1:1
            if ~any(X(i)==roundvec)
                DIFFs=abs(roundvec-X(i));
                if X(i)>=0
                    [~,ind]=min(DIFFs(:,end:-1:1));
                    newnums(i)=roundvec(length(DIFFs)-ind+1);
                elseif X(i)<0
                    [~,ind]=min(DIFFs);
                    newnums(i)=roundvec(ind);
                end
            end
        end
    elseif strcmpi(type,'fix') %to nearest value toward 0
        roundvec=reshape(unique([roundvec 0]),1,[]);
        for i=numel(X):-1:1
            if ~any(X(i)==roundvec)
                if X(i)>0
                    if X(i)>min(roundvec)
                        newnums(i)=roundvec(find(X(i)>roundvec,1,'last'));
                    else
                        newnums(i)=0;
                    end
                elseif X(i)<0
                    if X(i)<max(roundvec)
                        newnums(i)=roundvec(find(X(i)<roundvec,1,'first'));
                    else
                        newnums(i)=0;
                    end
                end
            end
        end
    elseif strcmpi(type,'ceil') %nearest value toward inf
        roundvec=reshape(unique(roundvec),1,[]);
        for i=numel(X):-1:1
            if ~isnan(X(i)) && ~any(X(i)==roundvec)
                if X(i)<max(roundvec)
                    newnums(i)=roundvec(find(X(i)<roundvec,1,'first'));
                else
                    newnums(i)=inf;
                end
            end
        end
    elseif strcmpi(type,'floor') %nearest value toward -inf
        roundvec=reshape(unique(roundvec),1,[]);
        for i=numel(X):-1:1
            if ~isnan(X(i)) && ~any(X(i)==roundvec)
                if X(i)>min(roundvec)
                    newnums(i)=roundvec(find(X(i)>roundvec,1,'last'));
                else
                    newnums(i)=-inf;
                end
            end
        end
    elseif strcmpi(type,'away') %nearest value away from 0
        roundvec=reshape(unique(roundvec),1,[]);
        for i=numel(X):-1:1
            if ~any(X(i)==roundvec)
                if X(i)>0
                    if X(i)<max(roundvec)
                        newnums(i)=roundvec(find(X(i)<roundvec,1,'first'));
                    else
                        newnums(i)=inf;
                    end
                elseif X(i)<0
                    if X(i)>min(roundvec)
                        newnums(i)=roundvec(find(X(i)>roundvec,1,'last'));
                    else
                        newnums(i)=-inf;
                    end
                elseif X(i)==0
                    DIFFs=abs(roundvec-X(i));
                    [~,ind]=min(DIFFs(:,end:-1:1));
                    newnums(i)=roundvec(length(DIFFs)-ind+1);
                end
            end
        end
    else
        error(sprintf('Round type not recognized. Options are:\n''round'' - round to nearest value\n''floor'' - round toward -Inf\n''ceil''  - round toward Inf\n''fix''   - round toward 0\n''away''  - round away from 0')) %#ok<SPERR>
    end
end
end

author: Sable

repo: mcbench-benchmarks


text:
function value = p27_r8vec ( action, name, dim_num, value )

%*****************************************************************************80
%
%% P27_R8VEC sets or gets R8VEC parameters for problem 27.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 November 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string ACTION, the action.
%    'D' sets the internal value of the object to a default value.
%    If NAME = '*', then all variables are defaulted.
%    'G' means the current value of the object should be returned.
%    'R' means randomize the object and return it.
%    'S' means the input values of the object and its dimension should
%    be stored.
%
%    Input, string NAME, the name of the parameter.
%    'C' is the coefficient vector.
%
%    Input, integer DIM_NUM, the dimension of the object.
%
%    Input/output, real VALUE(DIM_NUM,1), the value of the object.
%
  persistent c
  persistent dim_num_save

  if ( size ( dim_num_save ) == 0 )
    dim_num_save = 0;
  end

  if ( dim_num_save ~= dim_num )
    dim_num_save = 0;
    c = [];
  end

  if ( dim_num_save == 0 )
    dim_num_save = dim_num;
  end

  if ( action == 'D' | action == 'd' )

    if ( name == 'C' | name == 'c' | name == '*' )
      c(1:dim_num,1) = 1.0 / dim_num;
    end

  elseif ( action == 'G' | action == 'g' )

    if ( name == 'C' | name == 'c' )
      value = zeros ( 1:dim_num, 1 );
      value(1:dim_num,1) = c(1:dim_num,1);
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P27_R8VEC - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized name = "%s".\n', name );
      error ( 'P27_R8VEC - Fatal error!' );
    end

  elseif ( action == 'R' | action == 'r' )

    if ( name == 'C' | name == 'c' )
      c = rand ( 1, dim_num );
      c(1:dim_num,1) = c(1:dim_num,1) * 25.0 / dim_num;
      value(1:dim_num,1) = c(1:dim_num,1);
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P27_R8VEC - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized name = "%s".\n', name );
      error ( 'P27_R8VEC - Fatal error!' );
    end

  elseif ( action == 'S' | action == 's' )

    if ( name == 'C' | name == 'c' )
      c(1:dim_num,1) = value(1:dim_num,1);
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P27_R8VEC - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized name = "%s".\n', name );
      error ( 'P27_R8VEC - Fatal error!' );
    end

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'P27_R8VEC - Fatal error!\n' );
    fprintf ( 1, '  Unrecognized action = "%s".\n', action );
    error ( 'P27_R8VEC - Fatal error!' );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function pde = Stokesdata1
%% STOKESDATA1 data for Stokes equations
%
% A simple model of colliding flow. The force f = 0, the velocity  u1 =
% 20xy^3, u2 = 5x^4-5y^4, p = 60x^2y - 20y^3.
%
% Dirichlet boundary condition is imposed.
%
% Reference: page 237 in Finite Elements and Fast Iterative Solvers with
% Applications in Incompressible Fluid Dynamics. by Howard C. Elman, David
% J. Silvester, and Andrew J. Wathen.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

pde = struct('f', 0, 'exactp', @exactp, 'exactu', @exactu,'g_D',@exactu, 'exactw', @exactw);

    % exact velocity
    function z = exactu(p)
        x = p(:,1); y = p(:,2);
        z(:,1) = 20*x.*y.^3;
        z(:,2) = 5*x.^4-5*y.^4;
    end
    % exact pressure
    function z = exactp(p)
        x = p(:,1); y = p(:,2);
        z = 60*x.^2.*y-20*y.^3-5;
    end
    % exact vorticity
    function z = exactw(p)
        x = p(:,1); y = p(:,2);
        z = 20*x.^3-60*x.*y.^2;
    end
end

repo: ifem


text:
function substituteDefaultPropNamesWithCustomNamesInDepVarListbox(hDepVarListbox, propNames)
    vars = get(hDepVarListbox,'String');
%     propNames = maData.spacecraft.propellant.names;
    defaultNames = {'Liquid Fuel/Ox Mass','Monopropellant Mass','Xenon Mass'};
    propNames = sort(propNames);

    for(i=1:length(vars)) %#ok<*NO4LP>
        for(j=1:length(defaultNames))
            vars{i} = strrep(vars{i}, defaultNames{j}, [propNames{j}, ' Mass']);
        end
    end
    set(hDepVarListbox,'String',vars);
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%
%  Convert a coefficient matrix to a (sparse) coefficient vector
%
%    Syntax:  >> v = mvPolynCoefMat2Vec(F)
%             >> v = mvPolynCoefMat2Vec(F,d)
%
%    Input:      F -- multivariate polynomials in coeff. matrix
%                d -- (optional) degree bound for the polynomial vector space
%                     if not provided, the tuple degree of F is used
%
%   Output:      v -- (sparse vector) the coefficient vector
%
%   Example:
%      >> F = PolynString2CoefMat('3 + 5.5*x^2 - 7*x*y^3',{'x','y'})
%
%      F =
%
%               0    2.0000    1.0000
%               0         0    3.0000
%          3.0000    5.5000   -7.0000
%
%      >> f = mvPolynCoefMat2Vec(F)
%
%      ans =
%
%         (1,1)       3.0000
%         (3,1)       5.5000
%        (11,1)      -7.0000
%

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [A, b, Aeq, beq, lb, ub] = l_constraints(mpcModel, k, varargin)

        A   = []; %A*u<=b
        b   = [];
        Aeq = [1 1];
        beq = mpcModel.net_load(k);

        ub  = [10, mpcModel.battery.power(1,1)];
        lb  = [-5, mpcModel.battery.power(1,2)];

end

author: juchengquan


text:
function model = lleCreate(inputDim, outputDim, Y, options)

% LLECREATE Locally linear embedding model.
% FORMAT
% DESC creates a structure for a locally linear embedding.
% ARG latentDimension : dimension of latent space.
% ARG outputDim : dimension of data.
% ARG Y : the data to be modelled in design matrix format (as many
% rows as there are data points).
% ARG options : options structure as returned by lleOptions.
% RETURN model : model structure containing LLE model.
%
% COPYRIGHT : Neil D. Lawrence, 2008, 2009
%
% SEEALSO : lleOptions, modelCreate

% MLTOOLS

model.type = 'lle';

if size(Y, 2) ~= outputDim
  error(['Input matrix Y does not have dimension ' num2str(d)]);
end
model.isNormalised = options.isNormalised;
model.regulariser = options.regulariser;
model.acyclic = options.acyclic;
model.k = options.numNeighbours;
model.Y = Y;
model.d = outputDim;
model.q = inputDim;
model.N = size(Y, 1);

if isfield(model, 'acyclic') && model.acyclic
  model.indices = findAcyclicNeighbours(model.Y, model.k);
else
  model.indices = findNeighbours(model.Y, model.k);
end
model.W = spalloc(model.N, model.N, model.N*model.k);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function out=imNormalize99(im)

im=double(im);
temp=sort(im(:),'descend');
th1=temp(round(length(im(:))/100));
th2=min(im(:));

out=(im-th2)/(th1-th2);
out(out>1)=1;

author: xiuyechen

repo: FishExplorer


text:
function [min_v,min_i,min_diff] = local_min(F,S,ep)
  % LOCAL_MIN  find values and indices of local minima of a scalar field
  % defined on a mesh
  %
  % [min_v,min_i,min_diff] = local_min(F,S,ep)
  %
  % Inputs:
  %   F #faces by 3 list of triangle indices
  %   S #vertices by 1 list of scalar values
  % Outputs:
  %   min_v  #minima list of values at local minima
  %   min_i  #minima list of indices into S of local minima
  %   min_diff #minima list of min difference between each minima and its
  %     neighbors
  %

  if ~exist('ep','var')
    ep = 0;
  end

  % number of vertices
  n = numel(S);
  if size(F,2) == 2
    E = F;
  else
    E = edges(F);
  end
  E = [E;fliplr(E)];
  A = sparse( ...
    E(:,1), ...
    E(:,2), ...
    S(E(:,1))-S(E(:,2)), ...
    n,n);
  minA = minnz(A)';
  [min_i,~] = find(minA>-ep);
  min_diff = minA(min_i);
  min_v = S(min_i);

end

author: alecjacobson

repo: gptoolbox


text:
%% --------------------------
% DRRN_B1U9
% -------------------------------
function test_DRRN_B1U9()
setenv('LC_ALL','C')
addpath /data2/taiying/MSU_Code/my_caffe/matlab; % change to your caffe path
setenv('GLOG_minloglevel','2')
addpath('../');
addpath('../evaluation_func/');
addpath('../evaluation_func/matlabPyrTools-master/');

%% parameters
gpu_id = 2;
up_scale = 3;
data_set_id = 1;

iter_interval = 6536;
startNUM = 71;
endNUM = 71;
thresh = 1;
thresh_hei = 150; % threshold patch size for inference, since too big image may cost too much memory
thresh_wid = 150;
rf = 16;

pathfolder = ['../../data/'];
if data_set_id == 1
    % Set5
    setTestCur = 'Set5';
    path = [pathfolder setTestCur '/'];
    d = dir([path '*.bmp']);
    filenum = 5;
end
if data_set_id == 2
    % Set14
    setTestCur = 'Set14';
    path = [pathfolder setTestCur '/'];
    d = dir([path '*.bmp']);
    filenum = 14;
end
if data_set_id == 3
    % B100
    setTestCur = 'B100';
    path = [pathfolder setTestCur '/'];
    d = dir([path '*.jpg']);
    filenum = 100;
end
if data_set_id == 4
    % Urban100
    setTestCur = 'Urban100';
    path = [pathfolder setTestCur '/'];
    d = dir([path '*.png']);
    filenum = 100;
end

savepath = ['./results/'];
folderResultCur = fullfile(savepath, [setTestCur,'_x',num2str(up_scale)]);
%%% folder to store results
if ~exist(folderResultCur,'file')
    mkdir(folderResultCur);
end

mean_bicubic = [];
mean_drrn = [];
% caffe.set_mode_cpu(); % for CPU
caffe.set_mode_gpu(); % for GPU
caffe.set_device(gpu_id);

for NUM = startNUM:thresh:endNUM
    epochNUM = NUM*iter_interval;
    weights = ['../../model/DRRN_B1U9_20C128_iter_' num2str(epochNUM) '.caffemodel'];
    model_path = './DRRN_B1U9_20C128_deploy';

    bicubic_set =[];
    drrn_set = [];
    im_h_set = cell(filenum,1);
    im_gnd_set = cell(filenum,1);
    im_other_set = cell(6,1);
    scaleset = [];
    gatingset = [];

    for iii = 1:1:length(d)
        disp(['NUM: ' num2str(NUM) '  id: ' num2str(iii)]);
        imageName = d(iii).name;
        imageName = imageName(1:end-4);
        im  = imread([path d(iii).name]);

        %% work on luminance only
        im_ycbcr= im;
        if size(im,3)>1
            im_ycbcr = rgb2ycbcr(im);
            im_cb = im2double(im_ycbcr(:, :, 2));
            im_cr = im2double(im_ycbcr(:, :, 3));
            im_cb_gnd = modcrop(im_cb, up_scale);
            im_cr_gnd = modcrop(im_cr, up_scale);
            im_cb_b = imresize(imresize(im_cb_gnd,1/up_scale,'bicubic'),up_scale,'bicubic');
            im_cb_b = single(im_cb_b);
            im_cr_b = imresize(imresize(im_cr_gnd,1/up_scale,'bicubic'),up_scale,'bicubic');
            im_cr_b = single(im_cr_b);
        end
        im_y = im2double(im_ycbcr(:, :, 1));
        im_y_gnd = modcrop(im_y, up_scale);

        [hei,wid] = size(im_y_gnd);
        im_y_b = imresize(imresize(im_y_gnd,1/up_scale,'bicubic'),up_scale,'bicubic');
        im_y_b = single(im_y_b);

        %% adaptively spilt
        % decide patch numbers
        hei_patch = ceil(hei/(thresh_hei+rf));
        wid_patch = ceil(wid/(thresh_wid+rf));
        hei_stride = ceil(hei/hei_patch);
        wid_stride = ceil(wid/wid_patch);
        use_start_x = 0;
        use_start_y = 0;
        use_end_x = 0;
        use_end_y = 0;

        ext_start_x = 0;
        ext_end_x = 0;
        ext_start_y = 0;
        ext_end_y = 0;

        posext_start_x = 0;
        posext_start_y = 0;
        posext_end_x = 0;
        posext_end_y = 0;

        % extract each patch for inference
        im_y_h = [];
        temp1 = [];
        temp2 = [];
        temp3 = [];
        temp4 = [];
        temp5 = [];
        temp6 = [];
        for x = 1 : hei_stride : hei
            for y = 1 : wid_stride : wid
                % decide the length of hei and wid for each patch
                use_start_x = x;
                use_start_y = y;
                if x - rf > 1 % add border
                    ext_start_x = x-rf;
                    posext_start_x = rf+1;
                else
                    ext_start_x = x;
                    posext_start_x = 1;
                end
                if y-rf > 1
                    ext_start_y = y-rf;
                    posext_start_y = rf+1;
                else
                    ext_start_y = y;
                    posext_start_y = 1;
                end

                use_end_x = use_start_x+hei_stride-1;
                use_end_y = use_start_y+wid_stride-1;

                if use_start_x+hei_stride+rf-1 <= hei
                    hei_length = hei_stride+rf;
                    ext_end_x = use_start_x+hei_length-1;
                    posext_end_x = hei_length-rf+posext_start_x-1;

                else
                    hei_length = hei-ext_start_x+1;
                    ext_end_x = ext_start_x+hei_length-1;
                    posext_end_x = hei_length;
                    use_end_x = ext_start_x+hei_length-1;
                end
                if use_start_y+wid_stride+rf-1 <= wid
                    wid_length = wid_stride+rf;
                    ext_end_y = use_start_y+wid_length-1;
                    posext_end_y = wid_length-rf+posext_start_y-1;

                else
                    wid_length = wid-ext_start_y+1;
                    ext_end_y = ext_start_y+wid_length-1;
                    posext_end_y = wid_length;
                    use_end_y = ext_start_y+wid_length-1;
                end

                subim_input = im_y_b(ext_start_x : ext_end_x, ext_start_y : ext_end_y);  % input
                data = permute(subim_input,[2, 1, 3]);
                model = [model_path '.prototxt'];
                subim_output = do_cnn(model,weights,data);
                subim_output = subim_output';
                subim_output = subim_output(posext_start_x:posext_end_x,posext_start_y:posext_end_y);

                % fill im_h with sub_output
                im_y_h(use_start_x:use_end_x,use_start_y:use_end_y) = subim_output;

            end
        end

        %% remove border
        im_y_h1 = shave(uint8(single(im_y_h) * 255), [up_scale, up_scale]);
        im_y_gnd1 = shave(uint8(single(im_y_gnd) * 255), [up_scale, up_scale]);
        im_y_b1 = shave(uint8(single(im_y_b) * 255), [up_scale, up_scale]);

        if size(im,3) > 1
            im_cb_b1 = shave(uint8(single(im_cb_b) * 255), [up_scale, up_scale]);
            im_cr_b1 = shave(uint8(single(im_cr_b) * 255), [up_scale, up_scale]);
            ycbcr_h = cat(3,(im_y_h1),(im_cb_b1),(im_cr_b1));
            im_h1 = ycbcr2rgb(ycbcr_h);

            im_cb_gnd1 = shave(uint8(single(im_cb_gnd) * 255), [up_scale, up_scale]);
            im_cr_gnd1 = shave(uint8(single(im_cr_gnd) * 255), [up_scale, up_scale]);
            ycbcr_gnd = cat(3,(im_y_gnd1),(im_cb_gnd1),(im_cr_gnd1));
            im_gnd1 = ycbcr2rgb(ycbcr_gnd);
        else
            im_h1 = im_y_h1;
            im_gnd1 = im_y_gnd1;
        end
        im_h_set{iii} = im_h1;
        im_gnd_set{iii} = im_gnd1;

        %% save images
        imwrite(im_h1,fullfile(folderResultCur,[imageName,'_x',num2str(up_scale),'.png']));

        %% compute PSNR and SSIM and IFC
        bic(1) = compute_psnr(im_y_gnd1,im_y_b1);
        drrn(1) = compute_psnr(im_y_gnd1,im_y_h1);
        bic(2) = ssim_index(im_y_gnd1,im_y_b1);
        drrn(2) = ssim_index(im_y_gnd1,im_y_h1);
        bic(3) = ifcvec(double(im_y_gnd1),double(im_y_b1));
        drrn(3) = ifcvec(double(im_y_gnd1),double(im_y_h1));

        bicubic_set = [bicubic_set; bic];
        drrn_set = [drrn_set; drrn];
    end
    mean_bicubic = [mean_bicubic; [mean(bicubic_set(:,1)) mean(bicubic_set(:,2)) mean(bicubic_set(:,3))]];
    mean_drrn = [mean_drrn; [mean(drrn_set(:,1)) mean(drrn_set(:,2)) mean(drrn_set(:,3))]];
end

%%% save PSNR and SSIM metrics
PSNR_set = drrn_set(:,1);
SSIM_set = drrn_set(:,2);
IFC_set = drrn_set(:,3);
save(fullfile(folderResultCur,['PSNR_',setTestCur,'_x',num2str(up_scale),'.mat']),['PSNR_set'])
save(fullfile(folderResultCur,['SSIM_',setTestCur,'_x',num2str(up_scale),'.mat']),['SSIM_set'])
save(fullfile(folderResultCur,['IFC_',setTestCur,'_x',num2str(up_scale),'.mat']),['IFC_set'])

count = 1;
for kk = startNUM:thresh:endNUM
    disp(['epoch: ' num2str(kk) '---- bic = ' num2str(mean_bicubic(count,:)) '---- DRRN = ' num2str(mean_drrn(count,:))]);
    count = count+1;
end

end

author: tyshiwo


text:
function handles = dtiLoadROI(hObject,handles,fileName);
%
%  handles = dtiLoadROI(hObject,handles,[fileName]);
%
%Author: Wandell, Dougherty
%Purpose:
%

persistent defaultPath;

if(~exist('fileName','var') | isempty(fileName))
    if(isempty(defaultPath))
        fn = handles.defaultPath;
        if(exist(fullfile(fn, 'ROIs'),'dir')) fn = fullfile(fn, 'ROIs', filesep); end
    else
        fn = defaultPath;
    end
    [f, p] = uigetfile({'*.mat';'*.*'}, 'Load ROI...', fn);
    if(isnumeric(f)), disp('Load ROI canceled.'); return; end
    defaultPath = p;
    fileName = fullfile(p,f);
    refreshFig = 1;
else
    refreshFig = 0;
end
if(isfield(handles, 't1NormParams'))
  roi = dtiReadRoi(fileName, handles.t1NormParams);
else
  roi = dtiReadRoi(fileName);
end

handles = dtiAddROI(roi,handles);
if(refreshFig) handles = dtiRefreshFigure(handles, 0); end
if(~isempty(hObject)) guidata(hObject, handles); end
return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [dpixc, dveccr, N] = generatedata(d, sizevec, psf, maxphot, offset, Nt, rs)

nx = sizevec(2)-sizevec(1);
ny = sizevec(4)-sizevec(3);
N = size(d,1);

blinkmat = rand(N, Nt);

for ii=1 : N
    dpix_ind(:,:,ii) = pixelize(d(ii,:), 1, sizevec, nx, ny, [],0);
    dpixc_ind(:,:,ii) = conv2(dpix_ind(:,:,ii),psf,'same'); %convolution of individual points
end

dvecc_ind = squeeze(reshape(dpixc_ind, nx*ny, 1, N)); % each indiv image to vector

dvecc_nonoise = (dvecc_ind*blinkmat); % set of Nt vectors
dpixc_nonoise = imresize(reshape(dvecc_nonoise, nx, ny, Nt), rs);
dpixc_nonoise_dip = dip_image(dpixc_nonoise);

maxdvec = max(dpixc_nonoise(:));
offset_abs = offset/(1-offset)*maxdvec;
dpixc_nonoise_dip = (dpixc_nonoise_dip + offset_abs)/(maxdvec+offset_abs)*maxphot;
dpixc_dip = noise(dpixc_nonoise_dip,'poisson');
dpixc = double(dpixc_dip);
dveccr = double(squeeze(reshape(dpixc, nx*ny*rs^2, 1, Nt))); % vectors of resized images

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function mm_size_print ( input_filename, rep, symm, nrow, ncol, nnz )

%*****************************************************************************80
%
%% MM_SIZE_PRINT prints size information from a Matrix Market file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 April 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character ( len = * ) INPUT_FILENAME, the input file name.
%
%    Input, character ( len = 10 ) REP, the Matrix Market 'representation'
%    indicator.  Possible values include:
%    'coordinate'   (for sparse data)
%    'array'        (for dense data)
%    'elemental'    (to be added)
%
%    Input, character ( len = 19 ) SYMM, the Matrix Market symmetry.
%    Possible values include:
%    'symmetric'
%    'hermitian'
%    'skew-symmetric'
%    'general'
%
%    Input, integer NROW, the number of rows in the matrix.
%
%    Input, integer NCOL, the number of columns in the matrix.
%
%    Input, integer NNZ, the number of nonzero entries required to store
%    the matrix.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MM_SIZE_PRINT:\n' );
  fprintf ( 1, '  Size information from Matrix Market file "%s"\n', ...
    input_filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of rows    NROW = %d\n', nrow );
  fprintf ( 1, '  Number of columns NCOL = %d\n', ncol );

  if ( s_eqi ( rep, 'coordinate' ) )
    fprintf ( 1, '  Declared number of nonzeros NNZ = %d\n', nnz );
  elseif ( s_eqi ( rep, 'array' ) )
    nnz = mm_nnz_set ( rep, symm, nrow, ncol );
    fprintf ( 1, '  Inferred number of nonzeros NNZ = %d\n', nnz );
  end

  return
end


text:
function varargout = AffineTrnasformation(varargin)
% AFFINETRNASFORMATION M-file for AffineTrnasformation.fig
% Affine_Transf computes and applies the geometric affine transformation to a 2-D image.
%
% The program main functions are:
% - Load Image: Load the image to be transformed.
% - Transform Image: Computes the transformation matrix from the transformation parameters specified by the user, then it applies the transformation to the loaded image.
% - Save Image: Save the transformed image.
% - Quit: Exit the program.
% The transformation parameters are specified by the user by editing the
% values of the rotation, scale (x and y), and shear (x and y). The
% translation is not icluded because it can be easily (normalized or) removed from image by
% subtracting the shape cetroid.
%
% The affine matrix and the transformation are computed using Matlab Image Processing Toolbox
% functions: maketform and imtransform, so this program is only demonstrates the computation and the applying of these functions to a 2-D image.
%
% Version 1.0
% Copyright Feb 2003, Ibrahim El Rube
% ielrube@yahoo.com

% Begin initialization code - DO NOT EDIT
global gg gc;
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @AffineTrnasformation_OpeningFcn, ...
                   'gui_OutputFcn',  @AffineTrnasformation_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin & isstr(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before AffineTrnasformation is made visible.
function AffineTrnasformation_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to AffineTrnasformation (see VARARGIN)

% Choose default command line output for AffineTrnasformation
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes AffineTrnasformation wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = AffineTrnasformation_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global gg
% Get default command line output from handles structure
varargout{1} = handles.output;

gg=2;
% --- Executes on button press in Load_image.
function Load_image_Callback(hObject, eventdata, handles)
% hObject    handle to Load_image (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global gg gc;

h=guidata(gcbo);
set(h.message_text,'String','');  % Clear Messages
gg=gg+2;
h.image_filename=1;
h.image_pathname=1;
[filename, pathname] = uigetfile( ...
       {'*.bmp;*.tif;*.jpg;*.pcx;*.png;*.hdf;*.xwd;*.ras;*.pbm;*.pgm;*.ppm;*.pnm', 'All MATLAB SUPPORTED IMAGE Files (*.bmp,*.tif,*.jpg,*.pcx,*.png,*.hdf,*.xwd,*.ras,*.pbm,.pgm,*.ppm,*.pnm)'} ...
        ,'Pick a file');     % Load Image file and path names
    if filename~=0
        gg=1;
        h.gg=1;
        h.image_filename=filename;  % Image file name
        h.image_pathname=pathname;  % Image path name
        set(h.Showaxes,'visible','off');
        axes(h.axes1);
        set(h.figure1_title,'Visible','on');
        image_1=imread([pathname filename]); % Load Image
        imshow(image_1);                     % Show Loaded Image
        h.image_1=image_1;
        set(h.axes2,'Visible','off');
        axes(h.axes2);
        cla;
        set(h.figure2_title,'Visible','off');
        gc=0;
        h.axis12=0;
        guidata(gcbo,h);
    end
% --- Executes during object creation, after setting all properties.
function xscale_CreateFcn(hObject, eventdata, handles)
% hObject    handle to xscale (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

function xscale_Callback(hObject, eventdata, handles)
% hObject    handle to xscale (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of xscale as text
%        str2double(get(hObject,'String')) returns contents of xscale as a
%        double

% --- Executes during object creation, after setting all properties.
function rot_CreateFcn(hObject, eventdata, handles)
% hObject    handle to rot (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

function rot_Callback(hObject, eventdata, handles)
% hObject    handle to rot (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of rot as text
%        str2double(get(hObject,'String')) returns contents of rot as a double

% --- Executes during object creation, after setting all properties.
function xshear_CreateFcn(hObject, eventdata, handles)
% hObject    handle to xshear (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

function xshear_Callback(hObject, eventdata, handles)
% hObject    handle to xshear (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of xshear as text
%        str2double(get(hObject,'String')) returns contents of xshear as a double

% --- Executes during object creation, after setting all properties.
function xtrans_CreateFcn(hObject, eventdata, handles)
% hObject    handle to xtrans (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

function xtrans_Callback(hObject, eventdata, handles)
% hObject    handle to xtrans (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of xtrans as text
%        str2double(get(hObject,'String')) returns contents of xtrans as a double

% --- Executes on button press in deg.
function deg_Callback(hObject, eventdata, handles)
% hObject    handle to deg (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of deg
h=guidata(gcbo);
%set(h.rad,'Value',0);
ch=get(hObject,'Value');
switch ch
case 0
   set(h.rad,'Value',1)
case 1
   set(h.rad,'Value',0)
end
guidata(gcbo,h);
% --- Executes on button press in rad.
function rad_Callback(hObject, eventdata, handles)
% hObject    handle to rad (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of rad
h=guidata(gcbo);
%set(h.rad,'Value',0);
ch=get(hObject,'Value');
switch ch
case 0
  set(h.deg,'Value',1)
case 1
  set(h.deg,'Value',0)
end
guidata(gcbo,h);
% --- Executes on button press in bshear.
function bshear_Callback(hObject, eventdata, handles)
% hObject    handle to bshear (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of bshear

% --- Executes on button press in checkbox2.
function checkbox2_Callback(hObject, eventdata, handles)
% hObject    handle to checkbox2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of checkbox2

% --- Executes during object creation, after setting all properties.
function ytrans_CreateFcn(hObject, eventdata, handles)
% hObject    handle to ytrans (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

function ytrans_Callback(hObject, eventdata, handles)
% hObject    handle to ytrans (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of ytrans as text
%        str2double(get(hObject,'String')) returns contents of ytrans as a double

% --- Executes on button press in btrans.
function btrans_Callback(hObject, eventdata, handles)
% hObject    handle to btrans (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of btrans

% --- Executes during object creation, after setting all properties.
function yscale_CreateFcn(hObject, eventdata, handles)
% hObject    handle to yscale (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

function yscale_Callback(hObject, eventdata, handles)
% hObject    handle to yscale (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of yscale as text
%        str2double(get(hObject,'String')) returns contents of yscale as a double

% --- Executes on button press in bscale.
function bscale_Callback(hObject, eventdata, handles)
% hObject    handle to bscale (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of bscale

% --- Executes on button press in Transform.
function Transform_Callback(hObject, eventdata, handles)
% hObject    handle to Transform (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global gg gc;
gc=0;
h=guidata(gcbo);
if  gg~=2
    h=guidata(gcbo);
    image_1=h.image_1;  % Get Image 1
    hb(1)=num_check(h.rot); % Rotation
    hb(2)=num_check(h.xscale); % Scale in x-direction
    hb(3)=num_check(h.yscale); % Scale in y-direction
    hb(4)=num_check(h.xshear); % Shear
    hb(5)=num_check(h.yshear);
    hc=sum(hb);
    switch hc
        case 0
            a=str2num(get(h.xscale,'String'));
            b=str2num(get(h.xshear,'String'));
            c=str2num(get(h.yshear,'String'));
            d=str2num(get(h.yscale,'String'));
            th=str2num(get(h.rot,'String'));  %2*pi/10;
            ch=get(h.deg,'Value');
            switch ch
                case 0
                    th=th;
                otherwise
                    th=th*pi/180;;
            end
            k=1;
            TT=k*[cos(th) -sin(th); sin(th) cos(th)]*[a b; c d];
            TTT=[ 1 0 0; 0 1 0; 0 0 1];
            TTT(1:2,1:2)=TT;
            set(h.message_text,'String','Transformaing .... Please Wait');
            tform=maketform('affine',TTT);
            [image_2] = imtransform(image_1,tform);
            h.image_2=image_2;
            axes(h.axes2);
            set(h.figure2_title,'Visible','on');
            imshow(image_2,[]);
            gc=1;
            axes(h.axes2);
            g2=h.axis12;
            switch g2
                case 0
                    axes(h.axes1);
                    axis on;
                    axes(h.axes2);
                    axis on;
                otherwise
                    axes(h.axes1);
                    axis off;
                    axes(h.axes2);
                    axis off;
            end
            set(h.message_text,'String','Done');
            set(h.Showaxes,'Visible','on');
            guidata(gcbo,h);
        otherwise
            set(h.message_text,'String','Invalid value: enter a real number only');
    end
end

% --- Executes on button press in Save_image.
function Save_image_Callback(hObject, eventdata, handles)
% hObject    handle to Save_image (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global gg gc;
h=guidata(gcbo);
set(h.message_text,'String','');
if   gc~=0
    dot='.';
    pathname=h.image_pathname;
    [filename, pathname] = uiputfile('*.*', 'Save Transformed Image');
    if filename~=0
        try_again=find(filename==dot);
        if ~isempty(try_again)
            image_ex=filename(try_again+1:end);
            valid_ex=['BMP'; 'TIF'; 'JPG'; 'PCX'; 'PNG'; 'HDF'; 'XWD'; 'RAS'; 'PBM'; 'PGM'; 'PPM'; 'PNM';];
            p=['\<' image_ex '\>'];
            g=regexpi(valid_ex,p);  % Match Exact extension, ignoring case.
            b=cell2mat(g);           % Convert cell to matrix
            if b~=0 & length(image_ex)==3
                  CMap=colormap(h.axes1);
                  h.savedimage_filename=filename;
                  h.savedimage_pathname=pathname;
                  imwrite(h.image_2,CMap,[pathname filename]);
                  savedto=['Image Saved to:    ' pathname filename];
                  set(h.message_text,'String',savedto);
             else
                  set(h.message_text,'String','Invalid file extension: Image NOT saved, use one of these extensions, BMP TIF JPG PCX PNG HDF XWD RAS PBM PGM PPM PNM');
            end
        else
            set(h.message_text,'String','Empty file extension: Image NOT saved, use one of these extensions, BMP TIF JPG PCX PNG HDF XWD RAS PBM PGM PPM PNM');
        end

    end
end
% --- Executes on button press in pushbutton4.
function pushbutton4_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton4 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
close;

% --- Executes during object creation, after setting all properties.
function yshear_CreateFcn(hObject, eventdata, handles)
% hObject    handle to yshear (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

function yshear_Callback(hObject, eventdata, handles)
% hObject    handle to yshear (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of yshear as text
%        str2double(get(hObject,'String')) returns contents of yshear as a double

function hb=num_check(ha);
% This function checks the validation of the entered transformation
% parameters
% ha is the handle to be checked
% hb is a flag, 1 for correct range, 0 for incorrect range.
h=guidata(gcbo);
ba=get(ha,'String');
bb=str2num(ba);
hb=0;
if isempty(bb) | ba=='i' | ba=='j'
    axes(h.axes2);
    cla;
    set(h.message_text,'String','Invalid value: enter a real number only');
    switch ha
        case h.xscale
            set(ha,'String','1');
        case h.yscale
            set(ha,'String','1');
        otherwise
            set(ha,'String','0');
    end
    hb=1;
    %break
end
%hb=ha;
%guidata(gcbo,h);

% --- Executes on button press in Showaxes.
function Showaxes_Callback(hObject, eventdata, handles)
% hObject    handle to Showaxes (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
h=guidata(gcbo);
sh=get(h.Showaxes,'String');
switch sh
    case 'Show Axes'
        set(h.Showaxes,'String','Hide Axes');
        h.axis12=0;
    otherwise
        set(h.Showaxes,'String','Show Axes');
        h.axis12=1;
end
g1=get(h.axes1,'visible');
axes(h.axes1);
switch g1
    case 'off'
        axis on;
        set(h.axes1,'visible','on');
    otherwise
        axis off;
        set(h.axes1,'visible','off');
end
g2=get(h.axes2,'visible');
axes(h.axes2);
switch g2
    case 'off'
        axis on;
        set(h.axes2,'visible','on');
    otherwise
        axis off;
        set(h.axes2,'visible','off');
end;
guidata(gcbo,h);


text:
function outputImage = shuffle(this, iDimArray, shuffledSamplingPointsIndexArray)
% Shuffles sampling points of specified dimensions, and dimInfo-description
% (?)
%
%   Y = MrDataNd()
%   Y.shuffle(inputs)
%
% This is a method of class MrDataNd.
%
% IN
%
% OUT
%
% EXAMPLE
%   shuffle({'y','z'}, {[nSamples:-1:1], [1:2:nSamples, 2:2:nSamples]}),
%
%   See also MrDataNd

% Author:   Saskia Bollmann & Lars Kasper
% Created:  2019-10-10
% Copyright (C) 2019 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public License (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

outputImage = this.copyobj;

nDims = numel(iDimArray);
for iDim = 1:nDims
    idxDim = this.dimInfo.get_dim_index(iDimArray{iDim});
    dataSelectorString = repmat({':'}, 1, this.dimInfo.nDims);
    dataSelectorStringNew = dataSelectorString;
    dataSelectorStringNew{idxDim} = shuffledSamplingPointsIndexArray{iDim};
    outputImage.data(dataSelectorString{:}) = ...
        outputImage.data(dataSelectorStringNew{:});
    % TODO: do we adapt the dimInfo as well?
end

author: translationalneuromodeling

repo: tapas


text:
function fem_basis_test06 ( )

%*****************************************************************************80
%
%% FEM_BASIS_TEST06 repeats TEST03 using FEM_BASIS_MD.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 October 2010
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FEM_BASIS_TEST06\n' );
  fprintf ( 1, '  FEM_BASIS_MD evaluates an arbitrary\n' );
  fprintf ( 1, '  basis function over an M-dimensional simplex.\n' );

  m = 3;

  i1 = [ 1, 0, 2, 1 ]';
  d = sum ( i1 );
  x1(1:m) = i1(1:m) / d;
  fprintf ( 1, '\n' );
  fprintf ( 1, '   I   J   K   L        X           Y           Z      L(I,J,K,L)(X,Y,Z)\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  %2d  %2d  %2d  %2d  %10.4f  %10.4f  %10.4f  %14.6g\n', ...
    i1(1:m+1), x1(1:m), 1.0 );
  fprintf ( 1, '\n' );
  for p3 = 0 : d
    i2(3) = p3;
    for p2 = 0 : d - p3
      i2(2) = p2;
      for p1 = 0 : d - p2 - p3
        i2(1) = p1;
        i2(m+1) = d - sum ( i2(1:m) );
        x2(1:m) = i2(1:m) / d;
        l = fem_basis_md ( m, i1, x2 );
        fprintf ( 1, '  %2d  %2d  %2d  %2d  %10.4f  %10.4f  %10.4f  %14.6g\n', ...
          i2(1:m+1), x2(1:m), l );
      end
    end
  end

  return
end


text:
function Y = AtA_Toeplitz(X,Lambda);
% AtA_Toeplitz: Gram Operator of SeparateAngle
%  Usage:
%     Y = AtA(X,Lambda)
%  Inputs:
%    X      n*n matrix (x,y)
%    Lambda Fourier multipliers
%  Outputs:
%    Y      n*n matrix (x,y)
%  Description:
%    Performs  A'A where A is the angular windowing using the
%    Toeplitz structure of A'A. AtA_Toeplitz uses FFTs.
%
% By Emmanuel Candes, 2003-2004

  n  = size(X,1); n2 = n/2;

  X1 = zeros(n); X2 = zeros(n);

  [ix,w] = DetailMeyerWindow([n2/4 n2/2],3);
  lx = reverse(-ix);

  F = ifft_mid0(X).*sqrt(n); % Take FT

  for r = 1:(n2/2),
    k = lx(r); t = n2 + k + 1;
    X1(:,t) = GUSFT_Toeplitz(F(:,t),Lambda(:,r,1));

    rsym = n2/2+1-r; t = n2 - k + 1;
    X1(:,t) = GUSFT_Toeplitz(F(:,t),Lambda(:,rsym,2));
  end

  X1 = fft_mid0(X1)/sqrt(n); % Invert FT

  F = ifft_mid0(X.').*sqrt(n); % Take FT

  for r = 1:(n2/2),
    k = lx(r); t = n2 + k + 1;
    X2(:,t) = GUSFT_Toeplitz(F(:,t),Lambda(:,r,1));

    rsym = n2/2+1-r;  t = n2 - k + 1;
    X2(:,t) = GUSFT_Toeplitz(F(:,t),Lambda(:,rsym,2));
  end

  X2 = fft_mid0(X2)/sqrt(n); % Invert FT

  Y = X1 + X2.';

author: thomaskuestner


text:
function [dLdp,g,x] = spm_mci_adjoint (Pr,M,U,Y)
% Gradient of log joint from adjoint method
% FORMAT [dLdp,g,x] = spm_mci_adjoint (Pr,M,U,Y)
%
% Pr        Parameters (vectorised and in M.V subspace)
% M         Model structure
% U         Inputs  [Nin x N]
% Y         Data
%
% dLdp      Gradient    [Np x 1]
% g         Outputs     [N x Nout]
% x         States      [N x Nstates]
%
% If M.adjlike=1 this function returns gradient of log likelihood
%
% This function uses integrators from MATLAB's ODE Suite
%
% B. Sengupta, K. Friston and W. Penny (2014) Efficient Gradient
% Computation for Dynamical Models. Neuroimage,98, 521-527.
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Will Penny and Biswa Sengupta
% $Id: spm_mci_adjoint.m 6697 2016-01-27 14:57:28Z spm $

try, adjlike=M.adjlike; catch, adjlike=0; end

% Parameters in original space
P = M.V*Pr+M.vpE;

if isempty(U)
    U=zeros(1,M.N);
end

[g,x] = spm_mci_fwd(P,M,U);

% Gradient at computed times
% When computing output sensitivities, assume dydx=L, ie not a
% function of x. Generalise later
[tmp,L]=feval(M.g,M.x0,U(:,1),P,M);
e=Y-g;
djdx=e*M.iCe*L;

% integrate adjoint equation
lambda = spm_mci_adjoint_int(U,P,M,x,djdx);

% If observation function becomes dependent on parameters
% djdp will need updating
djdp = zeros(1,M.Np);

dLdp=zeros(1,M.Np);
for n=1:M.N,
    % Evaluate parameter Jacobian
    if isfield(M,'dfdp')
        Fp = feval(M.dfdp,x(n,:)',U(:,n),P,M);
    else
        Fp = spm_diff(M.f,x(n,:)',U(:,n),P,M,3);
    end
    dLdp=dLdp+djdp-lambda(n,:)*Fp;
end

if ~adjlike
    dlogpriordp = spm_mci_gprior_deriv (Pr,M);
    dLdp=dLdp+dlogpriordp;
end

author: spm


text:
function [fMshift, fProbability, fBic, mProblikelihood, fAICc] = calc_loglikelihood_dM(mCat1, mCat2)
% function [fMshift, fProbability, fBic, mProblikelihood] = calc_loglikelihood_dM(mCat1, mCat2);
% ----------------------------------------------------------------------------------------------
% Calculate log-likelihood estimation of magnitude shift dM between to periods
%
% Incoming variable
% mCat1 : EQ catalog period 1 (Catalog to be modified)
% mCat2 : EQ catalog period 2 (Observed catalog)
%
% Outgoing variable
% fProbability : log-likelihood probabilty
% fMshift      : Magnitude shift with the lowest  max. lieklihood score
% fBic         : Bayesian Information Criterion value
% mProblikelihood : Solution matrix shift and likelihood score
%
% Author: J. Woessner, woessner@seismo.ifg,.ethz.ch
% last update: 28.10.02

% Initialize
mProblikelihood = [];
vfProbability = [];
vMshift = [];

% Determine exact time period
fPeriod1 = max(mCat1(:,3)) - min(mCat1(:,3));
fPeriod2 = max(mCat2(:,3)) - min(mCat2(:,3));

mCat1Mod = mCat1;

for fMshift = -0.5:0.1:0.5
    % Apply shift
    mCat1Mod(:,6) = mCat1Mod(:,6)+fMshift;
    % Initialize values
    fMinMag = min([min(mCat1Mod(:,6)) min(mCat2(:,6))]);
    fMaxMag = max([max(mCat1Mod(:,6)) max(mCat2(:,6))]);

    [vPredFMD,vBin1] = hist(mCat1Mod(:,6),0:0.1:fMaxMag);
    [vObsFMD,vBin2] = hist(mCat2(:,6),0:0.1:fMaxMag);
    % Time normalization and round due to Poisson distribution calculation in calc_log10poisspdf
    vPredFMD = ceil(vPredFMD./fPeriod1);
    vObsFMD = ceil(vObsFMD./fPeriod2);
    % Calculate the likelihoods for both models
    vProb_ = calc_log10poisspdf(vObsFMD',vPredFMD');
    % Sum the probabilities
    fProbability = (-1) * sum(vProb_);
    vfProbability = [vfProbability; fProbability];
    vMshift = [vMshift; fMshift];
    mCat1Mod = mCat1;
end

%%% Find the minimum loglikelihodd score: if the minimum score is obtained several times, calculate MEAN
%%% of the magnitude shift
vdMloglikeli = [vfProbability vMshift];
vSel = (vdMloglikeli == min(vdMloglikeli(:,1)));
vdMloglikeli = vdMloglikeli(vSel,:);
if length(vdMloglikeli(:,1)) > 1
    fProbability = min(vdMloglikeli(:,1));
    fMshift = mean(vdMloglikeli(:,2));
else
    fProbability = vdMloglikeli(:,1);
    fMshift = vdMloglikeli(:,2);
end
% Solution matrix
mProblikelihood = [vMshift vfProbability];

nDegFree = 1; % Magnitude shift is the degree of freedom
n_samples = length(mCat1(:,6))+length(mCat2(:,6));
%% Bayesian Information Criterion (BIC)
fBic = 2*fProbability + 2*log(n_samples)*nDegFree;
%% Corrected Akaike Information Criterion (AICc)
fAICc = -2*(-fProbability)+2*nDegFree+2*nDegFree*(nDegFree+1)/(n_samples-nDegFree-1);

author: CelsoReyes


text:
%% This is for testing the Homogeneous Transformation functions in the robotics Toolbox

function tests = PluckerTest
  tests = functiontests(localfunctions);

  clc
end

function teardownOnce(tc)
    close all
end

% Primitives
function constructor1_test(tc)

    % construct from 6-vector
    L = Plucker([1 2 3 4 5 6]);
    tc.verifyTrue( isa(L, 'Plucker') );
    tc.verifyEqual(L.v, [1 2 3]','absTol',1e-10);
    tc.verifyEqual(L.w, [4 5 6]','absTol',1e-10);

    % construct from object
    L2 = Plucker(L)
    tc.verifyTrue( isa(L, 'Plucker') );
    tc.verifyEqual(L2.v, [1 2 3]','absTol',1e-10);
    tc.verifyEqual(L2.w, [4 5 6]','absTol',1e-10);

    % construct from point and direction
    L = Plucker.pointdir([1 2 3], [4 5 6]);
    tc.verifyTrue(L.contains([1 2 3]'));
    tc.verifyEqual(L.uw, unit([4 5 6])', 'absTol', 1e-10);
end

function double_test(tc)
    % verify double
    L = Plucker([1 2 3 4 5 6]);
    tc.verifyEqual(double(L), [1 2 3 4 5 6]', 'absTol', 1e-10);
end

function constructor2_test(tc)
    % 2 point constructor
    P = [2 3 7]; Q = [2 1 0];
    L = Plucker(P, Q);
    tc.verifyEqual(L.w', P-Q,'absTol',1e-10);
    tc.verifyEqual(L.v', cross(P-Q,Q),'absTol',1e-10);

    % test all possible input shapes
    L2 = Plucker(P', Q);
    tc.verifyEqual(double(L2), double(L), 'absTol',1e-10);
    L2 = Plucker(P, Q');
    tc.verifyEqual(double(L2), double(L), 'absTol',1e-10);
    L2 = Plucker(P', Q');
    tc.verifyEqual(double(L2), double(L), 'absTol',1e-10);
    L2 = Plucker(P, Q);
    tc.verifyEqual(double(L2), double(L), 'absTol',1e-10);

    % planes constructor
    P = [10 11 12]'; w = [1 2 3];
    L = Plucker.pointdir(P, w);
    tc.verifyEqual(double(L), [cross(w,P) w]','absTol',1e-10); %FAIL
    L2 = Plucker.pointdir(P', w);
    tc.verifyEqual(double(L2), double(L), 'absTol',1e-10);
    L2 = Plucker.pointdir(P, w');
    tc.verifyEqual(double(L2), double(L), 'absTol',1e-10);
    L2 = Plucker.pointdir(P', w');
    tc.verifyEqual(double(L2), double(L), 'absTol',1e-10);
end

function pp_test(tc)
    % validate pp and ppd
    L = Plucker([-1 1 2], [1 1 2]);
    tc.verifyEqual(L.pp, [0 1 2]', 'absTol',1e-10);
    tc.verifyEqual(L.ppd, sqrt(5), 'absTol',1e-10);

    % validate pp
    tc.verifyTrue( L.contains(L.pp) );
end

function contains_test(tc)
    P = [2 3 7]; Q = [2 1 0];
    L = Plucker(P, Q);

    % validate contains
    tc.verifyTrue( L.contains([2 3 7]') );
    tc.verifyTrue( L.contains([2 1 0]') );
    tc.verifyFalse( L.contains([2 1 4]') );
end

function closest_test(tc)
    P = [2 3 7]; Q = [2 1 0];
    L = Plucker(P, Q);

    [p,d] = L.closest(P);
    tc.verifyEqual(p, P', 'absTol',1e-10);
    tc.verifyEqual(d, 0, 'absTol',1e-10);

     % validate closest with given points and origin
    [p,d] = L.closest(Q);
    tc.verifyEqual(p, Q', 'absTol',1e-10);
    tc.verifyEqual(d, 0, 'absTol',1e-10);

    L = Plucker([-1 1 2], [1 1 2]);
    [p,d] = L.closest([0 1 2]);
    tc.verifyEqual(p, [0 1 2]', 'absTol',1e-10);
    tc.verifyEqual(d, 0, 'absTol',1e-10);

    [p,d] = L.closest([5 1 2]);
    tc.verifyEqual(p, [5 1 2]', 'absTol',1e-10);
    tc.verifyEqual(d, 0, 'absTol',1e-10);

    [p,d] = L.closest([0 0 0]);
    tc.verifyEqual(p, L.pp, 'absTol',1e-10);
    tc.verifyEqual(d, L.ppd, 'absTol',1e-10);

    [p,d] = L.closest([5 1 0]);
    tc.verifyEqual(p, [5 1 2]', 'absTol',1e-10);
    tc.verifyEqual(d, 2, 'absTol',1e-10);
end

function plot_test(tc)

    P = [2 3 7]; Q = [2 1 0];
    L = Plucker(P, Q);

    axis(10*[-1 1 -1 1 -1 1]);
    plot(L, 'r', 'LineWidth', 2);
end

function eq_test(tc)
    w = [1 2 3];
    P = [-2 4 3];

    L1 = Plucker(P, P+w);
    L2 = Plucker(P+2*w, P+5*w);
    L3 = Plucker(P+[1 0 0], P+w);

    tc.verifyTrue(L1 == L2);
    tc.verifyFalse(L1 == L3);

    tc.verifyFalse(L1 ~= L2);
    tc.verifyTrue(L1 ~= L3);
end

function skew_test(tc)

    P = [2 3 7]; Q = [2 1 0];
    L = Plucker(P, Q);

    m = L.skew;

    tc.verifyEqual(size(m), [4 4]);
    tc.verifyEqual(m+m', zeros(4,4), 'abstol', 1e-10);
end

function mtimes_test(tc)
    P = [1 2 0]; Q = [1 2 10];  % vertical line through (1,2)
    L = Plucker(P, Q);

    % check pre/post multiply by matrix
    M = rand(4,10);

    a = L * M;
    tc.verifyEqual(a, L.skew*M);

    M = rand(10,4);
    a = M * L;
    tc.verifyEqual(a, M*L.skew);

    % check transformation by SE3

    L2 = SE3 * L;
    tc.verifyEqual(L.double, L2.double);

    L2 = SE3(2, 3, 1) * L; % shift line in the xy directions
    pxy = L2.intersect_plane([0 0 1 0]);
    tc.verifyEqual(pxy, [1+2 2+3 0]');

end

function parallel_test(tc)

    L1 = Plucker.pointdir([4 5 6], [1 2 3]);
    L2 = Plucker.pointdir([5 5 6], [1 2 3]);
    L3 = Plucker.pointdir([4 5 6], [3 2 1]);

    % L1 || L2 but doesnt intersect
    % L1 intersects L3

    tc.verifyTrue( isparallel(L1, L1) );
    tc.verifyTrue(L1 | L1);

    tc.verifyTrue( isparallel(L1, L2) );
    tc.verifyTrue(L1 | L2);
    tc.verifyTrue( isparallel(L2, L1) );
    tc.verifyTrue(L2 | L1);
    tc.verifyFalse( isparallel(L1, L3) );
    tc.verifyFalse(L1 | L3);
end

function intersect_test(tc)

    L1 = Plucker.pointdir([4 5 6], [1 2 3]);
    L2 = Plucker.pointdir([5 5 6], [1 2 3]);
    L3 = Plucker.pointdir( [4 5 6], [0 0 1]);
    L4 = Plucker.pointdir([5 5 6], [1 0 0]);

    % L1 || L2 but doesnt intersect
    % L3 intersects L4
    tc.verifyFalse( L1^L2 );

    tc.verifyTrue( L3^L4 );

end

function commonperp_test(tc)
    L1 = Plucker.pointdir([4 5 6], [0 0 1]);
    L2 = Plucker.pointdir([6 5 6], [0 1 0]);

    tc.verifyFalse( L1|L2 );
    tc.verifyFalse( L1^L2 );

    tc.verifyEqual( distance(L1, L2), 2);

    L = commonperp(L1, L2);  % common perp intersects both lines

    tc.verifyTrue( L^L1 );
    tc.verifyTrue( L^L2 );
end

function line_test(tc)

    % mindist
    % intersect
    % char
    % intersect_volume
    % mindist
    % mtimes
    % or
    % side
end

function point_test(tc)
    P = [2 3 7]; Q = [2 1 0];
    L = Plucker(P, Q);

    tc.verifyTrue( L.contains(L.point(0)) );
    tc.verifyTrue( L.contains(L.point(1)) );
    tc.verifyTrue( L.contains(L.point(-1)) );
end

function char_test(tc)
    P = [2 3 7]; Q = [2 1 0];
    L = Plucker(P, Q);

    s = char(L);
    tc.verifyClass(s, 'char');
    tc.verifyEqual(size(s,1), 1);

    s = char([L L]);
    tc.verifyClass(s, 'char');
    tc.verifyEqual(size(s,1), 2);
    s = char([L L]');
    tc.verifyClass(s, 'char');
    tc.verifyEqual(size(s,1), 2);
end

function plane_test(tc)

    xyplane = [0 0 1 0];
    xzplane = [0 1 0 0];
    L = Plucker.planes(xyplane, xzplane); % x axis
    tc.verifyEqual(double(L), [0 0 0 -1 0 0]','absTol',1e-10);

    L = Plucker([-1 2 3], [1 2 3]);  %line at y=2,z=3
    x6 = [1 0 0 -6]; % x = 6

    % plane_intersect
    [p,lambda] = L.intersect_plane(x6)
    tc.verifyEqual(p, [6 2 3]','absTol',1e-10);
    tc.verifyEqual(L.point(lambda), [6 2 3]','absTol',1e-10);

    x6s.n = [1 0 0];
    x6s.p = [6 0 0];
    [p,lambda] = L.intersect_plane(x6s);
    tc.verifyEqual(p, [6 2 3]','absTol',1e-10);
    tc.verifyEqual(L.point(lambda), [6 2 3]','absTol',1e-10);
end

function methods_test(tc)
    % intersection
    px = Plucker([0 0 0], [1 0 0]);  % x-axis
    py = Plucker([0 0 0], [0 1 0]);  % y-axis
    px1 = Plucker([0 1 0], [1 1 0]); % offset x-axis

    verifyEqual(tc, px.ppd, 0);
    verifyEqual(tc, px1.ppd, 1);
    verifyEqual(tc, px1.pp, [0 1 0]');

    px.intersects(px)
    px.intersects(py)
    px.intersects(px1)
end

% function intersect_test(tc)
%     px = Plucker([0 0 0], [1 0 0]);  % x-axis
%     py = Plucker([0 0 0], [0 1 0]);  % y-axis
%
%     plane.d = [1 0 0]; plane.p = 2; % plane x=2
%
%     px.intersect_plane(plane)
%     py.intersect_plane(plane)
% end

author: petercorke

repo: spatialmath-matlab


text:
function y = perms(x)
% Linearizes non-vector input, so that it gives the same result as in Matlab; except for order,
% which is different even for vector input
if ~isvector(x)
    x = reshape(x, 1, []);
end
y = builtin('perms', x);
end

author: lmendo

repo: MATL

save_path: github-repos/MATLAB/lmendo-MATL


text:
function [pE, pC, x] = tapas_rdcm_spm_dcm_fmri_priors(A, B, C, D)
% Returns the priors for a two-state DCM for fMRI.
% FORMAT:[pE,pC,x] = spm_dcm_fmri_priors(A,B,C,D,options)
%
%   options.two_state:  (0 or 1) one or two states per region
%   options.endogenous: (0 or 1) exogenous or endogenous fluctuations
%
% INPUT:
%    A,B,C,D - constraints on connections (1 - present, 0 - absent)
%
% OUTPUT:
%    pE     - prior expectations (connections and hemodynamic)
%    pC     - prior covariances  (connections and hemodynamic)
%    x      - prior (initial) states
%__________________________________________________________________________
%
% References for state equations:
% 1. Marreiros AC, Kiebel SJ, Friston KJ. Dynamic causal modelling for
%    fMRI: a two-state model.
%    Neuroimage. 2008 Jan 1;39(1):269-78.
%
% 2. Stephan KE, Kasper L, Harrison LM, Daunizeau J, den Ouden HE,
%    Breakspear M, Friston KJ. Nonlinear dynamic causal models for fMRI.
%    Neuroimage 42:649-662, 2008.
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_dcm_fmri_priors.m 4146 2010-12-23 21:01:39Z karl $

% number of regions
%--------------------------------------------------------------------------
n = length(A);

% check options and D (for nonlinear coupling)
%--------------------------------------------------------------------------
try, options.two_state;  catch, options.two_state  = 0; end
try, options.endogenous; catch, options.endogenous = 0; end
try, D;                  catch, D = zeros(n,n,0);       end

% prior (initial) states and shrinkage priors on A for endogenous DCMs
%--------------------------------------------------------------------------
if options.two_state,  x = sparse(n,6); else, x = sparse(n,5); end
if options.endogenous, a = 128;         else, a = 8;           end

% connectivity priors
%==========================================================================
if options.two_state

    % enforce optimisation of intrinsic (I to E) connections
    %----------------------------------------------------------------------
    A     = (A + eye(n,n)) > 0;

    % prior expectations and variances
    %----------------------------------------------------------------------
    pE.A  =  A*32 - 32;
    pE.B  =  B*0;
    pE.C  =  C*0;
    pE.D  =  D*0;

    % prior covariances
    %----------------------------------------------------------------------
    pC.A  =  A/4;
    pC.B  =  B/4;
    pC.C  =  C*4;
    pC.D  =  D/4;

else

    % enforce self-inhibition
    %----------------------------------------------------------------------
    A     =  A > 0;
    A     =  A - diag(diag(A));

    % prior expectations
    %----------------------------------------------------------------------
    pE.A  =  A/(64*n) - eye(n,n)/2;
    pE.B  =  B*0;
    pE.C  =  C*0;
    pE.D  =  D*0;

    % prior covariances
    %----------------------------------------------------------------------
    pC.A  =  A*a/n + eye(n,n)/(8*n);
    pC.B  =  B;
    pC.C  =  C;
    pC.D  =  D;

end

% and add hemodynamic priors
%==========================================================================
pE.transit = sparse(n,1);
pE.decay   = sparse(n,1);
pE.epsilon = sparse(1,1);

pC.transit = sparse(n,1) + exp(-6);
pC.decay   = sparse(n,1) + exp(-6);
pC.epsilon = sparse(1,1) + exp(-6);

return


text:
function features = feature_extraction( x_flt, CSP, oldSMC, verbose )
% x_flt: cell(2, numsamples) : 2-class
%           x_flt{1, i} = [nChannels, nSamples, nTrials]

if verbose == 1
%     fprintf( '\tFeature Extraction...\n\t\t' );
end

[nclass, nparticles] = size( x_flt );
features = cell( nclass, nparticles );

for i=1:oldSMC.numBands
    if verbose == 1
        if mod(i, 5) == 0
%             fprintf( '%d', i );
        else
%             fprintf( '.' );
        end
        if mod(i, 100) == 0
%             fprintf( '\n' );
        end
    end

    for k=1:nclass
        [nc, ns, nt] = size( x_flt{k, i} );
        temp = reshape( x_flt{k, i}, [nc, ns*nt] );
        temp = CSP{i}.W' * temp;
        temp = reshape( temp, [size(temp, 1), ns, nt] );
        temp = permute( temp, [2 1 3] );
        features{k, i} = squeeze( log(var(temp, 0, 1)) );
    end
end
% fprintf( '\n' );

author: PatternRecognition

repo: OpenBMI


text:
function test_bug1770

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY read_neuralynx_dma read_neuralynx_ncs

filenameA = dccnpath('/home/common/matlab/fieldtrip/data/test/bug1770/2012-07-14_15-33-09');
filenameB = dccnpath('/home/common/matlab/fieldtrip/data/test/bug1770/DigitalLynxRawDataFile.nrd');

hdrA = ft_read_header(filenameA);
hdrB = ft_read_header(filenameB);

begsample = 1;
endsample = 100000;

datA = ft_read_data(filenameA, 'begsample', begsample, 'endsample', endsample);
datB = ft_read_data(filenameB, 'begsample', begsample, 'endsample', endsample);

ch1A = double(datA( 1,:));
ch1B = double(datB(18,:));

figure
subplot(2,1,1);
plot(ft_preproc_standardize(ch1A))
subplot(2,1,2);
plot(-ft_preproc_standardize(ch1B))

ch1As =  ft_preproc_standardize(ch1A);
ch1Bs = -ft_preproc_standardize(ch1B);

figure
hold on
plot(ch1As, 'b.');
plot(ch1Bs, 'r.');

% sofar there seems no reason for concern w.r.t. the reading function
% let us look at the filter kernel and the delay

hA = fft(ch1As);
hB = fft(ch1Bs);
H = hA./hB;
figure
plot(ifft(H));

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function h=m_gshhs(resolution,varargin)
% M_GSHHS Add a coastline to a given map using
%           the Global Self-consistant Hierarchical High-resolution
%           Shorelines, Rivers, and Borders
%
%         M_GSHHS(RES, (standard line option,...,...) ) draws the coastline
%         river network, or borders as  simple lines.
%
%         M_GSHHS(RES,'patch' ( ,standard patch options,...,...) ) draws the
%         coastline as a number of patches (rivers and borders are not
%         arranged so patches can be drawn).
%
%         M_GSHHS(RES,'save',FILENAME) saves the extracted coastline data
%         for the current projection in a file FILENAME. This allows
%         speedier replotting using M_USERCOAST(FILENAME).
%
%         RES: A one-char string (optionally 2 or 3)
%
%         First char: resolution - one of
%                      'c'  crude
%                      'l'  low
%                      'i'  intermediate
%                      'h'  high
%                      'f'  full
%
%         Second char: type - one of
%                      'c' GSHHS coastline (default)
%                      'b' WDB Border
%                      'r' WDB River
%
%         Third char - if 2nd char is 'b':
%                      '1' Country borders
%                      '2' State/Province and Country borders
%                    - if 2nd char is 'r': '1','2','3','4'
%                      add successively more tributaries
%
%         (also maintained is this optional format:
%
%         RES - selections resolution
%                  1  or 'crude'
%                  2  or 'low'
%                  3  or 'intermediate'
%                  4  or 'high'
%                  5  or 'full
%
%         but please don't use this).
%
%         See also M_PROJ, M_GRID, M_COAST, M_GSHHS_L, M_GSHHS_H, M_GSHHS_C
%         M_USERCOAST

% Rich Pawlowicz (rich@ocgy.ubc.ca) 15/June/98
%
%
% This software is provided "as is" without warranty of any kind. But
% it's mine, so you can't sell it.
%
%  16/Dec/2005
%*********************************************************************
%  Modified after code provided by Bruce Lipphardt (brucel@udel.edu) to
%  reduce the hierarchy of M_GSHHS_* routines to a single routine with a
%  variable resolution input:
% 20/Jan/2008 - added borders and rivers from gshhs v1.10
% 4/DEc/11 - isstr to ischar
% Sep/14 - added hierarchy to borders
% Aug/18 - fixed error that occurred when called m_gshhs_X with 'save'
%          option (Thanks to H. Grant for pointing this out).

% Root of directories where all the gshhs_X.b, wdb_borders-X.b and wdb_rivers_X.b
% files live
FILNAME='private/';

%-------------don't change below here----------------------------

res_list = {'c','l','i','h','f'};
typ_list={'c','b','r'};
typ_names={'gshhs_','wdb_borders_','wdb_rivers_'};

typ=1;
flaglim='9';

if ischar(resolution)
 if length(resolution)>=2
   typ = find(strcmpi(resolution(2),typ_list));
 end
 if length(resolution)>=3
   flaglim = resolution(3);
 end
 resolution =  find(strcmpi(resolution(1),res_list));
end

if isempty(resolution) || resolution<1 || resolution> length(res_list)
  error('**Don''t recognize the specified resolution');
end
if isempty(typ) || typ<1 || typ> length(res_list)
  error('**Don''t recognize the specified type');
end

res_char = res_list{resolution} ;
file     = [FILNAME,sprintf('%s%s.b',typ_names{typ},res_char)] ;
tag_name = sprintf('m_%s%s',typ_names{typ},res_char) ;

% Set current projection to geographic
Currentmap=m_coord('set');
m_coord('geographic');

if length(varargin)>1 && strcmp(varargin{1},'save')
  [ncst,Area,k]=mu_coast(res_char,file);
  save(varargin{2},'ncst','k','Area');
  h=varargin{2};   % Error if you call m_gshhs_i with 'save' option - thanks HG, Aug/1/2018
else
  h=mu_coast([res_char flaglim],file,varargin{:},'tag',tag_name);
end

m_coord(Currentmap.name);

if nargout==0
    clear h
end

author: goGPS-Project


text:
function [xout,yout] = coords2normfig(x,y,h)

if ~exist('h','var')
  h = gca;
end

if ~strcmpi(get(h,'type'),'axes'),
  error('Handle must correspond to an axis');
end

xax  = get(h,'xlim');
yax = get(h,'ylim');

pos = get(h,'position');

xout = (x - pos(1))/pos(3)*(xax(2)-xax(1)+1) + xax(1);
yout = (y - pos(2))/pos(4)*(yax(2)-yax(1)+1) + yax(1);

yisreversed = strcmpi(get(h,'ydir'),'reverse');
xisreversed = strcmpi(get(h,'xdir'),'reverse');
if xisreversed,
  xout = (1-(x - xax(1))/(xax(2)-xax(1)))*pos(3) + pos(1);
else
  xout = (x - xax(1))/(xax(2)-xax(1))*pos(3) + pos(1);
end
if yisreversed,
  yout = (1-(y - yax(1))/(yax(2)-yax(1)))*pos(4) + pos(2);
else
  yout = (y - yax(1))/(yax(2)-yax(1))*pos(4) + pos(2);
end

author: kristinbranson

repo: JAABA


text:
function theta = vectorAngle3d(v1, v2)
%VECTORANGLE3D Angle between two 3D vectors.
%
%   THETA = vectorAngle3d(V1, V2)
%   Computes the angle between the 2 3D vectors V1 and V2. The result THETA
%   is given in radians, between 0 and PI.
%
%
%   Example
%   % angle between 2 orthogonal vectors
%   vectorAngle3d([1 0 0], [0 1 0])
%   ans =
%       1.5708
%
%   % angle between 2 parallel vectors
%   v0 = [3 4 5];
%   vectorAngle3d(3*v0, 5*v0)
%   ans =
%       0
%
%   See also
%   vectors3d, vectorNorm3d, crossProduct3d
%

% ------
% Author: David Legland
% E-mail: david.legland@inra.fr
% Created: 2010-10-04, using Matlab 7.9.0.529 (R2009b)
% Copyright 2010-2022 INRA - Cepia Software Platform

% compute angle using arc-tangent to get better precision for angles near
% zero, see the discussion in:
% http://www.mathworks.com/matlabcentral/newsreader/view_thread/151925#381952
theta = atan2(vectorNorm3d(crossProduct3d(v1, v2)), sum(bsxfun(@times, v1, v2),2));

% equivalent to:
% v1 = normalizeVector3d(v1);
% v2 = normalizeVector3d(v2);
% theta = acos(dot(v1, v2, 2));

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function vcell = lift_ape_v1(r,t,theta,tBound,depthBound,FOV)
    r       = r(:);
    t       = t(:);
    x       = [r;t];
    theta   = theta(:);
    v       = [1;x;theta;kron(theta,x)];
    if tBound^2 < t'*t
        v1  = [1;theta] * 0;
    else
        v1  = [1;theta] * sqrt(tBound^2 - t'*t);
    end

    if t(end) < depthBound
        v2  = [1;theta] * 0;
    else
        v2  = [1;theta] * sqrt(t(end)-depthBound);
    end

    scale = tan(deg2rad(FOV)/2);
    if scale^2 * t(3)^2 < t(1)^2 + t(2)^2
        v3  = [1;theta] * 0;
    else
        v3  = [1;theta] * sqrt(scale^2 * t(3)^2 - t(1)^2 - t(2)^2);
    end

    vcell   = {v;v1;v2;v3};
end

author: MIT-SPARK

repo: CertifiablyRobustPerception


text:
% Script used to generate a visualisation of how well EKF, UKF and PF
% approximate a distribution
% =========================================================================

% Assume a target positioned at x = 1, travelling with speed v = 0.1
state = [2000; 0.1; 0; 0];

% Instantiate an Observation model
% obs = LinGaussObsModelX_2D('NumStateDims',4,'ObsErrVariance',0.2,'Mapping',[1 2]);
obs = RangeBearing2CartesianX('NumStateDims',4,'MeasurementErrVariance',[(pi/90)^2,(1)^2],...
                              'Mapping',[1 3]);

% View the transition matrix and process covariance matrices
R = obs.covar();

% Predict the target's position and velocity after the interval has passed
measurement  = obs.feval(state);

%% EKF
[ekf_mean, H_jac] = ExtendedKalmanFilterX.computeJac_(@(x)obs.finv(x),measurement);
ekf_covar = H_jac*R*H_jac';

%% UKF
alpha = 0.5;
kappa = 0;
beta  = 2;
% Calculate unscented transformation parameters
[c, Wmean, Wcov, OOM] = matlabshared.tracking.internal.calcUTParameters(alpha,beta,kappa,2);
% Form the sigma points
X = formSigmaPoints(measurement, R, c);
% Perform Unscented Transform to get predicted measurement mean,
% covariance and cross-covariance
[ukf_mean,ukf_covar,Pxy] = unscentedTransform(@(x)obs.finv(x),X,Wmean,Wcov,OOM);

%% PF
% Generate 50 random noise samples from the dynamic model
noise = obs.random(50000000);

% Add noise to the measurements
Yk1 = noise + measurement;
Xk1 = obs.finv(Yk1);

Xk = [Xk1];%,Xk2];
Yk = [Yk1];%,Yk2];

figure;
[bandwidth,density,X,Y]=kde2d(Xk([1,3],:)');
xLim = [1985,2005];
yLim = [-200,200];

% EKF
subplot(1,3,1);
%contour(X,Y,density,1);
h = surf(X,Y,density);
z = get(h,'ZData');
% set(h,'ZData',z-10^6)
hold on;
shading interp
x1 = xLim(1):0.1:xLim(2); x2 = yLim(1):0.1:yLim(2);
[X1,X2] = meshgrid(x1,x2);
F = mvnpdf([X1(:) X2(:)],ekf_mean([1,3])',ekf_covar([1,3],[1,3]));
F = reshape(F,length(x2),length(x1));
contour(x1,x2,F,4,'r','LineWidth',2);
% plot_gaussian_ellipsoid(ekf_mean([1,3]),ekf_covar([1,3],[1,3]),'r',1);
% plot_gaussian_ellipsoid(ekf_mean([1,3]),ekf_covar([1,3],[1,3]),'r',2);
% shading interp
%colormap(jet(3000))
title('Taylor-Series Expansion (EKF)')
xlabel('X (m)');
ylabel('Y (m)');
xlim(xLim);
ylim(yLim);
view(0,-90);

% UKF
subplot(1,3,2);
h = surf(X,Y,density);
z = get(h,'ZData');
%set(h,'ZData',z-10^6)
hold on;
shading interp
x1 = xLim(1):0.1:xLim(2); x2 = yLim(1):0.1:yLim(2);
[X1,X2] = meshgrid(x1,x2);
F = mvnpdf([X1(:) X2(:)],ukf_mean([1,3])',ukf_covar([1,3],[1,3]));
F = reshape(F,length(x2),length(x1));
contour(x1,x2,F,4,'r','LineWidth',2);
%plot_gaussian_ellipsoid(ukf_mean([1,3]),ukf_covar([1,3],[1,3]),'r',1);
%plot_gaussian_ellipsoid(ukf_mean([1,3]),ukf_covar([1,3],[1,3]),'r',2);
%colormap(jet(3000))
title('Unscented Transform (UKF)')
xlabel('X (m)');
ylabel('Y (m)');
xlim(xLim);
ylim(yLim);
view(0,-90);

% PF
subplot(1,3,3);
%contour(X,Y,density,1);
h = surf(X,Y,density);
z = get(h,'ZData');
%set(h,'ZData',z-10^6)
hold on;
shading interp
%plot(Xk(1,1:1000),Xk(3,1:1000),'r.');
%[bandwidth,density,X,Y]=kde2d(Xk([1,3],1:10000)');
contour(X,Y,density,4,'r','LineWidth',2)
%colormap(jet(3000))
title('Sampling (PF)')
%set(h,'ZData',z-10)
xlabel('X (m)');
ylabel('Y (m)');
%set(gca,'Xdir','reverse')
set(gca,'Ydir','reverse')
xlim(xLim);
ylim(yLim);
view(0,-90);

% model = PositionalObsModelX(config);%PositionalObsModelX(config);%Polar2CartGaussianModelX(config);
%
% xkm1 = [0.5; 0.3; sqrt(2); sqrt(2)];
% pkm1 = [0 0.5; 0 0.3; sqrt(2) sqrt(2); sqrt(2) sqrt(2)];
% Pkm1 = [0 1 0 1; 0 1 0 1; 0 1 0 1; 1 1 0 1];
% yk = model.obs(1,xkm1);
% Yk = model.sample(1, yk, 500000);
% pk = model.obs_cov();
% Pk = model.eval(1,Yk,xkm1);
% [bandwidth,density,X,Y]=kde2d(Yk');
% figure;
% %contour3(X,Y,density,50);
% surf(X,Y,density);
% hold on;

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function model = Initialization_Cluster_Basis(Seqs, ClusterNum, ...
                        baseType, bandwidth, landmark)

N = length(Seqs);
D = zeros(N,1);
for i = 1:N
    D(i) = max(Seqs(i).Mark);
end
D = max(D);
model.K = ClusterNum;
model.D = D;

switch nargin
    case 2
        sigma = zeros(length(Seqs), 1);
        Tmax = zeros(length(Seqs), 1);

        for i = 1:length(Seqs)
            sigma(i) = ((4*std(Seqs(i).Time)^5)/(3*length(Seqs(i).Time)))^0.2;
            Tmax(i) = Seqs(i).Time(end) + eps;
        end
        Tmax = mean(Tmax);

        model.kernel = 'gauss';
        model.w = mean(sigma);
        model.landmark = model.w*(0:ceil(Tmax/model.w));

    case 3
        model.kernel = baseType;
        model.w = 1;
        model.landmark = 0;
    case 4
        model.kernel = baseType;
        model.w = bandwidth;
        model.landmark = 0;
    otherwise
        model.kernel = baseType;
        model.w = bandwidth;
        model.landmark = landmark;
end

model.alpha = 1;
M = length(model.landmark);
model.beta = ones(D, M, model.K, D)./(M*D^2);
% hyperparameter of Rayleigh prior of basic intensity (upper bound)
model.b = ones(D, model.K)./(D);

% initialize label and responsibility randomly
label = ceil(model.K*rand(1,N));
model.R = full(sparse(1:N,label,1,N,model.K,N));

author: HongtengXu

repo: Hawkes-Process-Toolkit


text:
function value = p10_i4 ( action, name, value )

%*****************************************************************************80
%
%% P10_I4 sets or gets I4 parameters for problem 10.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string ACTION,
%    'D' to set a parameter to its default value.
%    'S' to set a parameter,
%    'G' to get a parameter,
%    'I' to increment a parameter.
%
%    Input, string NAME, the name of the variable.
%    '#' is the number of calls to the integrand routine.
%
%    Input/output, integer VALUE.
%    * If ACTION = 'I', then VALUE is an input quantity, and is the
%      new value to be added to NAME.
%    * If ACTION = 'S', then VALUE is an input quantity, and is the
%      new value to be assigned to NAME.
%    * If ACTION = 'G', then VALUE is an output quantity, and is the
%      current value of NAME.
%
  persistent calls;

  if ( action == 'D' | action == 'd' )

    if ( name == '#' | name == '*' )
      calls = 0;
    end

  elseif ( action == 'G' | action == 'g' )

    if ( name == '#' )
      value = calls;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P10_I4 - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized name = "%s".\n', name );
      error ( 'P10_I4 - Fatal error!' );
    end

  elseif ( action == 'I' | action == 'i' )

    if ( name == '#' )
      calls = calls + value;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P10_I4 - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized name = "%s".\n', name );
      error ( 'P10_I4 - Fatal error!' );
    end

  elseif ( action == 'S' | action == 's' )

    if ( name == '#' )
      calls = value;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P10_I4 - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized name = "%s".\n', name );
      error ( 'P10_I4 - Fatal error!' );
    end

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'P10_I4 - Fatal error!\n' );
    fprintf ( 1, '  Unrecognized action = "%s".\n', action );
    error ( 'P10_I4 - Fatal error!' );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% FLOATREAD - Read matrix from float file ssuming four byte floating point number
%               Can use FSEEK to read an arbitrary (contiguous) submatrix.
%
% Usage:        >> a = floatread(filename,size,'format',offset)
%
% Inputs:
%   filename - name of the file
%   size     - determine the number of float elements to be read and
%              the dimensions of the resulting matrix. If the last element
%              of 'size' is Inf, the size of the last dimension is determined
%              by the file length. If size is 'square,' FLOATREAD attempts
%              to read a square 2-D matrix.
%
% Optional inputs:
%  'format'  - the option FORMAT argument specifies the storage format as
%              defined by FOPEN. Default format ([]) is 'native'.
%  offset    - either the number of first floats to skip from the beginning of the
%              float file, OR a cell array containing the dimensions of the original
%              data matrix and a starting position vector in that data matrix.
%
%              Example: % Read a [3 10] submatrix of a four-dimensional float matrix
%                >> a = floatread('mydata.fdt',[3 10],'native',{[[3 10 4 5],[1,1,3,4]});
%              % Note: The 'size' and 'offset' arguments must be compatible both
%              % with each other and with the size and ordering of the float file.
%
% Author: Sigurd Enghoff, CNL / Salk Institute, La Jolla, 7/1998
%
% See also: FLOATWRITE, FOPEN

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

% 04-26-99  modified by Sigurd Enghoff to handle variable-sized and
%           multi-dimensional data.
% 07-08-99  modified by Sigurd Enghoff, FORMAT argument added.
% 02-08-00  help updated for toolbox inclusion -sm
% 02-14-00  added segment arg -sm
% 08-14-00  added size 'square' option -sm
% 01-25-02  reformated help & license, added links -ad

function A = floatread(fname,Asize,fform,offset)

if nargin<2
  help floatread
  return
end

if ~exist('fform') || isempty(fform) || fform(1)==0
	fform = 'native';
end

if ~exist('offset')
	offset = 0;
end

fid = fopen(fname,'rb',fform);
if fid>0
 if exist('offset')
   if iscell(offset)
     if length(offset) ~= 2
        error('offset must be a positive integer or a 2-item cell array');
     end
     datasize = offset{1};
     startpos = offset{2};
     if length(datasize) ~= length(startpos)
        error('offset must be a positive integer or a 2-item cell array');
     end
     for k=1:length(datasize)
       if startpos(k) < 1 || startpos(k) > datasize(k)
          error('offset must be a positive integer or a 2-item cell array');
       end
     end
     if length(Asize)> length(datasize)
        error('offset must be a positive integer or a 2-item cell array');
     end
     for k=1:length(Asize)-1
         if startpos(k) ~= 1
            error('offset must be a positive integer or a 2-item cell array');
         end
     end
     sizedim = length(Asize);
     if Asize(sizedim) + startpos(sizedim) - 1 > datasize(sizedim)
        error('offset must be a positive integer or a 2-item cell array');
     end
     for k=1:length(Asize)-1
         if Asize(k) ~= datasize(k)
            error('offset must be a positive integer or a 2-item cell array');
         end
     end

     offset = 0;
     jumpfac = 1;
     for k=1:length(startpos)
           offset = offset + jumpfac * (startpos(k)-1);
           jumpfac = jumpfac * datasize(k);
     end

   elseif length(offset) > 1
     error('offset must be a positive integer or a 2-item cell array');
   end

   % perform the fseek() operation
   % -----------------------------
   stts = fseek(fid,4*offset,'bof');

   if stts ~= 0
     error('floatread(): fseek() error.');
     return
   end
 end

% determine what 'square' means
% -----------------------------
 if ischar('Asize')
   if iscell(offset)
         if length(datasize) ~= 2 || datasize(1) ~= datasize(2)
              error('size ''square'' must refer to a square 2-D matrix');
         end
         Asize = [datsize(1) datasize(2)];
   elseif strcmp(Asize,'square')
         fseek(fid,0,'eof'); % go to end of file
         bytes = ftell(fid); % get byte position
         fseek(fid,0,'bof'); % rewind
         bytes = bytes/4; % nfloats
         froot = sqrt(bytes);
         if round(froot)*round(froot) ~= bytes
              error('floatread(): filelength is not square.')
         else
              Asize = [round(froot) round(froot)];
         end
   end
 end
 A = fread(fid,prod(Asize),'float');
else
 error('floatread() fopen() error.');
 return
end

% fprintf('   %d floats read\n',prod(size(A)));

% interpret last element of Asize if 'Inf'
% ----------------------------------------
if Asize(end) == Inf
	Asize = Asize(1:end-1);
	A = reshape(A,[Asize length(A)/prod(Asize)]);
else
	A = reshape(A,Asize);
end

fclose(fid);

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function PlotInterElem(vert,intpt)

scatter3(vert(1,1),vert(1,2),vert(1,3),'k');
hold on

scatter3(vert(2,1),vert(2,2),vert(2,3),'k');
scatter3(vert(3,1),vert(3,2),vert(3,3),'k');
scatter3(vert(4,1),vert(4,2),vert(4,3),'k');

if size(intpt,1) == 3
    scatter3(intpt(1,1),intpt(1,2),intpt(1,3),'r');
    scatter3(intpt(2,1),intpt(2,2),intpt(2,3),'b');
    scatter3(intpt(3,1),intpt(3,2),intpt(3,3),'m');
elseif size(intpt,1) == 4
    scatter3(intpt(1,1),intpt(1,2),intpt(1,3),'k');
    scatter3(intpt(2,1),intpt(2,2),intpt(2,3),'b');
    scatter3(intpt(3,1),intpt(3,2),intpt(3,3),'k');
    scatter3(intpt(4,1),intpt(4,2),intpt(4,3),'y');
end

axis equal

hold off

repo: ifem


text:
function [ce_mean_p, ce_mean_ps, ce_mean_s, ce_mean_sn, ce_mean_n] = interpolateElectrolyteConcentration(ce,param)
%	interpolateElectrolyteConcentration interpolates the value of electrolyte concentration at the edges of control volumes using harmonic mean.

%   This file is part of the LIONSIMBA Toolbox
%
%	Official web-site: 	http://sisdin.unipv.it/labsisdin/lionsimba.php
% 	Official GitHUB: 	https://github.com/lionsimbatoolbox/LIONSIMBA
%
%   LIONSIMBA: A Matlab framework based on a finite volume model suitable for Li-ion battery design, simulation, and control
%   Copyright (C) 2016-2018 :Marcello Torchio, Lalo Magni, Davide Raimondo,
%                            University of Pavia, 27100, Pavia, Italy
%                            Bhushan Gopaluni, Univ. of British Columbia,
%                            Vancouver, BC V6T 1Z3, Canada
%                            Richard D. Braatz,
%                            Massachusetts Institute of Technology,
%                            Cambridge, Massachusetts 02142, USA
%
%   Main code contributors to LIONSIMBA 2.0:
%                           Ian Campbell, Krishnakumar Gopalakrishnan,
%                           Imperial college London, London, UK
%
%   LIONSIMBA is a free Matlab-based software distributed with an MIT
%   license.

%% Electrolyte concentration interpolation

% Interpolation within the positive electrode
beta_ce_p = 0.5;
ce_mean_p = ce(1:param.Np-1).*ce(2:param.Np)./ (beta_ce_p*ce(2:param.Np) + (1-beta_ce_p)*ce(1:param.Np-1));

% Interpolation on the interface between separator and positive electrode
beta_ce_ps = param.deltax_p*param.len_p/2 / (param.deltax_p*param.len_p/2 + param.deltax_s*param.len_s/2);
ce_mean_ps = ce(param.Np)*ce(param.Np+1)/(beta_ce_ps*ce(param.Np+1) + (1-beta_ce_ps)*ce(param.Np));

% Interpolation within the separator
beta_ce_s = 0.5;
ce_mean_s = ce(param.Np+1:param.Np+param.Ns-1).*ce(param.Np+2:param.Np+param.Ns)./ (beta_ce_s*ce(param.Np+2:param.Np+param.Ns) + (1-beta_ce_s)*ce(param.Np+1:param.Np+param.Ns-1));

% Interpolation on the interface between separator and negative electrode
beta_ce_sn = param.deltax_s*param.len_s/2 / (param.deltax_n*param.len_n/2 + param.deltax_s*param.len_s/2);
ce_mean_sn = ce(param.Np+param.Ns)*ce(param.Np+param.Ns+1)/(beta_ce_sn*ce(param.Np+param.Ns+1) + (1-beta_ce_sn)*ce(param.Np+param.Ns));

% Interpolation within the negative electrode
beta_ce_n = 0.5;
ce_mean_n = ce(param.Np+param.Ns+1:end-1).*ce(param.Np+param.Ns+2:end)./ (beta_ce_n*ce(param.Np+param.Ns+2:end) + (1-beta_ce_n)*ce(param.Np+param.Ns+1:end-1));

end

author: lionsimbatoolbox

repo: LIONSIMBA


text:
function [phase,amplitude,unwrapped] = Phase(samples,times)

%Phase - Compute instantaneous phase in signal.
%
% Compute instantaneous phase (wrapped and unwrapped) and amplitude
% in signal.
%
%  USAGE
%
%    [phase,amplitude,unwrapped] = Phase(samples,times)
%
%    samples        signal, e.g. filtered local field potential samples
%    times          optional timestamps where phase should be interpolated
%
%  NOTE
%
%    Angles are returned in radians.
%
%  SEE
%
%    See also FilterLFP, PhasePrecession, PhaseMap.

% Copyright (C) 2004-2011 by Michaël Zugaro
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.

% Check number of parameters
if nargin < 1,
	error('Incorrect number of parameters (type ''help Phase'' for details).');
end

% Check parameter sizes
if size(samples,2) < 2,
	error('Parameter ''samples'' is not a matrix (type ''help Phase'' for details).');
end
if nargin == 2,
	if ~isdvector(times),
		error('Parameter ''times'' is not a vector (type ''help Phase'' for details).');
	end
	if size(times,2) ~= 1, times = times'; end
end

phase = zeros(size(samples));
phase(:,1) = samples(:,1);
unwrapped = zeros(size(samples));
unwrapped(:,1) = samples(:,1);
amplitude = zeros(size(samples));
amplitude(:,1) = samples(:,1);

for j = 2:size(samples,2),
	% Compute phase and amplitude using Hilbert transform
	h = hilbert(samples(:,j));
	phase(:,j) = mod(angle(h),2*pi);
	amplitude(:,j) = abs(h);
	unwrapped(:,j) = unwrap(phase(:,j));
end

if nargin == 2,
%  	phase = Interpolate([phase(:,1) exp(i*phase(:,2))],times(:,1));
%  	phase(:,2) = angle(phase(:,2));
	phase = Interpolate(phase,times(:,1),'type','circular');
	amplitude = Interpolate(amplitude,times(:,1));
	unwrapped = Interpolate(unwrapped,times(:,1));
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function tetra_node2 = tet_mesh_order10_to_order4_compute ( tetra_num1, ...
  tetra_node1, tetra_num2 )

%*****************************************************************************80
%
%% TET_MESH_ORDER10_TO_ORDER4_COMPUTE linearizes a quadratic tet mesh.
%
%  Discussion:
%
%    A quadratic tet mesh is assumed to consist of 10-node
%    tetrahedrons.
%
%    This routine rearranges the information so as to define a 4-node
%    tet mesh.
%
%    The same nodes are used, but there are 8 times as many
%    tetrahedrons.
%
%    The node ordering for the quadratic tetrahedron is somewhat
%    arbitrary.  In the current scheme, the vertices are listed
%    first, followed by the 6 midside nodes.  Each midside node
%    may be identified by the two vertices that bracket it.  Thus,
%    the node ordering may be suggested by:
%
%      1  2  3  4 (1+2) (1+3) (1+4) (2+3) (2+4) (3+4)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 January 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Anwei Liu, Barry Joe,
%    Quality Local Refinement of Tetrahedral Meshes Based
%    on 8-Subtetrahedron Subdivision,
%    Mathematics of Computation,
%    Volume 65, Number 215, July 1996, pages 1183-1200.
%
%  Parameters:
%
%    Input, integer TETRA_NUM1, the number of tetrahedrons in the quadratic
%    tet mesh.
%
%    Input, integer TETRA_NODE1(10,TETRA_NUM1), the quadratic tet mesh.
%
%    Input, integer TETRA_NUM2, the number of tetrahedrons in the linear
%    tet mesh.  TETRA_NUM2 = 8 * TETRA_NUM1.
%
%    Output, integer TETRA_NODE2(4,TETRA_NUM2), the linear tet mesh.
%
  tetra2 = 0;

  for tetra1 = 1 : tetra_num1

    n1 = tetra_node1(1,tetra1);
    n2 = tetra_node1(2,tetra1);
    n3 = tetra_node1(3,tetra1);
    n4 = tetra_node1(4,tetra1);
    n5 = tetra_node1(5,tetra1);
    n6 = tetra_node1(6,tetra1);
    n7 = tetra_node1(7,tetra1);
    n8 = tetra_node1(8,tetra1);
    n9 = tetra_node1(9,tetra1);
    nx = tetra_node1(10,tetra1);

    tetra2 = tetra2 + 1;
    tetra_node2(1:4,tetra2 ) = [ n1, n5, n6, n7 ]';
    tetra2 = tetra2 + 1;
    tetra_node2(1:4,tetra2 ) = [ n2, n5, n8, n9 ]';
    tetra2 = tetra2 + 1;
    tetra_node2(1:4,tetra2 ) = [ n3, n6, n8, n9 ]';
    tetra2 = tetra2 + 1;
    tetra_node2(1:4,tetra2 ) = [ n4, n7, n9, nx ]';
    tetra2 = tetra2 + 1;
    tetra_node2(1:4,tetra2)  = [ n5, n6, n7, n9 ]';
    tetra2 = tetra2 + 1;
    tetra_node2(1:4,tetra2 ) = [ n5, n6, n8, n9 ]';
    tetra2 = tetra2 + 1;
    tetra_node2(1:4,tetra2 ) = [ n6, n7, n9, nx ]';
    tetra2 = tetra2 + 1;
    tetra_node2(1:4,tetra2 ) = [ n6, n8, n9, nx ]';

  end

  return
end


text:
function [p alpha] = circ_vmpdf(alpha, thetahat, kappa)

% [p alpha] = circ_vmpdf(alpha, w, p)
%   Computes the circular von Mises pdf with preferred direction thetahat
%   and concentration kappa at each of the angles in alpha
%
%   The vmpdf is given by f(phi) =
%   (1/(2pi*I0(kappa))*exp(kappa*cos(phi-thetahat)
%
%   Input:
%     alpha     angles to evaluate pdf at, if empty alphas are chosen to
%               100 uniformly spaced points around the circle
%     [thetahat preferred direction, default is 0]
%     [kappa    concentration parameter, default is 1]
%
%   Output:
%     p         von Mises pdf evaluated at alpha
%     alpha     angles at which pdf was evaluated
%
%
%   References:
%     Statistical analysis of circular data, Fisher
%
% Circular Statistics Toolbox for Matlab

% By Philipp Berens and Marc J. Velasco, 2009
% velasco@ccs.fau.edu

% if no angles are supplied, 100 evenly spaced points around the circle are
% chosen
if nargin < 1 || isempty(alpha)
    alpha = linspace(0, 2*pi, 101)';
    alpha = alpha(1:end-1);
end
if nargin < 3
    kappa = 1;
end
if nargin < 2
    thetahat = 0;
end

alpha = alpha(:);

% evaluate pdf
C = 1/(2*pi*besseli(0,kappa));
p = C * exp(kappa*cos(alpha-thetahat));