text:
function [xi,w]=sixthOrderTriangleCubPoints()
%%SIXTHORDERTRIANGLECUBPOINTS Obtain sixth-order cubature points for
%   integration over a triangle in 2D. The points and weights are for the
%   triangle with vertices (1,0), (0,1), (0,0), but can be transformed to
%   any triangle using transformSimplexTriPoints.
%
%INPUTS: None
%
%OUTPUTS: xi A 2XnumCubPoints set of points for the standard triangle.
%          w A 1XnumCubPoints set of cubature weights. This sums to the
%            volume of the triangle (1/2).
%
%This function implements the points given in [1] (12 points).
%
%EXAMPLE:
%Given the vertices of the simplex, we compare a sixth-order moment
%computed using these cubature points to one computed using
%monomialIntSimplex. The results are the same within typical finite
%precision limits.
% [xi,w]=sixthOrderTriangleCubPoints();
% alpha=[4;2];
% theMoment=findMomentFromSamp(alpha,xi,w)
% intVal=monomialIntSimplex(alpha)
%
%REFERENCES:
%[1] F. D. Witherden and P. E. Vincent, "On the identification of symmetric
%    quadrature rules for finite element methods," Computer and Mathematics
%    with Applications, vol. 69, no. 10, pp. 1232-1241, May 2015.
%
%October 2022 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

M=[-0.87382197101699554331933679425836168532,    0.74764394203399108663867358851672337064,    0.10168981274041363384187361821373796809;
    0.74764394203399108663867358851672337064,   -0.87382197101699554331933679425836168532,    0.10168981274041363384187361821373796809;
   -0.87382197101699554331933679425836168532,   -0.87382197101699554331933679425836168532,    0.10168981274041363384187361821373796809;
   -0.50142650965817915741672289378596184782,  0.0028530193163583148334457875719236956481,    0.23357255145275873205057922277115888265;
  0.0028530193163583148334457875719236956481,   -0.50142650965817915741672289378596184782,    0.23357255145275873205057922277115888265;
   -0.50142650965817915741672289378596184782,   -0.50142650965817915741672289378596184782,    0.23357255145275873205057922277115888265;
   -0.89370990031036610529350065673720370601,    0.27300499824279729446028518882409939961,    0.16570215123674715038710691284088490796;
    0.27300499824279729446028518882409939961,   -0.89370990031036610529350065673720370601,    0.16570215123674715038710691284088490796;
   -0.37929509793243118916678453208689569359,    0.27300499824279729446028518882409939961,    0.16570215123674715038710691284088490796;
    0.27300499824279729446028518882409939961,   -0.37929509793243118916678453208689569359,    0.16570215123674715038710691284088490796;
   -0.37929509793243118916678453208689569359,   -0.89370990031036610529350065673720370601,    0.16570215123674715038710691284088490796;
   -0.89370990031036610529350065673720370601,   -0.37929509793243118916678453208689569359,    0.16570215123674715038710691284088490796];
w=M(:,3);
xi=M(:,1:2)';
%Transform the points to the standard triangle.
v1=[-1,-1, 1;
    -1, 1,-1];
v2=[1,0,0;
    0,1,0];
[A,d]=affineTransBetweenTriangles(v1,v2);
xi=bsxfun(@plus,A*xi,d);
w=w/4;
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function wavelet_test03 ( )

%*****************************************************************************80
%
%% WAVELET_TEST03 tests DAUB6_TRANSFORM and DAUB6_TRANSFORM_INVERSE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 July 2011
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'WAVELET_TEST03\n' );
  fprintf ( 1, '  DAUB6_TRANSFORM computes the DAUB6 transform of a vector.\n' );
  fprintf ( 1, '  DAUB6_TRANSFORM_INVERSE inverts it.\n' );
%
%  Random data.
%
  n = 16;
  seed = 123456789;
  [ u, seed ] = r8vec_uniform_01 ( n, seed );

  v = daub6_transform ( n, u );

  w = daub6_transform_inverse ( n, v );

  fprintf ( 1, '\n' );
  fprintf ( 1, '   i      U(i)       D6(U)(i)  D6inv(D6(U))(i)\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %2d  %10.4f  %10.4f  %10.4f\n', i, u(i), v(i), w(i) );
  end
%
%  Constant signal.
%
  n = 8;
  u(1:n) = 1.0;

  v = daub6_transform ( n, u );

  w = daub6_transform_inverse ( n, v );

  fprintf ( 1, '\n' );
  fprintf ( 1, '   i      U(i)        D6(U)(i)  D6inv(D6(U))(i)\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %2d  %10.4f  %10.4f  %10.4f\n', i, u(i), v(i), w(i) );
  end
%
%  Linear signal.
%
  n = 16;
  a_first = 1.0;
  a_last = n;
  u = linspace ( a_first, a_last, n );

  v = daub6_transform ( n, u );

  w = daub6_transform_inverse ( n, v );

  fprintf ( 1, '\n' );
  fprintf ( 1, '   i      U(i)        D6(U)(i)  D6inv(D6(U))(i)\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %2d  %10.4f  %10.4f  %10.4f\n', i, u(i), v(i), w(i) );
  end
%
%  Quadratic data.
%
  n = 8;
  u(1) = 25.0;
  u(2) = 16.0;
  u(3) = 9.0;
  u(4) = 4.0;
  u(5) = 1.0;
  u(6) = 0.0;
  u(7) = 1.0;
  u(8) = 4.0;

  v = daub6_transform ( n, u );

  w = daub6_transform_inverse ( n, v );

  fprintf ( 1, '\n' );
  fprintf ( 1, '   i      U(i)        D6(U)(i)  D6inv(D6(U))(i)\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %2d  %10.4f  %10.4f  %10.4f\n', i, u(i), v(i), w(i) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [infos] = canlab_glm_getinfo(modeldir,varargin)
% :SUBJECT LEVEL input:
% ::
%
%    INFO = canlab_glm_getinfo(spm_subject_level_directory, option, [n])
%
% Get information out of an spm subject level analysis.
%
% :Options: (each option can be called by the listed letter or word + a number, when noted)
%
%   **i' 'input':**
%        number of volumes and (first) volume name for each run
%
%   **'b' 'betas' [n]:**
%        beta names (for nth session)
%
%   **'B' 'taskbetas' [n]:**
%        beta names (that didn't come from multiple regressors)  (for nth session)
%
%   **'c' 'cons' [n]:**
%        contrast names (for nth contrast)
%
%   **'C' 'conw' [n]:**
%        beta names and weights for contrasts (for nth con)
%
%   **'v' 'image' [n]:**
%        create figure of design matrix (for nth session)
%                            (design matrix is multiplied by 100 for visibility)
%                            (works well for multiple runs)
%
%   **'V' 'taskimage' [n]:**
%        same as 'image', but only for task betas
%
%   **'imagesc':**
%        same as 'image', but uses imagesc
%                            (works well for single runs)
%
%
% :GROUP LEVEL input:
% ::
%
%    INFO = canlab_glm_getinfo(robfit_group_level_directory, option, [n])
%
% Get information out of a robfit group level analysis.
%
% :Options: (each option can be called by the listed word or letter + a number, when noted)
%
%   Any of the subject level options can be used on a group level robfit
%   analysis by prefixing '1i' (output is generated based on the first input
%   to the first robust analysis).
%
%   Ex:
%   ::
%
%      canlab_glm_getinfo('second_level/model3','1iconw')
%
%
%   **'i' 'input' [n]:**
%        input contrasts by number and name (for nth analysis)
%
%   **'I' 'allinput' [n]:**
%        input images (for nth analysis)
%
%   **'m' 'model':**
%        weights by subject (i.e., directory containing input contrast images)
%
%   **'M' 'allmodels' [n]:**
%        weights and input images (for nth analysis)
%
% :Assumptions:
%    In some options, the first contrasts and group level analysis
%    directories are assumed to represent the rest, which may not be the
%    case.
%
% :Note:
%    group level options do not yet return a usable INFO struct.

if nargin ~= 2 && nargin ~= 3
    error('USAGE: INFO = canlab_glm_getinfo(spm_dir|robfit_dir,option,[n])')
end

opt = varargin{1};
if nargin == 3, s = varargin{2}; else s = 0; end

if ~exist(modeldir,'dir')
    error('No such directory: %s',modeldir)
else
    spmmat = fullfile(modeldir,'SPM.mat');
    setupmats = filenames(fullfile(modeldir,'robust[0-9][0-9][0-9][0-9]','SETUP.mat'));
    if exist(spmmat,'file')
        infos = get_subject_level_info(modeldir,opt,s);
    elseif ~isempty(setupmats)
        if regexp(opt,'^1i')
            load(setupmats{1});
            firstinput = fileparts(deblank(SETUP.files(1,:)));
            infos = get_subject_level_info(firstinput,regexprep(opt,'^1i',''),s);
        else
            infos = get_group_level_info(modeldir,opt,s);
        end
    else
        error('%s is neither a subject level SPM directory nor a group level robfit directory',modeldir);
    end
end

end

function [infos] = get_subject_level_info(modeldir,opt,s)

load(fullfile(modeldir,'SPM'));

switch opt
    case {'c' 'con' 'cons'}
        if s
            infos.contrast_name{1} = SPM.xCon(s).name;
        else
            infos.contrast_name = cellstr(strvcat(SPM.xCon.name));  %#ok
        end

        fprintf('CONTRASTS:\n')
        for i=1:numel(infos.contrast_name)
            fprintf('%4d %s\n',i,infos.contrast_name{i});
        end

    case {'C' 'conw' 'consw'}
        if s
            infos.contrast_name{1} = SPM.xCon(s).name;
            infos.beta_numbers{1} = find(SPM.xCon(s).c);
            infos.beta_names{1} = cellstr(strvcat(SPM.xX.name{infos.beta_numbers{1}})); %#ok
            infos.beta_weights{1} = SPM.xCon(s).c(infos.beta_numbers{1});
        else
            infos.contrast_name = cellstr(strvcat(SPM.xCon.name)); %#ok
            infos.contrast_num = 1:numel(infos.contrast_name);
            for i = 1:numel(infos.contrast_name);
                infos.beta_numbers{i} = find(SPM.xCon(i).c);
                infos.beta_names{i} = cellstr(strvcat(SPM.xX.name{infos.beta_numbers{i}})); %#ok
                infos.beta_weights{i} = SPM.xCon(i).c(infos.beta_numbers{i});
            end
        end

        for s = 1:numel(infos.contrast_name)
            fprintf('%3d: %s:\n',s,infos.contrast_name{s})
            for i=1:numel(infos.beta_names{s})
                fprintf('\t%10.5f\t%s\n',infos.beta_weights{s}(i),infos.beta_names{s}{i});
            end
        end

    case {'b' 'betas'}
        if s
            infos.beta_number = find(~cellfun('isempty',regexp(SPM.xX.name,sprintf('^Sn\\(%d\\)',s))))';
            infos.beta_name = cellstr(strvcat(SPM.xX.name{infos.beta_number})); %#ok
        else
            infos.beta_name = cellstr(strvcat(SPM.xX.name)); %#ok
            infos.beta_number = [1:numel(infos.beta_name)]'; %#ok
        end

        for i = 1:numel(infos.beta_name)
            fprintf('%4d %s\n',infos.beta_number(i),infos.beta_name{i});
        end

    case {'B' 'taskbetas'}
        n = cellfun('isempty',regexp(SPM.xX.name,' R[0-9]*$'))';
        if s
            n = n & ~cellfun('isempty',regexp(SPM.xX.name,sprintf('^Sn\\(%d\\)',s)))';
        end
        n = n & cellfun('isempty',regexp(SPM.xX.name,'constant$'))';
        infos.beta_number = find(n);
        infos.beta_name = cellstr(strvcat(SPM.xX.name{infos.beta_number})); %#ok

        for i = 1:numel(infos.beta_number)
            fprintf('%4d %s\n',infos.beta_number(i),infos.beta_name{i});
        end

    case 'nsess'
        infos.nsess = numel(SPM.Sess);
        fprintf('%d\n',infos.nsess);

    case {'i' 'input' 'sess'}
        infos.nsess = numel(SPM.nscan);
        infos.nscan = SPM.nscan;

        n=1; for i=1:infos.nsess, first(i)=n; n=n+infos.nscan(i); end %#ok

        fprintf('%4s %5s  %s\n','sess','nscan','first_volume_name')
        for i = 1:numel(first)
            infos.first_volume_name{i} = SPM.xY.VY(first(i)).fname;
            fprintf('%4d %5d  %s\n',i,infos.nscan(i),infos.first_volume_name{i})
        end

    case {'v' 'image'}
        if s
            canlab_glm_getinfo(modeldir,'betas',s);
            infos.X = SPM.xX.X(SPM.Sess(s).row,SPM.Sess(s).col);
        else
            canlab_glm_getinfo(modeldir,'betas');
            infos.X = SPM.xX.X;
        end

        figure; image(infos.X * 40);
        colormap('Bone');

    case {'V' 'taskimage'}
%         n = cellfun('isempty',regexp(SPM.xX.name,' R[0-9]*$'));
        if s
            n = canlab_glm_getinfo(modeldir,'taskbetas',s);
            %n = n .* ~cellfun('isempty',regexp(SPM.xX.name,['^Sn\(' num2str(s) '\)']));
            infos.X = SPM.xX.X(SPM.Sess(s).row,n.beta_number);
        else
            n = canlab_glm_getinfo(modeldir,'taskbetas');
            infos.X = SPM.xX.X(:,n.beta_number);
        end
        figure; image(infos.X * 400)
        colormap('Bone');

    case {'vsc' 'imagesc'}
        spmlowerinfo(modeldir,'betas')
        if s
            canlab_glm_getinfo(modeldir,'betas',s)
            figure; imagesc(SPM.xX.X(SPM.Sess(s).row,SPM.Sess(s).col))
            colormap('Bone');
        else
            canlab_glm_getinfo(modeldir,'betas')
            figure; imagesc(SPM.xX.X)
            colormap('Bone');
        end

    otherwise
        error('Unrecognized option for spm lower level analyses: %s',opt)
end

end

function [infos] = get_group_level_info(modeldir,opt,s)

infos = [];

% load SETUP structs
f = filenames(fullfile(modeldir,'robust[0-9][0-9][0-9][0-9]'));
for i=1:numel(f),
    load(fullfile(f{i},'SETUP'));
    SETUPS(i) = SETUP; %#ok
    clear SETUP
end

switch opt
    case {'i' 'input'}
        % load lower levels
        evalc('tmp = canlab_glm_getinfo(fileparts(SETUPS(i).V(1).fname),''cons'');');
        connames = tmp.contrast_name;
        if ~s, s = 1:numel(SETUPS); end
        fprintf('%-12s %-12s %-s\n','dir','input_con','con_name')
        for i = s
            [ignore finf] = fileparts(SETUPS(i).V(1).fname); %#ok

            fprintf('%-12s',SETUPS(i).dir)
            fprintf(' %-12s',finf)
            if regexp(finf,'^con_')
                fprintf(' "%-s"',connames{str2num(regexprep(finf,'^con_',''))}) %#ok
            end
            fprintf('\n')
        end
    case {'I' 'allinput'}
        if ~s, s = 1:numel(SETUPS); end
        for i = s
            fprintf('%s:\n',SETUPS(i).dir)
            disp(SETUPS(i).files)
        end
    case {'m' 'model'}
        fprintf('weight input\n')
        for i = 1:numel(SETUPS(1).X)
            fprintf('%-7.3f %s\n',SETUPS(1).X(i),regexprep(SETUPS(1).V(i).fname,'/[^/]*$',''))
        end
    case {'M' 'allmodels'}
        if ~s, s = 1:numel(SETUPS); end
        for m = s
            fprintf('%s\n',SETUPS(m).dir)
            fprintf('weight input\n')
            for i = 1:numel(SETUPS(m).X)
                fprintf('%-7.3f %s\n',SETUPS(m).X(i),SETUPS(m).V(i).fname)
            end
        end

    otherwise
        error('Unrecognized option for robfit analyses: %s',opt)
end

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [gf,logE,X]=FE_GT(sam,fs,NbCh,FrameLen,tGtMat)

Efloor=exp(-50);
Gfloor=exp(-50);
PreEmp=0.97;

if ~exist('fs', 'var')
   fs=16000; % sampling frequency
end
if ~exist('FrameLen', 'var')
   FrameLen=0.025*fs;
end
FrameShift=0.01*fs;
Nfft=2^ceil(log(FrameLen)/log(2));

if ~exist('NbCh', 'var')
    NbCh=64; % number of GT filter banks
end
if ~exist('tGtMat', 'var')
    tGtMat=gammatone_matrix(Nfft,fs,NbCh)';
end
% truncate as in ReadWave
NbFr=floor( (length(sam)-FrameLen+FrameShift)/FrameShift);
sam=sam(1:NbFr*FrameShift+FrameLen-FrameShift);

% DC removal
sam=filter([1 -1],[1 -0.999],[0 reshape(sam,1,length(sam))]);

% framing
ind1=1:FrameShift:length(sam)-1-FrameLen+FrameShift;
ind2=(1:FrameLen)';
x=sam(ind1(ones(FrameLen,1),:)+ind2(:,ones(1,NbFr)));

% logE
logE=log(max(sum(x.^2,1),Efloor));

% preemphasis & windowing
T=length(sam);
sam=[0 sam(2:T)-PreEmp*sam(1:T-1)];
win=hamming(FrameLen);
% Fourier Transform
X=fft(win(:,ones(1,NbFr)).*sam(ind1(ones(FrameLen,1),:)+ind2(:,ones(1,NbFr))),Nfft);
X(1,:)=0;
X=abs(X(1:Nfft/2,:));
gf=max(tGtMat'*X,Gfloor).^(1/3);

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
classdef TRAFO < handle
    % TRAFO Transformation class to overload mathematical operations
    %
    % (c) Thomas Kuestner
    % ---------------------------------------------------------------------

    properties
        data;
        meta;
    end

    methods
        function obj = TRAFO(data,meta,paraTrafo,initialize)
            if(nargin < 3)
                obj.data = data;
                obj.meta = meta;
            else
                if(nargin < 4)
                    initialize = 'default';
                end
                [obj.data, obj.meta] = initTRAFO(data,paraTrafo,meta,initialize);
%                 if(isempty(obj.meta))
%                     obj.meta = meta;
%                 end
            end
        end

        function data = getData(obj)
            data = obj.data;
        end

        function meta = getMeta(obj)
            meta = obj.meta;
        end

    end

end

author: thomaskuestner


text:
function locpts = glob2loc_3d ( cospitch, cosroll, cosyaw, sinpitch, ...
  sinroll, sinyaw, globas, glopts )

%*****************************************************************************80
%
%% GLOB2LOC_3D converts from a global to a local coordinate system in 3D.
%
%  Discussion:
%
%    A global coordinate system is given.
%
%    A local coordinate system has been translated to the point with
%    global coordinates GLOBAS, and rotated through a yaw, a pitch, and
%    a roll.
%
%    A point has global coordinates GLOPTS, and it is desired to know
%    the point's local coordinates LOCPTS.
%
%    The transformation may be written as
%
%      LOC = M_ROLL * M_PITCH * M_YAW * ( GLOB - GLOBAS )
%
%    where
%
%               (       1            0            0      )
%    M_ROLL =   (       0        cos(Roll)    sin(Roll)  )
%               (       0      - sin(Roll)    cos(Roll)  )
%
%               (   cos(Pitch)       0      - sin(Pitch) )
%    M_PITCH =  (       0            1            0      )
%               (   sin(Pitch)       0        cos(Pitch) )
%
%               (   cos(Yaw)     sin(Yaw)         0      )
%    M_YAW    = ( - sin(Yaw)     cos(Yaw)         0      )
%               (       0            0            1      )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real COSPITCH, COSROLL, COSYAW, the cosines of
%    the pitch, roll and yaw angles.
%
%    Input, real SINPITCH, SINROLL, SINYAW, the sines of the pitch,
%    roll and yaw angles.
%
%    Input, real GLOBAS(3), the X, Y, and Z coordinates of the
%    global base vector.
%
%    Input, real GLOPTS(3), the global ( X, Y and Z ) coordinates
%    of the point whose coordinates are to be transformed.
%
%    Output, real LOCPTS(3), the local coordinates of the point
%    whose global coordinates were given in GLOPTS.
%
  locpts(1) = ( cosyaw * cospitch ) * ( glopts(1) - globas(1) ) ...
            + ( sinyaw * cospitch ) * ( glopts(2) - globas(2) ) ...
            -   sinpitch            * ( glopts(3) - globas(3) );

  locpts(2) = ( cosyaw * sinpitch * sinroll - sinyaw * cosroll ) ...
    * ( glopts(1) - globas(1) ) ...
    + ( sinyaw * sinpitch * sinroll + cosyaw * cosroll ) ...
    * ( glopts(2) - globas(2) ) ...
    +   cospitch * sinroll * ( glopts(3) - globas(3) );

  locpts(3) = ( cosyaw * sinpitch * cosroll + sinyaw * sinroll ) ...
    * ( glopts(1) - globas(1) ) ...
    + ( sinyaw * sinpitch * cosroll - cosyaw * sinroll  ) ...
    * ( glopts(2) - globas(2) ) ...
    + ( cospitch * cosroll ) * ( glopts(3) - globas(3) );

  return
end


text:
function a = diag(a,k)
%DIAG         Implements  diag(a,k)  for Taylor
%
%   c = diag(a,k)
%
% functionality as Matlab function diag for matrices
%

% written  05/21/09     S.M. Rump
% modified 08/26/12     S.M. Rump  global variables removed
%

  if nargin==1
    k = 0;
  end

  if any(a.size==1)
    K1 = getappdata(0,'INTLAB_TAYLOR_ORDER') + 1;
    index = ( diag(1:prod(a.size),k)~= 0 );
    a.size = size(index);
    at = a.t;
    a.t = zeros(K1,prod(size(index)));
    if isa(at,'intval')
      a.t = intval(a.t);
    end
    a.t(:,index) = at;
  else
    index = 1:prod(a.size);
    index = diag(reshape(index,a.size),k);
    a.size = [length(index) 1];
    a.t = a.t(:,index);
  end

repo: OpenMAS


text:
function colors = custom_colors(startrgb, endrgb, k)
% Define cell array of custom colors in a specific range of the spectrum,
% to be used in visualization.
%
% :Usage:
% ::
%     colors = custom_colors(startrgb, endrgb, k [length])
% ..
%     Author and copyright information:
%
%     Copyright (C) 2016 Tor Wager
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
% ..
%
% :Inputs:
%
%   **startrgb**
%        3-color RGB triplet, values 0-1, for starting color in range
%
%   **endrgb**
%        3-color RGB triplet, values 0-1, for starting color in range
%
%   **k**
%        how many colors you want, from 1 - 64. Colors will be sampled
%        evenly along the range you specify
%
% :Outputs:
%
%   **colors_cell**
%        Cell array of k colors
%
% :Examples:
% ::
%
%    % give examples of code here
%    param1 = abc();
%    param2 = xyz();
%    [out1,out2] = func_call(param1, param2)
%
% :See also:
%   colormap_tor, bucknerlab_colors, scn_standard_colors
%

% define 64-color map in range you specify
% ---------------------------------------------
cm = colormap_tor(startrgb, endrgb); % custom colormap

% turn it into a cell array and sample only the number desired
% ---------------------------------------------

whcolors = round(linspace(1, size(cm, 1), k));

colors = mat2cell(cm(whcolors, :), ones(k, 1));

end % function


text:
function test_bug3343

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

%% this used to return ctf
elec = [];
elec.elecpos = randn(3,2);
elec.chanpos = elec.elecpos;
elec.label = {'BR1';'BR2'};

assert(~strcmp(ft_senstype(elec), 'ctf'))
assert( strcmp(ft_senstype(elec), 'eeg')) % because it has elecpos

%% this used to return ctf
elec = [];
elec.elecpos = randn(3,2);
elec.chanpos = elec.elecpos;
elec.label = {'BR1';'BR2'};
elec.chantype = {'eeg';'eeg'};

assert(~strcmp(ft_senstype(elec), 'ctf'))
assert( strcmp(ft_senstype(elec), 'eeg')) % because of the chantype

%% use boolean in detection
elec.elecpos = randn(3,2);
elec.chanpos = elec.elecpos;
elec.label = {'BR1';'BR2'};
elec.chantype = {'seeg';'ecog'};

assert( ft_senstype(elec, 'eeg'))
assert( ft_senstype(elec, 'ieeg'))
assert(~ft_senstype(elec, 'seeg')) % not all channels are like this
assert(~ft_senstype(elec, 'ecog')) % not all channels are like this

%% use boolean in detection
elec.elecpos = randn(3,2);
elec.chanpos = elec.elecpos;
elec.label = {'BR1';'BR2'};
elec.chantype = {'seeg';'seeg'};

assert( ft_senstype(elec, 'eeg'))
assert( ft_senstype(elec, 'ieeg'))
assert( ft_senstype(elec, 'seeg'))
assert(~ft_senstype(elec, 'ecog'))

%% use boolean in detection
elec.elecpos = randn(3,2);
elec.chanpos = elec.elecpos;
elec.label = {'BR1';'BR2'};
elec.chantype = {'ecog';'ecog'};

assert( ft_senstype(elec, 'eeg'))
assert( ft_senstype(elec, 'ieeg'))
assert( ft_senstype(elec, 'ecog'))
assert(~ft_senstype(elec, 'seeg'))

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function point_num = sparse_grid_ofn_size ( dim_num, level_max )

%*****************************************************************************80
%
%% SPARSE_GRID_OFN_SIZE sizes a sparse grid using Open Fully Nested rules.
%
%  Discussion:
%
%    The grid is defined as the sum of the product rules whose LEVEL
%    satisfies:
%
%      0 <= LEVEL <= LEVEL_MAX.
%
%    This calculation is much faster than a previous method.  It simply
%    computes the number of new points that are added at each level in the
%    1D rule, and then counts the new points at a given DIM_NUM dimensional
%    level vector as the product of the new points added in each dimension.
%
%    This approach will work for nested families, and may be extensible
%    to other families, and to mixed rules.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 December 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    A Sparse Grid Stochastic Collocation Method for Partial Differential
%    Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2309-2345.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer LEVEL_MAX, the maximum value of LEVEL.
%
%    Output, integer POINT_NUM, the total number of unique
%    points in the grids.
%

%
%  Special case.
%
  if ( level_max < 0 )
    point_num = 0;
    return
  end

  if ( level_max == 0 )
    point_num = 1;
    return
  end
%
%  Construct the vector that counts the new points in the 1D rule.
%
  new_1d = zeros ( level_max+1, 1 );

  new_1d(0+1) = 1;

  for l = 1 : level_max
    new_1d(l+1) = 2 * new_1d(l);
  end
%
%  Count the number of points by counting the number of new points
%  associated with each level vector.
%
  level_1d = zeros ( dim_num, 1 );

  point_num = 0;

  for level = 0 : level_max

    more = 0;
    h = 0;
    t = 0;

    while ( 1 )

      [ level_1d, more, h, t ] = comp_next ( level, dim_num, level_1d, more, h, t );

      point_num = point_num + prod ( new_1d(level_1d(1:dim_num)+1) );

      if ( ~more )
        break
      end

    end

  end

  return
end


text:
clear
clc
G1=[1.15,1.46,1.67,1.82,1.94,2.03,2.11,2.18,2.23,2.29,2.33,2.37,2.41,2.44,2.47,2.50,2.53,2.56]
G2=[1.15,1.49,1.75,1.94,2.10,2.22,2.32,2.41,2.48,2.55,2.61,2.66,2.70,2.74,2.78,2.82,2.85,2.85]
grubbstab=input('Please input measuring datas:')
g=input('please input confidence level ')
if g==0.99
    G=G2
else
    G=G1
end
meangrubbstab =mean(grubbstab)
stdgrubbstab=std(grubbstab)
len=length(grubbstab)
for i=1:len
  A(i)=abs(grubbstab(i)-meangrubbstab)
end
Max1=max(A)
location1=find(A==Max1)
if    abs(grubbstab(location1)-meangrubbstab)< G(len-2)*stdgrubbstab
      disp('no wide error!')
else
       grubbstab(location1)=[]
       A(location1)=[]
       meangrubbstab =mean(grubbstab)
       stdgrubbstab=std(grubbstab)
       Max2=max(A)
       location2=find(A==Max2)
       len=length(grubbstab)
       if    abs(grubbstab(location2)-meangrubbstab)< G(len-2)*stdgrubbstab
             disp('no more wide error!')
       else
       grubbstab(location2)=[]
       end
end

author: sfvsfv

repo: Mathematical-modeling


text:
function jed = yjf_to_jed_english ( y, j, f )

%*****************************************************************************80
%
%% YJF_TO_JED_ENGLISH converts an English YJF date to a JED.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, J, real F, the YJF date.
%
%    Output, real JED, the Julian Ephemeris Date.
%

%
%  Copy the input.
%
  y1 = y;
  j1 = j;
  f1 = f;
%
%  Check the input.
%
  [ y1, j1, ierror ] = yj_check_english ( y1, j1 );

  if ( ierror ~= 0 )
    jed = -1.0;
    return
  end
%
%  Convert the input.
%
  [ y2, m2, d2, f2 ] = yjf_to_ymdf_english ( y1, j1, f1 );

  jed = ymdf_to_jed_english ( y2, m2, d2, f2 );

  return
end


text:
function [posterior] = tapas_vlinear_estimate_interface(data, model, ...
    inference)
%%
%
% aponteeduardo@gmail.com
% copyright (C) 2016
%

% Verify input
tapas_validate_data(data);
tapas_validate_model(model);
tapas_validate_inference(inference);

[data] = tapas_linear_prepare_data(data, model, inference);
[model] = tapas_vlinear_prepare_model(data, model, inference);
[inference] = tapas_vlinear_prepare_inference(data, model, inference);

[posterior] = inference.estimate_method(data, model, inference);

end

author: translationalneuromodeling

repo: tapas


text:
% ARG model: the hsvargplvm model
% ARG X: the test latent points of the TOP layer (parent)
% ARG varX: the variance associated with X. Can be left empty ([])

% TODO: Replace vargplvmPosteriorMeanVarHier2 with
% vargplvmPosteriorMeanVarHier which is an external function (check what
% changes were made)

% This function is very similar to hsvargplvmPosteriorMeanVar, but without
% the extra arguments.

function [mu, varsigma] = hsvargplvmPosteriorMeanVarSimple(model, X, varX)

if nargin < 3
    varX = [];
end

H = model.H;

Xall{H} = X;
varXall{H} = varX;

for h=H-1:-1:1
    if ~isempty(varX)
        [Xall{h} varXall{h}] = vargplvmPosteriorMeanVarHier2(model.layer{h+1}.comp{1}, Xall{h+1}, varXall{h+1});
    else
        [Xall{h} varXall{h}] = vargplvmPosteriorMeanVarHier2(model.layer{h+1}.comp{1}, Xall{h+1});
    end
end

if ~isempty(varX)
    [mu, varsigma] = vargplvmPosteriorMeanVarHier2(model.layer{1}.comp{1}, Xall{1}, varXall{1});
else
    [mu, varsigma] = vargplvmPosteriorMeanVarHier2(model.layer{1}.comp{1}, Xall{1});
end

% h=1;sc = vargplvmRetainedScales(model.layer{h}.comp{1}); close all; plotDistr(varXall{h}(:,sc));

% The same as vargplvmPosteriorMeanVar, but a small change in the
% calculations involving model.m, because in the indermediate layers
% model.m is replaced by the expectation varmu*varmu'+Sn wrt the X of the
% bottom layer.
function [mu, varsigma] = vargplvmPosteriorMeanVarHier2(model, X, varX)

if nargin < 3
  vardistX.covars = repmat(0.0, size(X, 1), size(X, 2));%zeros(size(X, 1), size(X, 2));
else
  vardistX.covars = varX;
end
vardistX.latentDimension = size(X, 2);
vardistX.numData = size(X, 1);
%model.vardist.covars = 0*model.vardist.covars;
vardistX.means = X;
%model = vargplvmUpdateStats(model, model.X_u);

Ainv = model.P1' * model.P1; % size: NxN

if ~isfield(model,'alpha')
    if isfield(model, 'mOrig')
        model.alpha = Ainv*model.Psi1'*model.mOrig; % size: 1xD
    else
        model.alpha = Ainv*model.Psi1'*model.m; % size: 1xD
    end
end
Psi1_star = kernVardistPsi1Compute(model.kern, vardistX, model.X_u);

% mean prediction
mu = Psi1_star*model.alpha; % size: 1xD

if nargout > 1
   %
   % precomputations
   vard = vardistCreate(zeros(1,model.q), model.q, 'gaussian');
   Kinvk = (model.invK_uu - (1/model.beta)*Ainv);
   %
   for i=1:size(vardistX.means,1)
      %
      vard.means = vardistX.means(i,:);
      vard.covars = vardistX.covars(i,:);
      % compute psi0 term
      Psi0_star = kernVardistPsi0Compute(model.kern, vard);
      % compute psi2 term
      Psi2_star = kernVardistPsi2Compute(model.kern, vard, model.X_u);

      vars = Psi0_star - sum(sum(Kinvk.*Psi2_star));

      for j=1:model.d
         %[model.alpha(:,j)'*(Psi2_star*model.alpha(:,j)), mu(i,j)^2]
         varsigma(i,j) = model.alpha(:,j)'*(Psi2_star*model.alpha(:,j)) - mu(i,j)^2;
      end
      varsigma(i,:) = varsigma(i,:) + vars;
      %
   end
   %
   if isfield(model, 'beta')
      varsigma = varsigma + (1/model.beta);
   end
   %
end

% Rescale the mean
mu = mu.*repmat(model.scale, size(vardistX.means,1), 1);

% Add the bias back in
mu = mu + repmat(model.bias, size(vardistX.means,1), 1);

% rescale the variances
if nargout > 1
    varsigma = varsigma.*repmat(model.scale.*model.scale, size(vardistX.means,1), 1);
end

author: SheffieldML

repo: deepGP

save_path: github-repos/MATLAB/SheffieldML-deepGP


text:
function err = getH1errorQ1(node,elem,Du,uh,K,quadOrder)
%% GETH1ERRORQ1 H1 norm of the approximation error on quad mesh.
%
%  err = getH1errorQ1(node,elem,@Du,uh) computes the H1 norm of the
%  error between the exact solution Du and finite element approximation
%  uh on a quad mesh described by node and elem.
%
%  The input parameter Du is a function handle uh is a column array with
%  length N representing a piecewise linear function on the mesh.
%
% Author: Huayi Wei < huayiwei1984@gmail.com>
% Modified by Long Chen. The original computation is not quite right. The
% diffusion coefficient cannot be computed separately.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

[N,  dim] = size(node);
[NT, nv] = size(elem);
if ~exist('quadOrder','var')
    quadOrder = 3;
end

%% compute gradient of basis.
[pts, weight] = quadquadpts(quadOrder);
nQuad = size(pts,1);
err = zeros(size(elem,1),1);
for p = 1:nQuad
    [phi, Dphi, J] = quadbasis(node,elem, pts(p,:));
    pxy = zeros(NT, dim);
    for i = 1:dim
        xi = node(:,i);
        pxy(:,i) = xi(elem)*phi;
    end
    Duh = zeros(NT,dim);
    for i = 1:dim
        Dphi_i(:,1:nv) = Dphi(:,i,:);
        Duh(:,i) = sum(uh(elem).*Dphi_i,2);
    end
    if exist('K','var') && ~isempty(K) && ~isnumeric(K) % K is a function
        err = err + weight(p)*K(pxy).*sum((Du(pxy)-Duh).^2,2);
    else
        err = err + weight(p)*sum((Du(pxy)-Duh).^2,2);
    end
end
if exist('K','var') && ~isempty(K) && isnumeric(K) && size(K,1) == NT
    err = K.*err;    % K is piecewise constant
end
err = J.*err;
err(isnan(err)) = 0; % singular values are excluded
err = sqrt(sum(err));
end
%% TODO: add more help information.

repo: Aero-matlab


text:
%% DEMO_febio_0022_multigen_interface_band
% Below is a demonstration for:
%
% * Building geometry for limb-like segment with an elastic band wrapped around it
% * Defining the boundary conditions
% * Coding the febio structure
% * Running the model
% * Importing and visualizing the displacement results

%% Keywords
%
% * febio_spec version 3.0
% * febio, FEBio
% * indentation
% * contact, sliding, sticky, friction
% * rigid body constraints
% * tetrahedral elements, tet4
% * static, solid
% * multi-generational
% * hyperelastic, Ogden
% * displacement logfile
% * stress logfile

%%

clear; close all; clc;

%% Plot settings
fontSize=15;
faceAlpha1=0.8;
faceAlpha2=0.3;
markerSize=40;
lineWidth=3;
plotColors=gjet(9);
lineWidth1=2;
markerSize1=25;

%% Control parameters

% Path names
defaultFolder = fileparts(fileparts(mfilename('fullpath')));
savePath=fullfile(defaultFolder,'data','temp');

% Defining file names
febioFebFileNamePart='tempModel';
febioFebFileName=fullfile(savePath,[febioFebFileNamePart,'.feb']); %FEB file name
febioLogFileName=[febioFebFileNamePart,'.txt']; %FEBio log file name
febioLogFileName_disp=[febioFebFileNamePart,'_disp_out.txt']; %Log file name for exporting displacement
febioLogFileName_strainEnergy=[febioFebFileNamePart,'_energy_out.txt']; %Log file name for exporting strain energy density

%Material parameter set
c1=1e-3; %Shear-modulus-like parameter
m1=2; %Material parameter setting degree of non-linearity
k_factor=1e2; %Bulk modulus factor
k=c1*k_factor; %Bulk modulus

%Multi-generational properties of band
c1_g=[c1/1000 c1*100];
k_g=c1_g*k_factor;

% FEA control settings
numTimeSteps=6; %Number of time steps desired
max_refs=25; %Max reforms
max_ups=0; %Set to zero to use full-Newton iterations
opt_iter=10; %Optimum number of iterations
max_retries=5; %Maximum number of retires
dtmin=(1/numTimeSteps)/100; %Minimum time step size
dtmax=1/numTimeSteps; %Maximum time step size

appliedPressure=3e-3; %pressure value

% Geometry parameters
tissueRadius=35;
tissueHeight=150;
boneRadius=10;
wrapHeight=24;
wrapThickness=5;
pointSpacing=6; % Aproximate node spacing

%% Build tissue skin surface top

%Sketching profile
ns=150;
t=linspace(0,2*pi,ns);
t=t(1:end-1);

x=tissueRadius*cos(t);
y=tissueRadius*sin(t);
z=zeros(size(x));
Vc=[x(:) y(:) z(:)];
np=ceil(max(pathLength(Vc))./pointSpacing);
[Vc]=evenlySampleCurve(Vc,np,'pchip',1);

% Extruding model
h=(tissueHeight/2)-wrapHeight/2;
cPar.numSteps=round(h/pointSpacing);
cPar.numSteps=cPar.numSteps+iseven(cPar.numSteps);
cPar.depth=h;
cPar.patchType='tri';
cPar.dir=-1;
cPar.closeLoopOpt=1;
[Fg1,Vg1]=polyExtrude(Vc,cPar);
Vg1(:,3)=Vg1(:,3)+tissueHeight/2;

[T,R] = cart2pol(Vg1(:,1),Vg1(:,2));
R=R-(tissueRadius/2*((Vg1(:,3)-(wrapHeight/2))/(tissueHeight/2)).^2);
[Vg1(:,1),Vg1(:,2)] = pol2cart(T,R);

Vg1b=Vg1(cPar.numSteps:cPar.numSteps:end,:);
Vg1t=Vg1(1:cPar.numSteps:end,:);

%% Build tissue skin surface bottom

Fg2=Fg1;
Fg2=fliplr(Fg2);
Vg2=Vg1;
Vg2(:,3)=-Vg2(:,3);
Vg2b=Vg2(cPar.numSteps:cPar.numSteps:end,:);
Vg2t=Vg2(1:cPar.numSteps:end,:);

%% Build tissue skin surface middle

% Extruding model
h=wrapHeight;
cPar.numSteps=round(h/pointSpacing);
cPar.numSteps=cPar.numSteps+iseven(cPar.numSteps);
cPar.depth=h;
cPar.patchType='tri';
cPar.closeLoopOpt=1;

Vc_start=Vg1b;
Vc_end=Vg2b;
[Fg3,Vg3]=polyLoftLinear(Vc_start,Vc_end,cPar);

% Vg3(:,3)=Vg3(:,3)+tissueHeight/2;
Vg3b=Vg3(cPar.numSteps:cPar.numSteps:end,:);
Vg3t=Vg3(1:cPar.numSteps:end,:);

%% Build wrap outer surface

ns=150;
t=linspace(0,2*pi,ns);
t=t(1:end-1);

x=(tissueRadius+wrapThickness)*cos(t);
y=(tissueRadius+wrapThickness)*sin(t);
z=zeros(size(x));
Vc=[x(:) y(:) z(:)];
np=ceil(max(pathLength(Vc))./pointSpacing);
[Vwt]=evenlySampleCurve(Vc,np,'pchip',1);
Vwt(:,3)=mean(Vg3t(:,3));
Vwb=Vwt;
Vwb(:,3)=Vwt(:,3)-wrapHeight;

h=wrapHeight;
cPar.numSteps=round(h/pointSpacing);
cPar.numSteps=cPar.numSteps+iseven(cPar.numSteps);
cPar.depth=h;
cPar.patchType='tri';
cPar.closeLoopOpt=1;

Vc_start=Vwt;
Vc_end=Vwb;
[Fw1,Vw1]=polyLoftLinear(Vc_start,Vc_end,cPar);

%% Build bone surface

x=boneRadius*cos(t);
y=boneRadius*sin(t);
z=zeros(size(x));
Vc=[x(:) y(:) z(:)];
np=ceil(max(pathLength(Vc))./pointSpacing);
[Vc]=evenlySampleCurve(Vc,np,'pchip',1);

% Extruding model
cPar.numSteps=round(tissueHeight/pointSpacing);
cPar.numSteps=cPar.numSteps+iseven(cPar.numSteps);
cPar.depth=tissueHeight;
cPar.patchType='tri';
cPar.dir=-1;
cPar.closeLoopOpt=1;
[Fb,Vb]=polyExtrude(Vc,cPar);
Fb=fliplr(Fb);
Vb(:,3)=Vb(:,3)+tissueHeight/2;

Vbb=Vb(cPar.numSteps:cPar.numSteps:end,:);
Vbt=Vb(1:cPar.numSteps:end,:);

%% Capping tissue top

regionCell={Vg1t(:,[1 2]),Vbt(:,[1 2])};
[Ft,Vt]=regionTriMesh2D(regionCell,pointSpacing,0,0);
Vt(:,3)=mean(Vg1t(:,3));
Ft=fliplr(Ft);

%% Capping tissue bottom

regionCell={Vg2t(:,[1 2]),Vbb(:,[1 2])};
[Fgb,Vgb]=regionTriMesh2D(regionCell,pointSpacing,0,0);
Vgb(:,3)=mean(Vg2t(:,3));

%% Capping wrap top

regionCell={Vwt(:,[1 2]),Vg1b(:,[1 2])};
[Fwtt,Vwtt]=regionTriMesh2D(regionCell,pointSpacing,0,0);
Vwtt(:,3)=mean(Vwt(:,3));
Fwtt=fliplr(Fwtt);

%% Capping wrap bottom

regionCell={Vwb(:,[1 2]),Vg2b(:,[1 2])};
[Fwbb,Vwbb]=regionTriMesh2D(regionCell,pointSpacing,0,0);
Vwbb(:,3)=mean(Vwb(:,3));

%% Visualizing surface geometry

cFigure;
xlabel('X','FontSize',fontSize); ylabel('Y','FontSize',fontSize); zlabel('Z','FontSize',fontSize)
hold on;

gpatch(Fg1,Vg1,plotColors(1,:),'k');
patchNormPlot(Fg1,Vg1);
plotV(Vg1t,'r.-','lineWidth',lineWidth1,'MarkerSize',markerSize1);
plotV(Vg1b,'y.-','lineWidth',lineWidth1,'MarkerSize',markerSize1);

gpatch(Fg2,Vg2,plotColors(2,:),'k');
patchNormPlot(Fg2,Vg2);
plotV(Vg2t,'r.-','lineWidth',lineWidth1,'MarkerSize',markerSize1);
plotV(Vg2b,'y.-','lineWidth',lineWidth1,'MarkerSize',markerSize1);

gpatch(Fg3,Vg3,plotColors(3,:),'k');
patchNormPlot(Fg3,Vg3);
plotV(Vg3t,'r.-','lineWidth',lineWidth1,'MarkerSize',markerSize1);
plotV(Vg3b,'y.-','lineWidth',lineWidth1,'MarkerSize',markerSize1);

gpatch(Fw1,Vw1,plotColors(4,:),'k');
patchNormPlot(Fw1,Vw1);
plotV(Vwt,'g.-','lineWidth',lineWidth1,'MarkerSize',markerSize1);
plotV(Vwb,'g.-','lineWidth',lineWidth1,'MarkerSize',markerSize1);

gpatch(Fb,Vb,plotColors(5,:),'k');
patchNormPlot(Fb,Vb);
plotV(Vbt,'r.-','lineWidth',lineWidth1,'MarkerSize',markerSize1);
plotV(Vbb,'y.-','lineWidth',lineWidth1,'MarkerSize',markerSize1);

gpatch(Ft,Vt,plotColors(6,:),'k');
patchNormPlot(Ft,Vt);

gpatch(Fgb,Vgb,plotColors(7,:),'k');
patchNormPlot(Fgb,Vgb);

gpatch(Fwtt,Vwtt,plotColors(8,:),'k');
patchNormPlot(Fwtt,Vwtt);

gpatch(Fwbb,Vwbb,plotColors(9,:),'k');
patchNormPlot(Fwbb,Vwbb);

axisGeom;
colormap(plotColors); colorbar;
drawnow;

%% Joining and merging geometry sets

%Joining sets

%Creating color information
Cg1=1*ones(size(Fg1,1),1); %Tissue top cylinder
Cg2=2*ones(size(Fg2,1),1); %Tissue bottom cylinder
Cg3=3*ones(size(Fg3,1),1); %Tissue middle cylinder
Cw1=4*ones(size(Fw1,1),1); %Wrap outer cylinder
Cb=5*ones(size(Fb,1),1); %Bone cylinder
Ct=6*ones(size(Ft,1),1); %Tissue top
Cgb=7*ones(size(Fgb,1),1); %Tissue bottom
Cwtt=8*ones(size(Fwtt,1),1); %Wrap top
Cwbb=9*ones(size(Fwbb,1),1); %Wrap bottom

[F,V,C]=joinElementSets({Fg1,Fg2,Fg3,Fw1,Fb,Ft,Fgb,Fwtt,Fwbb},{Vg1,Vg2,Vg3,Vw1,Vb,Vt,Vgb,Vwtt,Vwbb},{Cg1,Cg2,Cg3,Cw1,Cb,Ct,Cgb,Cwtt,Cwbb}); %joining sets together

%merging sets
[F,V]=mergeVertices(F,V);

%%

cFigure;
xlabel('X','FontSize',fontSize); ylabel('Y','FontSize',fontSize); zlabel('Z','FontSize',fontSize)
hold on;

gpatch(F,V,C,'none',0.5);
% patchNormPlot(F,V);

axisGeom;
colormap(plotColors); icolorbar;
drawnow;

%% Find solid mesh region interior points

logicRegion=ismember(C,[3 4 8 9]);
[V_in_1]=getInnerPoint(F(logicRegion,:),V);

logicRegion=ismember(C,[1 2 3 5 6 7 ]);
[V_in_2]=getInnerPoint(F(logicRegion,:),V);

V_regions=[V_in_1;V_in_2];

%%

cFigure;
xlabel('X','FontSize',fontSize); ylabel('Y','FontSize',fontSize); zlabel('Z','FontSize',fontSize)
hold on;

gpatch(F,V,C,'none',0.2);
plotV(V_in_1,'r.','MarkerSize',25);
plotV(V_in_2,'b.','MarkerSize',25);

axisGeom;
colormap(plotColors); colorbar;
drawnow;

%% Mesh solid using tetgen

% Create tetgen meshing input structure
modelName=fullfile(savePath,'tempModel');

% Regional mesh volume parameter
[regionA]=tetVolMeanEst(F,V); %Volume for a regular tet based on edge lengths
volumeFactors=(regionA.*ones(size(V_regions,1),1));

inputStruct.stringOpt='-pq1.2AaY';
inputStruct.Faces=F;
inputStruct.Nodes=V;
inputStruct.holePoints=[];
inputStruct.faceBoundaryMarker=C; %Face boundary markers
inputStruct.regionPoints=V_regions; %region points
inputStruct.regionA=volumeFactors; %Desired volume for tets
inputStruct.minRegionMarker=2; %Minimum region marker
inputStruct.modelName=modelName;

%%
% Mesh model using tetrahedral elements using tetGen (see:
% <http://wias-berlin.de/software/tetgen/>)

[meshOutput]=runTetGen(inputStruct); %Run tetGen

%% Visualizing mesh
meshView(meshOutput);

%%
% Access model element and patch data
F=meshOutput.faces;
V=meshOutput.nodes;
C=meshOutput.faceMaterialID;
E=meshOutput.elements;
CE=meshOutput.elementMaterialID;

Fb=meshOutput.facesBoundary;
Cb=meshOutput.boundaryMarker;

%%

cFigure;
xlabel('X','FontSize',fontSize); ylabel('Y','FontSize',fontSize); zlabel('Z','FontSize',fontSize)
hold on;

gpatch(Fb,V,Cb);

axisGeom;
colormap(plotColors)
drawnow;

%% Define boundary condition node sets

logicRigid=ismember(Cb,[5 6 7]);
bcSupportList=Fb(logicRigid,:);
bcSupportList=unique(bcSupportList(:));

%% Define pressure surface

F_pressure=fliplr(Fb(Cb==3,:));

%%
% Plot boundary condition nodes
cFigure;
xlabel('X','FontSize',fontSize); ylabel('Y','FontSize',fontSize); zlabel('Z','FontSize',fontSize)
hold on;

gpatch(Fb,V,Cb,'none',0.5);
plotV(V(bcSupportList,:),'k.','lineWidth',lineWidth1,'MarkerSize',markerSize1);
gpatch(F_pressure,V,0.5*ones(1,3),'k');

axisGeom;
colormap(plotColors);
drawnow;

%% Create element sets and material indices

%Create material indices
elementMaterialIndices=CE;
elementMaterialIndices(elementMaterialIndices==-3)=1;
elementMaterialIndices(elementMaterialIndices==-2)=2;

%Order material sets
E1=E(elementMaterialIndices==1,:); %Tissue material
E2=E(elementMaterialIndices==2,:); %Wrap material
E=[E1;E2];

%% Defining the FEBio input structure
% See also |febioStructTemplate| and |febioStruct2xml| and the FEBio user
% manual.

%Get a template with default settings
[febio_spec]=febioStructTemplate;

%febio_spec version
febio_spec.ATTR.version='3.0';

%Module section
febio_spec.Module.ATTR.type='solid';

%Create control structure for use by all steps
stepStruct.Control.analysis='STATIC';
stepStruct.Control.time_steps=numTimeSteps;
stepStruct.Control.step_size=1/numTimeSteps;
stepStruct.Control.solver.max_refs=max_refs;
stepStruct.Control.solver.max_ups=max_ups;
stepStruct.Control.time_stepper.dtmin=dtmin;
stepStruct.Control.time_stepper.dtmax=dtmax;
stepStruct.Control.time_stepper.max_retries=max_retries;
stepStruct.Control.time_stepper.opt_iter=opt_iter;

%Add template based default settings to proposed control section
[stepStruct.Control]=structComplete(stepStruct.Control,febio_spec.Control,1); %Complement provided with default if missing

%Remove control field (part of template) since step specific control sections are used
febio_spec=rmfield(febio_spec,'Control');

febio_spec.Step.step{1}.Control=stepStruct.Control;
febio_spec.Step.step{1}.ATTR.id=1;
febio_spec.Step.step{2}.Control=stepStruct.Control;
febio_spec.Step.step{2}.ATTR.id=2;

%Material section
materialName1='Normal_material';
febio_spec.Material.material{1}.ATTR.name=materialName1;
febio_spec.Material.material{1}.ATTR.id=1;
febio_spec.Material.material{1}.ATTR.type='Ogden unconstrained';
febio_spec.Material.material{1}.c1=c1;
febio_spec.Material.material{1}.m1=m1;
febio_spec.Material.material{1}.c2=c1;
febio_spec.Material.material{1}.m2=-m1;
febio_spec.Material.material{1}.cp=k;

materialName2='Multigen_material';
febio_spec.Material.material{2}.ATTR.name=materialName2;
febio_spec.Material.material{2}.ATTR.id=2;
febio_spec.Material.material{2}.ATTR.type='multigeneration';

febio_spec.Material.material{2}.generation{1}.ATTR.id=1;
febio_spec.Material.material{2}.generation{1}.start_time=0;
febio_spec.Material.material{2}.generation{1}.solid{1}.ATTR.type='Ogden unconstrained';
febio_spec.Material.material{2}.generation{1}.solid{1}.c1=c1_g(1);
febio_spec.Material.material{2}.generation{1}.solid{1}.m1=m1;
febio_spec.Material.material{2}.generation{1}.solid{1}.c2=c1_g(1);
febio_spec.Material.material{2}.generation{1}.solid{1}.m2=-m1;
febio_spec.Material.material{2}.generation{1}.solid{1}.cp=k_g(1);

febio_spec.Material.material{2}.generation{2}.ATTR.id=2;
febio_spec.Material.material{2}.generation{2}.start_time=1;
febio_spec.Material.material{2}.generation{2}.solid{1}.ATTR.type='Ogden unconstrained';
febio_spec.Material.material{2}.generation{2}.solid{1}.c1=c1_g(2);
febio_spec.Material.material{2}.generation{2}.solid{1}.m1=m1;
febio_spec.Material.material{2}.generation{2}.solid{1}.c2=c1_g(2);
febio_spec.Material.material{2}.generation{2}.solid{1}.m2=-m1;
febio_spec.Material.material{2}.generation{2}.solid{1}.cp=k_g(2);

%Mesh section
% -> Nodes
febio_spec.Mesh.Nodes{1}.ATTR.name='nodeSet_all'; %The node set name
febio_spec.Mesh.Nodes{1}.node.ATTR.id=(1:size(V,1))'; %The node id's
febio_spec.Mesh.Nodes{1}.node.VAL=V; %The nodel coordinates

% -> Elements
partName1='Part1';
febio_spec.Mesh.Elements{1}.ATTR.name=partName1; %Name of this part
febio_spec.Mesh.Elements{1}.ATTR.type='tet4'; %Element type
febio_spec.Mesh.Elements{1}.elem.ATTR.id=(1:1:size(E1,1))'; %Element id's
febio_spec.Mesh.Elements{1}.elem.VAL=E1;

partName2='Part2';
febio_spec.Mesh.Elements{2}.ATTR.name=partName2; %Name of this part
febio_spec.Mesh.Elements{2}.ATTR.type='tet4'; %Element type
febio_spec.Mesh.Elements{2}.elem.ATTR.id=size(E1,1)+(1:1:size(E2,1))'; %Element id's
febio_spec.Mesh.Elements{2}.elem.VAL=E2;

% -> NodeSets
nodeSetName1='bcSupportList';
febio_spec.Mesh.NodeSet{1}.ATTR.name=nodeSetName1;
febio_spec.Mesh.NodeSet{1}.node.ATTR.id=bcSupportList(:);

% -> Surfaces
surfaceName1='LoadedSurface';
febio_spec.Mesh.Surface{1}.ATTR.name=surfaceName1;
febio_spec.Mesh.Surface{1}.tri3.ATTR.id=(1:1:size(F_pressure,1))';
febio_spec.Mesh.Surface{1}.tri3.VAL=F_pressure;

%MeshDomains section
febio_spec.MeshDomains.SolidDomain{1}.ATTR.name=partName1;
febio_spec.MeshDomains.SolidDomain{1}.ATTR.mat=materialName1;

febio_spec.MeshDomains.SolidDomain{2}.ATTR.name=partName2;
febio_spec.MeshDomains.SolidDomain{2}.ATTR.mat=materialName2;

%Boundary condition section
% -> Fix boundary conditions
febio_spec.Boundary.bc{1}.ATTR.type='fix';
febio_spec.Boundary.bc{1}.ATTR.node_set=nodeSetName1;
febio_spec.Boundary.bc{1}.dofs='x,y,z';

%Loads section
% -> Surface load
febio_spec.Loads.surface_load{1}.ATTR.type='pressure';
febio_spec.Loads.surface_load{1}.ATTR.surface=surfaceName1;
febio_spec.Loads.surface_load{1}.pressure.ATTR.lc=1;
febio_spec.Loads.surface_load{1}.pressure.VAL=appliedPressure;
febio_spec.Loads.surface_load{1}.symmetric_stiffness=1;

%LoadData section
% -> load_controller
febio_spec.LoadData.load_controller{1}.ATTR.id=1;
febio_spec.LoadData.load_controller{1}.ATTR.type='loadcurve';
febio_spec.LoadData.load_controller{1}.interpolate='LINEAR';
febio_spec.LoadData.load_controller{1}.points.point.VAL=[0 0; 1 1; 2 0];

%Output section
% -> log file
febio_spec.Output.logfile.ATTR.file=febioLogFileName;
febio_spec.Output.logfile.node_data{1}.ATTR.file=febioLogFileName_disp;
febio_spec.Output.logfile.node_data{1}.ATTR.data='ux;uy;uz';
febio_spec.Output.logfile.node_data{1}.ATTR.delim=',';

febio_spec.Output.logfile.element_data{1}.ATTR.file=febioLogFileName_strainEnergy;
febio_spec.Output.logfile.element_data{1}.ATTR.data='sed';
febio_spec.Output.logfile.element_data{1}.ATTR.delim=',';

%% Quick viewing of the FEBio input file structure
% The |febView| function can be used to view the xml structure in a MATLAB
% figure window.

%%
% febView(febio_spec); %Viewing the febio file

%% Exporting the FEBio input file
% Exporting the febio_spec structure to an FEBio input file is done using
% the |febioStruct2xml| function.

febioStruct2xml(febio_spec,febioFebFileName); %Exporting to file and domNode

%% Running the FEBio analysis
% To run the analysis defined by the created FEBio input file the
% |runMonitorFEBio| function is used. The input for this function is a
% structure defining job settings e.g. the FEBio input file name. The
% optional output runFlag informs the user if the analysis was run
% succesfully.

febioAnalysis.run_filename=febioFebFileName; %The input file name
febioAnalysis.run_logname=febioLogFileName; %The name for the log file
febioAnalysis.disp_on=1; %Display information on the command window
febioAnalysis.runMode='external';%'internal';

[runFlag]=runMonitorFEBio(febioAnalysis);%START FEBio NOW!!!!!!!!

%% Import FEBio results

if runFlag==1 %i.e. a succesful run

    %%
    % Importing nodal displacements from a log file
    dataStruct=importFEBio_logfile(fullfile(savePath,febioLogFileName_disp),1,1);

    %Access data
    N_disp_mat=dataStruct.data; %Displacement
    timeVec=dataStruct.time; %Time

    %Create deformed coordinate set
    V_DEF=N_disp_mat+repmat(V,[1 1 size(N_disp_mat,3)]);

    %%
    % Importing element stress from a log file
    dataStruct=importFEBio_logfile(fullfile(savePath,febioLogFileName_strainEnergy),1,1);

    %Access data
    E_energy=dataStruct.data;

    %%
    % Plotting the simulated results using |anim8| to visualize and animate
    % deformations

    c1_plot=c1*ones(size(timeVec));
    cg_plot=c1_g(1)*ones(size(timeVec));
    cg_plot(timeVec>=1)=c1_g(2);

    [CV]=faceToVertexMeasure(E,V,E_energy(:,:,end));

    % Create basic view and store graphics handle to initiate animation
    hf=cFigure; %Open figure
    gtitle([febioFebFileNamePart,': Press play to animate']);

    subplot(1,2,1); hold on;
    title('Ogden parameter c_1');
    xlabel('Time'); ylabel('c_1');
    plot(timeVec,c1_plot,'b-','lineWidth',2);
    plot(timeVec,cg_plot,'r-','lineWidth',2);
    hp1=plot(timeVec(1),c1_plot(1),'b.','MarkerSize',50);
    hp2=plot(timeVec(1),cg_plot(1),'r.','MarkerSize',50);
    legend([hp1 hp2],'Material 1','Material 2');
    axis tight; axis square; set(gca,'fontsize',fontSize);
    grid on;

    subplot(1,2,2); hold on;
    hp3=gpatch(Fb,V_DEF(:,:,end),CV,'k',1); %Add graphics object to animate
    hp3.FaceColor='interp';
    colormap(gjet(250)); hc=colorbar;
    caxis([0 max(E_energy(:))]/10);
    axisGeom(gca,fontSize);
    axis(axisLim(V_DEF)); %Set axis limits statically
    axis manual;
    camlight headlight;
    drawnow;

    % Set up animation features
    animStruct.Time=timeVec; %The time vector
    for qt=1:1:size(N_disp_mat,3) %Loop over time increments
        DN=N_disp_mat(:,:,qt); %Current displacement
        DN_magnitude=sqrt(sum(DN.^2,2)); %Current displacement magnitude

        [CV]=faceToVertexMeasure(E,V,E_energy(:,:,qt));

        %Set entries in animation structure
        animStruct.Handles{qt}=[hp3 hp3 hp1 hp1 hp2 hp2]; %Handles of objects to animate
        animStruct.Props{qt}={'Vertices','CData','XData','YData','XData','YData'}; %Properties of objects to animate
        animStruct.Set{qt}={V_DEF(:,:,qt),CV,timeVec(qt),c1_plot(qt),timeVec(qt),cg_plot(qt)}; %Property values for to set in order to animate
    end
    anim8(hf,animStruct); %Initiate animation feature
    drawnow;

end

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [Lk]=logicErodeDilate(L,k,erodeDilateOption)

%Normalise structural kernel
k=k./sum(k(:));

%Convolve logic with kernel
M=double(L);
p=size(k)-1;
M_rep=zeros(size(M)+p);
M_rep(p(1)-1+(1:size(M,1)),p(2)-1+(1:size(M,2)),p(3)-1+(1:size(M,3)))=M;
Mk = convn(M_rep,k,'valid');
Mk=Mk./max(Mk(:)); %Scale max to 1

epsMax=max(eps(Mk(:)));

%Erode or dilate logic
switch erodeDilateOption
    case 'erode'
        Lk=(Mk>=(1-epsMax)); %Stayed nearly 1
    case 'dilate'
        Lk=Mk>=(0+epsMax); %Became higher than 0
end

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This script is used to generate the training and test samples which
% are saved as hdf5 format for the caffe input
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all;
clear all;
clc;

addpath('drtoolbox');
addpath('drtoolbox/gui');
addpath('drtoolbox/techniques');
addpath('caffe-master/matlab/hdf5creation');

%%%%%%%%%%%%%%%%%  Load dataset and groundtruth %%%%%%%%%%%%%%%%%%
%% for Indian Pines, the size of input is 25*25*3 %%%%%%%%%%%%%%%%
% load datasets/Indian_pines_gt.mat;
% load datasets/Indian_pines_corrected.mat;
% dir='caffe-master/data/indian_pines/';
% savepath_train = strcat(dir,'/train.h5');
% savepath_test  = strcat(dir,'/');
%
% pad_size=12;
% num_bands=3;
% input_size=2*pad_size+1;
% no_classes=16;
%
% img=indian_pines_corrected;
% [I_row,I_line,I_high] = size(img);
% img=reshape(img,I_row*I_line,I_high);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% for University of Pavia, the size of input is 23*23*5  %%%%%
load datasets/PaviaU.mat;
load datasets/PaviaU_gt.mat;
dir='caffe-master/data/paviau/';
savepath_train = strcat(dir,'/train.h5');
savepath_test  = strcat(dir,'/');

count_batch_test=5;
pad_size=11;
num_bands=5;
input_size=2*pad_size+1;
no_classes=9;

img=paviaU;
[I_row,I_line,I_high] = size(img);
img=reshape(img,I_row*I_line,I_high);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% for Salinas, the size of input is 27*27*10  %%%%%
% load datasets/Salinas_corrected.mat;
% load datasets/Salinas_gt.mat;
% dir='caffe-master/data/salinas/';
% savepath_train = strcat(dir,'/train.h5');
% savepath_test  = strcat(dir,'/');
%
% count_batch_test=5;
% pad_size=13;
% num_bands=10;
% input_size=2*pad_size+1;
% no_classes=16;
%
% img=salinas_corrected;
% [I_row,I_line,I_high] = size(img);
% img=reshape(img,I_row*I_line,I_high);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%  Dimension-reducing by PCA  %%%%%%%%%%%%%%%%%%%%%%%
im=img;
im=compute_mapping(im,'PCA',num_bands);
im = mat2gray(im);
im =reshape(im,[I_row,I_line,num_bands]);
[I_row,I_line,I_high] = size(im);

%%%%%  scale the image from -1 to 1
im=reshape(im,I_row*I_line,I_high);
[im ] = scale_func(im);
im =reshape(im,[I_row,I_line,I_high]);
%%%%%% extending the image %%%%%%%%
im_extend=padarray(im,[pad_size,pad_size],'symmetric');

[r,l,b]=size(im_extend);
im_extend=reshape(im_extend,[r*l,b]);
im_extend=im_extend-repmat(mean(im_extend),r*l,1);
im_extend=reshape(im_extend,[r,l,b]);

% [train_label,test_label,unlabeled_label,train_index,test_index,unlabeled_index] ...
%        = generating_index('indian_pines',indian_pines_gt,no_classes);

[train_label,test_label,unlabeled_label,train_index,test_index,unlabeled_index] ...
       = generating_index('paviau',paviaU_gt,no_classes);

% [train_label,test_label,unlabeled_label,train_index,test_index,unlabeled_index] ...
%        = generating_index('salinas',salinas_gt,no_classes);

[TRAIN_DATA, TRAIN_LABEL, TRAIN_INDEX] = write_save_h5('train', input_size,pad_size, ...
                     I_row, I_high,train_label,train_index, savepath_train, im_extend) ;

[TEST_DATA, TEST_LABEL, TEST_INDEX] = write_save_h5('test', input_size,pad_size, ...
                     I_row, I_high,test_label,test_index, savepath_test, im_extend) ;

save auxiliary_data/indian_pines/data.mat TRAIN_INDEX TEST_INDEX TEST_LABEL TRAIN_LABEL
%save auxiliary_data/paviau/data.mat TRAIN_INDEX TEST_INDEX TEST_LABEL TRAIN_LABEL
%save auxiliary_data/salinas/data.mat TRAIN_INDEX TEST_INDEX TEST_LABEL TRAIN_LABEL

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
%PLOTE Plot error curves
%
%   H = PLOTE(E,LINEWIDTH,S,FONTSIZE,OPTIONS)
%
% INPUT
%   E          Structure containing error curves (see e.g. CLEVAL)
%   LINEWIDTH  Line width, < 5 (default 2)
%   S          Plot strings
%   FONTSIZE   Font size, >= 5 (default 16)
%   OPTIONS    Character strings:
%              'nolegend' suppresses the legend plot
%              'errorbar' add errorbars to the plot
%              'noapperror' suppresses the apparent error plot
%
% OUTPUT
%   H          Array of graphics handles
%
% DESCRIPTION
% Various evaluation routines like CLEVAL return an error curves packed in a
% structure E. PLOTE uses the information stored in E to plot the curves. The
% remaining parameters may be given in an arbitrary order.
%
% E may contain the following fields (E.ERROR is obligatory):
%   E.ERROR    C x N matrix of error values for C methods at N points
%              (typically errors estimated on an independent test set)
%   E.XVALUES  C x N matrix of measurement points; if 1 x N, it is used for
%                all C curves
%   E.TITLE    the title of the plot
%   E.XLABEL   the label for the x-axis
%   E.YLABEL   the label for the y-axis
%   E.NAMES    a string array of C names used for creating a LEGEND
%   E.PLOT     the plot command in a string: 'plot', 'semilogx', 'semilogy'
%                or 'loglog'
%   E.STD      C x N matrix with standard deviations of the mean error values
%                which are plotted if ERRBAR == 1
%                Note that this is the st. dev. in the estimate of the mean
%                and not the std. dev. of the error itself.
%   E.APPERROR C x N matrix of error values for C methods at N points
%              (typically errors estimated on the training set)
%   E.APPSTD   C x N matrix with standard deviations of the mean
%              APPERROR values which are plotted if ERRBAR == 1
%
% These fields are automatically set by a series of commands like CLEVAL,
% CLEVALF, PRROC and REJECT.
%
% The legend generated by PLOTE can be removed by LEGEND OFF. A new legend
% may be created by the LEGEND command using the handles stored in H.
%
% E may be a cell array of structures. These structures are combined
% vertically, assuming multiple runs of the same method and
% horizontally, assuming different methods.
%
% EXAMPLES
% See PREX_CLEVAL
%
% SEE ALSO (PRTools Guide)
% CLEVAL, CLEVALF, PRROC, REJECT

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

% $Id: plote.m,v 1.5 2009/11/26 10:45:43 davidt Exp $

function handle = plote(varargin)

	% Set default parameter values.

	e         = [];
	s         = [];
	linewidth = 2;
	nolegend  = 0;
	fontsize  = 16;
	errbar    =  1;
	noapperror = 0;
	ss = char('k-','r-','b-','m-','k--','r--','b--','m--');
	ss = char(ss,'k-.','r-.','b-.','m-.','k:','r:','b:','m:');
	ss_app = char('k--','r--','b--','m--','k:','r:','b:','m:');
	ss_app = char(ss_app,'k--','r--','b--','m--','k-','r-','b-','m-');

	% The input is so flexible, that we have to do a lot of work...

	for j = 1:nargin
		p = varargin{j};
		if (isstruct(p)) | iscell(p)
			e = p;
		elseif (isstr(p)) & (strcmp(p,'errorbar') | strcmp(p,'ERRORBAR'))
			errbar = 1;
		elseif (isstr(p)) & (strcmp(p,'nolegend') | strcmp(p,'NOLEGEND'))
			nolegend = 1;
		elseif (isstr(p)) & (strcmp(p,'noapperror') | strcmp(p,'NOAPPERROR'))
			noapperror = 1;
		elseif (isstr(p))
			ss = p;
		elseif (length(p) == 1) & (p < 5)
			linewidth = p;
		elseif (length(p) == 1) & (p >= 5)
			fontsize = p;
		end
	end

	if iscell(e)
		if min(size(e)) > 1
			ee = cell(1,size(e,2));
			for j=1:size(e,2)
				ee{j} = vertcomb(e(:,j));
			end
			e = horzcomb(ee);
		elseif size(e,1) > 1
      for j=1:size(e,1)
        if j>1, figure; end
        plote(e{j});
      end
      showfigs
      return
			%e = vertcomb(e);
		elseif size(e,2) > 1
			e = horzcomb(e);
		else
			e = e{1};
		end
	end

	% Handle multiple plots

	if length(e) > 1
		names = [];
		hold_stat = ishold;
		h = [];
		ymax = 0;
    xmax = 0;
		for j = 1:length(e)
			if errbar & isfield(e,'std')
				if noapperror
					hh = plote(e(j),linewidth,ss(j,:),'nolegend','errorbar','noapperror');
				else
					hh = plote(e(j),linewidth,ss(j,:),'nolegend','errorbar');
				end
			else
				if noapperror
					hh = plote(e(j),linewidth,ss(j,:),'nolegend','noapperror');
				else
					hh = plote(e(j),linewidth,ss(j,:),'nolegend');
				end
			end
			V = axis;
      ymax = max(ymax,V(4));
      xmax = max(xmax,V(2));
			hold on
      if ~isfield(e(j),'names')
        e(j).names = ' ';
      end
			names = correct_(char(names,e(j).names));
			h = [h(:); hh(:)];
		end
		names(1,:) = [];
		V(4) = ymax;
    V(2) = xmax;
		axis(V);
		if ~nolegend
			legend(h,correct_(names),0);
		end
		if ~hold_stat, hold off; end
		if nargout > 0, handle = h; end
		return
	end

	% Check if we have the required data and data fields.

	if (isempty(e))
		error('Error structure not specified.')
	end

	if (~isfield(e,'error'))
		error('Input structure should contain the ''error''-field.');
	end

	n = size(e.error,1);

	if (~isfield(e,'xvalues'))
		e.xvalues = [1:length(e.error)];
	end
	if (size(e.xvalues,1) == 1)
		e.xvalues = repmat(e.xvalues,n,1);
	end

	if (isempty(s))
		if n > size(ss,1)
			nn = ceil(n/size(ss,1));
			ss = repmat(ss,nn,1);
			ss_app = repmat(ss_app,nn,1);
		end
		s = ss(1:n,:);
		s_app = ss_app(1:n,:);
	end
	if (size(s,1) == 1) & (n > 1)
		s = repmat(s,n,1);
		s_app = repmat(s_app,n,1);
	end
	if (size(s,1) < n)
		error('Insufficient number of plot strings.')
	end

	if (~isfield(e,'plot'))
		e.plot = 'plot';
	end

	if errbar & isfield(e,'std')
		ploterrorbar = 1;
	else
		ploterrorbar = 0;
	end

	plotapperror = (~noapperror && isfield(e,'apperror'));

	% We can now start making the plot.

	if ~ishold
		clf;
	end
	h = []; ha = []; % handles for true and apparent error
	for j = 1:n
		L = find(e.error(j,:) ~= NaN);
		if ploterrorbar
			hh = feval('errorbar',e.xvalues(j,L),e.error(j,L),e.std(j,L),deblank(s(j,:)));
		else
			hh = feval(e.plot,e.xvalues(j,L),e.error(j,L),deblank(s(j,:)));
		end
		set(hh,'linewidth',linewidth); hold on; h = [h hh(end)];

		% and the apparent error
		if plotapperror
			hh = errorbar(e.xvalues(j,L),e.apperror(j,L),e.appstd(j,L),...
			deblank(s_app(j,:)));
			ha = [ha hh(end)];
		end
	end

	% That was basically it, now we only have to beautify it.
	errmax = max(e.error(:));
	set(gca,'fontsize',fontsize);

	if (isfield(e,'xlabel')), xlabel(correct_(e.xlabel)); end
	if (isfield(e,'ylabel')), ylabel(correct_(e.ylabel)); end
	if (isfield(e,'title')),  title(correct_(e.title));   end
	if (isfield(e,'names')) & (~isempty(e.names) & (~nolegend))
		if plotapperror % take care for the legend in this case
			nrn = size(e.names,1);
      e.names = correct_(e.names);
			names = {};
			for j=1:nrn
				names{j} = ['\epsilon ' e.names(j,:)];
				%names{j} = ['true error ' e.names(j,:)];
			end
			for j=1:nrn
				names{nrn+j} = ['\epsilon_A ' e.names(j,:)];
				%names{nrn+j} = ['apparent error ' e.names(j,:)];
			end
			legend([h ha],strvcat(names),0);
		else
			legend(h,correct_(e.names),0);
		end
	end

	% A lot of work to make the scaling of the y-axis nice.

	if (errmax > 0.6)
		errmax = ceil(errmax*5)/5;
		yticks = [0:0.2:errmax];
	elseif (errmax > 0.3)
		errmax = ceil(errmax*10)/10;
		yticks = [0:0.1:errmax];
	elseif (errmax > 0.2)
		errmax = ceil(errmax*20)/20;
		yticks = [0:0.05:errmax];
	elseif (errmax > 0.1)
		errmax = ceil(errmax*30)/30;
		yticks = [0:0.03:errmax];
	elseif (errmax > 0.06)
		errmax = ceil(errmax*50)/50;
		yticks = [0:0.02:errmax];
	elseif (errmax > 0.03)
		errmax = ceil(errmax*100)/100;
		yticks = [0:0.01:errmax];
	else
		yticks = [0:errmax/3:errmax];
	end

		% atttempt to beautify plot
	if (e.xvalues(end) >= 2)
		%DXD
		%axis([e.xvalues(1)-1,e.xvalues(end)+1,0,errmax]);
		axis([min(min(e.xvalues)),max(max(e.xvalues)),0,errmax]);
	elseif (e.xvalues(1) == 0)
		axis([-0.003,e.xvalues(end),0,errmax]);
	end

	set(gca,'ytick',yticks);

	hold off; if (nargout > 0), handle = h; end;

	return

function e = vertcomb(e) % combine cell array
	e1 = e{1};
	for j=2:length(e);
		e2 = e{j};
		v = e1.n*(e1.n*(e1.std.^2) + e1.error.^2);
		v = v + e2.n*(e2.n*(e2.std.^2) + e2.error.^2);
		n = e1.n + e2.n;
		e1.error = (e1.n*e1.error + e2.n*e2.error)/n;
		e1.std = sqrt((v/n - e1.error.^2)/n);
		e1.n = n;
	end
	e = e1;
return

function ee = horzcomb(e) % combine cell array
  isx = false;
  isy = false;
  isd = false;
  for j=1:length(e)
    isx = (isx || strcmp(e{j}.plot,'semilogx'));
    isy = (isy || strcmp(e{j}.plot,'semilogy'));
    isd = (isd || strcmp(e{j}.plot,'loglog'));
  end
  if isx && isy, isd = true; end
	for j=1:length(e)
    if isd
      e{j}.plot = 'loglog';
    elseif isx
      e{j}.plot = 'semilogx';
    elseif isy
      e{j}.plot = 'semilogy';
    else
      e{j}.plot = 'plot';
    end
		ee(j) = e{j};
  end

function s = correct_(s) %correct string for underscores

  if size(s,1) > 1
    ss = cell(1,size(s,1));
    for j=1:size(s,1)
      ss{j} = correct_(s(j,:));
    end
    s = char(ss);
    return
  end

  J = find(s=='_');
  if ~isempty(J)
    J = fliplr(J(:)');
    for j=1:length(J)
      if s(J(j)-1) ~= '\'
        s = [s(1:J(j)-1) '\' s(J(j):end)];
      end
    end
  end

return

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
%%
% \documentclass[12pt]{article}
%
% \title{ODEbox: A Toolbox for Ordinary Differential Equations\\
% Example 2}
%
% \author{Matthew Harker and Paul O'Leary\\
% Institute for Automation\\
% University of Leoben\\
% A-8700 Leoben,
% Austria\\
% URL: automation.unileoben.ac.at\\
% \\
% Original: January 9, 2013\\
% $\copyright$ 2013\\
% \\
% Last Modified: \today}
%
%%
%\section{Initial Value Problem}
%
% This file demonstrated the solution of initial value problems using the
% DOPBox and ODEBox toolbaxes. The example being considered here is:
%
% \begin{equation}
%    2\,x^2\,\ddot{y} - x\,\dot{y} - 2 y = 0
%    \hspace{5mm}
%    \text{with}
%    \hspace{5mm}
%    y(1) = 5
%   \hspace{5mm}
%   \text{and}
%   \hspace{5mm}
%   \dot{y}(1) = 0
% \end{equation}
%
% The analytical solution to this problem is:
%
% \begin{equation}
%   y = x^2 + \frac{4}{\sqrt{x}}.
% \end{equation}
%
% This example has been taken from~\cite[pp. 394-395]{adams}.
%
% This M-file demonstrated the new approach to solving this problem and
% compared the result with a classical runga-kutta solution as provided by
% MATLAB.
%
%
%%
%\section{Clear up the Workspace}
close all;
clear all;
setUpGraphics(12) ;
%%
% Define the equation and solution as a strings for documentation
%
Eq = '$$2\,x^2\,\ddot{y} - x\,\dot{y} - 2 y = 0$$, with, $$y(1) = 5$$ and $$\dot{y}(1) = 0$$' ;
Sol = '$$y = x^2 + \frac{4}{\sqrt{x}}$$' ;
%
%%
% \section{Define the Matrix Linear Differential Operator}
%
% Define the number of points used in the solution
%
noPts = 73 ;
%
% Define the interval and use linearly spaced nodes.
%
xMin = 1 ;
xMax = 10 ;
x = linspace( xMin, xMax, noPts )';
%
% Setup the differentiating matrix with support length 13
%
ls = 13;
D = dopDiffLocal( x, ls, ls );
%
% Define the matrix C of constraints, i.e. the initial conditions
%
C = zeros(noPts,2) ;
C(1,1) = 1 ;
C(:,2) = D(1,:)' ;
%
d = [ 5 ; 0 ] ;
%
% Compute the linear differential operator.
%
L = 2*diag(x.^2)*D*D - diag(x)*D - 2*eye(noPts) ;
%%
% \section{Compute the solutions}
%
% Solve the algebraic system of equations
%
y = odeLSE(L,[],C,d) ;
%%
% For comparison solve using the MATLAB ode45 method.
% Runge-Kutta Solution:
%
[xm,Y] = ode45(@eulerODE02,[xMin xMax],[5 0]);
%
ym = Y(:,1);
%%
% the analytical solution for comparison.
%
f = inline('x.^2 + (4./(sqrt(x)))') ;
%
ya = f(x);
%%
% \section{Plot the results}
%
% for plot purposes solve beyond the ends
%
xp = linspace(xMin-1/4,xMax+1/4,1000) ;
yp = f(xp);
%
% Present the solutions
%
setUpGraphics(12)
FigureSize=[1 1 10 6];
set(0,'DefaultFigureUnits','centimeters');
set(0,'DefaultFigurePosition',FigureSize);
set(0,'DefaultFigurePaperUnits','centimeters');
set(0,'DefaultFigurePaperPosition',FigureSize);
MyAxesPosition=[0.13 0.15 0.86 0.82];
set(0,'DefaultaxesPosition',MyAxesPosition);
%
fig1 = figure;
plot( xp, yp, 'k' ) ;
hold on
plot( x, y, 'ko', 'MarkerFaceColor', 'w' ) ;
%
plot( xm, ym,'kv');
legend('Analytical','New','Runge-Kutta','Location','Northwest');
xlabel('$$x$$');
ylabel('$$y$$');
grid on;
range = axis;
plot( range(1:2), [0,0],'k');
axis([0.5,10.5,0,110]);
%\caption{Comparison of the analytical solution, the new numberical
% method and the Runga-Kutta solution.}
%%
%
setUpGraphics(12)
FigureSize=[1 1 10 6];
set(0,'DefaultFigureUnits','centimeters');
set(0,'DefaultFigurePosition',FigureSize);
set(0,'DefaultFigurePaperUnits','centimeters');
set(0,'DefaultFigurePaperPosition',FigureSize);
MyAxesPosition=[0.13 0.13 0.86 0.82];
set(0,'DefaultaxesPosition',MyAxesPosition);
%
fig2 = figure;
subplot(2,1,1)
plot( x, y - ya, 'k', 'MarkerFaceColor', 'w' ) ;
hold on
plot( x, y - ya, 'k.', 'MarkerFaceColor', 'w' ) ;
range = axis;
plot( range(1:2), [0,0],'k');
ylabel('$$e_{L}$$');
grid on;
%title(Sol);
%
subplot(2,1,2)
plot( xm, ym - f(xm), 'k', 'MarkerFaceColor', 'w' ) ;
hold on
plot( xm, ym - f(xm), 'k.', 'MarkerFaceColor', 'w' ) ;
hold on;
range = axis;
plot( range(1:2), [0,0],'k');
xlabel('$$x$$');
ylabel('$$e_{RK}$$');
grid on;
%
% \caption{Comparison the residual error for the new method (top) and
% for the Runga-Kutta method (Bottom). Note the residual is orders of magnitude smaller for the new method.}
%%
%\section{Save the figures to disk}
fileType = 'eps';
printFigure( fig1, 'ivpExp1Sol', fileType);
printFigure( fig2, 'ivpExp1Err', fileType);
%
%
%%
%% Define the Bibliography
%
% \bibliographystyle{IEEETran}
% \bibliography{odebib}

author: Sable

repo: mcbench-benchmarks


text:
% distance from radius of the head to the center
function [data,units] = compute_rhead_mm(trx,n)
larvae = trx.exp2flies{n};
numlarvae = numel(larvae);
data = cell(1,numlarvae);
for i = 1:numlarvae,
  larva = larvae(i);
  data{i} = bsxfun(@hypot,trx(larva).yhead_mm,trx(larva).xhead_mm);
end
units = parseunits('mm');

author: kristinbranson

repo: JAABA


text:
function p = set(p,varargin)
%PREAL/SET  The preal class SET method.

baseDims=7;

propertyArgIn=varargin;
while length(propertyArgIn)>=2,
    prop=propertyArgIn{1};
    val=propertyArgIn{2};
    propertyArgIn=propertyArgIn(3:end);
    switch prop
    case 'value'
        if isscalar(val)
            val=val*ones(size(p));
        end
        if ndims(val)~=ndims(p)
            error('Array dimensions must agree')
        end
        l=size(val)~=size(p);
        if any(l(:))
            error('Array dimensions must agree')
        end
        if isnumeric(val) & imag(val)==0
            for k=1:numel(val)
                p(k).value=val(k);
            end
        else
            error('Not a valid preal value')
        end
    case 'units'
        if isvector(val)&length(val)==baseDims&...
                isnumeric(val)&imag(val)==0
            for k=1:numel(p)
                p(k).units=val;
            end
        else
            error('Not valid preal units')
        end
    otherwise
        error('preal properties: value, units')
    end
end


text:
function [newSpikeTrain, newCalcium, newLL] = replaceSpike(oldSpikeTrain,oldCalcium,oldLL,filter,tau,obsCalcium,timetoRemove,indx,timetoAdd,Dt,A) %#codegen

% Replace a given spike with a new one in the existing spike train.

% Inputs:
% oldSpikeTrain:        current spike train
% oldCalcium:           current noiseless calcium trace
% oldLL:                current value of the log-likelihood function
% filters:              exponential rise and decay kernels for calcium transient
% tau:                  continuous time rise and decay time constants
% obsCalcium:           observed fluorescence trace
% timetoRemove:         time of the spike to be removed
% indx:                 place where the spike to be removed is in the existing spike train vector
% timetoAdd:            time of the spike to be added
% Dt:                   time-bin width
% A:                    spike amplitude

% Outputs:
% newSpikeTrain:        new vector of spike times
% newCalcium:           new noiseless calcium trace
% newLL:                new value of the log-likelihood function

% Author: Eftychios A. Pnevmatikakis and Josh Merel

    tau_h = tau(1);
    tau_d = tau(2);

    ef_h = filter{1,1};
    ef_d = filter{1,2};
    %ef_nh = filter{2,1};
    %ef_nd = filter{2,2};

    newSpikeTrain = oldSpikeTrain;
    newSpikeTrain(indx) = timetoAdd;

    %use infinite precision to scale the precomputed FIR approximation to the calcium transient
    wk_hr = A*exp((timetoRemove - Dt*ceil(timetoRemove/Dt))/tau_h);
    wk_dr = A*exp((timetoRemove - Dt*ceil(timetoRemove/Dt))/tau_d);

    wk_ha = A*exp((timetoAdd - Dt*ceil(timetoAdd/Dt))/tau_h);
    wk_da = A*exp((timetoAdd - Dt*ceil(timetoAdd/Dt))/tau_d);

    %%%%%%%%%%%%%%%%%
    %handle ef_h first
    newCalcium = oldCalcium;
    min_t = floor(min(timetoRemove,timetoAdd));
    if any(ef_h)
        tmp = 1+ (min_t:min((length(ef_h)+min_t-1),length(newCalcium)-1));
        if floor(timetoRemove) == floor(timetoAdd)
            wef_h = (wk_hr-wk_ha)*ef_h(1:length(tmp));
        elseif floor(timetoRemove) > floor(timetoAdd)
            wef_h = wk_hr*[zeros(1,min(length(tmp),floor(timetoRemove)-floor(timetoAdd))),ef_h(1:length(tmp)-(floor(timetoRemove)-floor(timetoAdd)))] - wk_ha*ef_h(1:length(tmp));
        else
            wef_h = wk_hr*ef_h(1:length(tmp)) - wk_ha*[zeros(1,min(length(tmp),floor(timetoAdd)-floor(timetoRemove))),ef_h(1:length(tmp)-(floor(timetoAdd)-floor(timetoRemove)))];
        end
        newCalcium(tmp) = newCalcium(tmp) - wef_h;
    end

    %%%%%%%%%%%%%%%%%
    %handle ef_d next
    tmp = 1+ (min_t:min((length(ef_d)+min_t-1),length(newCalcium)-1));
    if floor(timetoRemove) == floor(timetoAdd)
        wef_d = (wk_dr-wk_da)*ef_d(1:length(tmp));
    elseif floor(timetoRemove) > floor(timetoAdd)
        wef_d = wk_dr*[zeros(1,min(length(tmp),floor(timetoRemove)-floor(timetoAdd))),ef_d(1:length(tmp)-(floor(timetoRemove)-floor(timetoAdd)))] - wk_da*ef_d(1:length(tmp));
    else
        wef_d = wk_dr*ef_d(1:length(tmp)) - wk_da*[zeros(1,min(length(tmp),floor(timetoAdd)-floor(timetoRemove))),ef_d(1:length(tmp)-(floor(timetoAdd)-floor(timetoRemove)))];
    end
    newCalcium(tmp) = newCalcium(tmp) - wef_d;

    obstemp = obsCalcium(tmp);
    newLL = oldLL - norm(newCalcium(tmp)-obstemp)^2 + norm(oldCalcium(tmp)-obstemp)^2; %+ 2*(newCalcium(tmp)-oldCalcium
    %%%%%%%%%%%%%%%%

author: flatironinstitute

repo: CaImAn-MATLAB


text:
function pass = test_initialConditions(pref)
%TEST_INITIALCONDITIONS   Checks conditions for IVP/FVP problems.

%% Setup
if ( nargin == 0 )
    pref = cheboppref();
end
dom = [0, 1];
% [TODO]: This used to be 1e-12 until we had to loosen the absTol of the ODE
% solvers. Once/if that gets resolved, can we tighten the tolerance again?
tol = 5e-9;

%% First order IVP
N = chebop(@(x,u) diff(u) + sin(u), dom);
N.lbc = @(u) u - .5;
u = N\0;
pass(1) = abs(u(dom(1)) - .5) < tol;

%% First order FVP
N = chebop(@(x,u) diff(u) + sin(u), dom);
N.rbc = @(u) u - .5;
u = N\0;
pass(2) = abs(u(dom(end)) - .5) < tol;

%% Second order IVP
N = chebop(@(x,u) diff(u, 2) + sin(u), dom);
N.lbc = @(u) [u - .5; diff(u) - 1.3];
u = N\0;
up = diff(u);
pass(3) = abs(u(dom(1)) - .5) + abs(up(dom(1)) - 1.3) < tol;

%% Second order FVP
N = chebop(@(x,u) diff(u, 2) + sin(u), dom);
N.rbc = @(u) [u - .5; diff(u) - 1.3];
u = N\0;
up = diff(u);
pass(4) = abs(u(dom(end)) - .5) + abs(up(dom(end)) - 1.3) < tol;

%% First order IVP, no x
N = chebop(@(u) diff(u) + sin(u), dom);
N.lbc = @(u) u - .5;
u = N\0;
pass(5) = abs(u(dom(1)) - .5) < tol;

%% First order FVP, no x
N = chebop(@(u) diff(u) + sin(u), dom);
N.rbc = @(u) u - .5;
u = N\0;
pass(6) = abs(u(dom(end)) - .5) < tol;

%% Second order IVP, no x
N = chebop(@(u) diff(u, 2) + sin(u), dom);
N.lbc = @(u) [u - .5; diff(u) - 1.3];
u = N\0;
up = diff(u);
pass(7) = abs(u(dom(1)) - .5) + abs(up(dom(1)) - 1.3) < tol;

%% Second order FVP, no x
N = chebop(@(u) diff(u, 2) + sin(u), dom);
N.rbc = @(u) [u - .5; diff(u) - 1.3];
u = N\0;
up = diff(u);
pass(8) = abs(u(dom(end)) - .5) + abs(up(dom(end)) - 1.3) < tol;

%% First order coupled IVP
N = chebop(@(x,u, v, w) [diff(u) + v; diff(v) + w; diff(w) + u], dom);
N.lbc = @(u,v,w) [u - .5; v - 1.5; w + 2.5];
uvw = N\[0; 1; 2];
u = uvw{1}; v = uvw{2}; w = uvw{3};
pass(9) = abs(u(dom(1)) - .5) + abs(v(dom(1)) - 1.5) + abs(w(dom(1)) + 2.5) < tol;

%% First order coupled FVP
N = chebop(@(x,u, v, w) [diff(u) + v; diff(v) + w; diff(w) + u], dom);
N.rbc = @(u,v,w) [u - .5; v - 1.5; w + 2.5];
uvw = N\[0; 1; 2];
u = uvw{1}; v = uvw{2}; w = uvw{3};
pass(10) = abs(u(dom(end)) - .5) + abs(v(dom(end)) - 1.5) + abs(w(dom(end)) + 2.5) < tol;

%% Third order coupled IVP
N = chebop(@(x,u, v, w) [diff(u, 3) + v; diff(v, 3) + w; diff(w, 3) + u], dom);
N.lbc = @(u,v,w) [u - .5; v - 1.5; w + 2.5; ...
    diff(u) - 1; diff(v) - 2; diff(w) + 1; ...
    diff(u,2); diff(v,2) - .5; diff(w,2) + .5];
uvw = N\[0; 1; 2];
u = uvw{1}; v = uvw{2}; w = uvw{3};
up = diff(u); vp = diff(v); wp = diff(w);
upp = diff(u, 2); vpp = diff(v, 2); wpp = diff(w, 2);
pass(11) = abs(u(dom(1)) - .5) + abs(v(dom(1)) - 1.5) + abs(w(dom(1)) + 2.5) + ...
    abs(up(dom(1)) - 1) + abs(vp(dom(1)) - 2) + abs(wp(dom(1)) + 1) + ...
    abs(upp(dom(1))) + abs(vpp(dom(1)) - .5) + abs(wpp(dom(1)) + .5) < 100*tol;

%% Third order coupled FVP
N = chebop(@(x,u, v, w) [diff(u, 3) + v; diff(v, 3) + w; diff(w, 3) + u], dom);
N.rbc = @(u,v,w) [u - .5; v - 1.5; w + 2.5; ...
    diff(u) - 1; diff(v) - 2; diff(w) + 1; ...
    diff(u,2); diff(v,2) - .5; diff(w,2) + .5];
uvw = N\[0; 1; 2];
u = uvw{1}; v = uvw{2}; w = uvw{3};
up = diff(u); vp = diff(v); wp = diff(w);
upp = diff(u, 2); vpp = diff(v, 2); wpp = diff(w, 2);
pass(12) = abs(u(dom(end)) - .5) + abs(v(dom(end)) - 1.5) + abs(w(dom(end)) + 2.5) + ...
    abs(up(dom(end)) - 1) + abs(vp(dom(end)) - 2) + abs(wp(dom(end)) + 1) + ...
    abs(upp(dom(end))) + abs(vpp(dom(end)) - .5) + abs(wpp(dom(end)) + .5) < 15*tol;

%% Incorrect number of conditions

% Too few conditions.
N = chebop(@(x,u) diff(u, 2) + sin(u), dom);
N.lbc = @(u) u;
try
    mldivide(N, 0);
catch ME
    pass(13) = strcmp(ME.identifier, 'CHEBFUN:CHEBOP:solveivp:numConditions');
end

% Too many conditions
N = chebop(@(x,u) diff(u) + sin(u), dom);
N.lbc = @(u) [u; diff(u)];
try
    mldivide(N, 0);
catch ME
    pass(14) = strcmp(ME.identifier, 'CHEBFUN:CHEBOP:solveivp:numConditions');
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function obj = modelAddNewFolders(obj,varargin)
	% Add folders to the class.
	% Biafra Ahanonu
	% started: 2014.12.22 (probably before)
	% changelog
		% 2019.10.09 [18:28:12] - use inputdlgcol to allow re-sizing of window
		% 2019.11.18 [15:15:47] - Add the ability for users to use GUI to add folders as alternative option.
		% 2020.01.16 [12:28:29] - Choosing how to enter files and manual enter list of files now uses uicontrol and figure to reduce number of pop-ups and increase flexibility.
		% 2020.05.07 [17:22:20] - Adding option to quickly add all the example downloaded folders.
		% 2021.01.24 [14:03:44] - Added support for direct input of method type, useful for command-line or unit testing.
		% 2021.02.02 [11:27:21] - 'Add CIAtah example folders.' now adds the absolute path to avoid path errors if user changes Matlab current directory.
		% 2021.06.01 [??15:43:11] - Add 'next' button to initial menu, in part to add MATLAB Online support.
		% 2021.08.10 [09:57:36] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package.
		% 2022.04.10 [09:23:32] - Add additional check for when user adds blank folder paths (NOT empty folders) and remove them.
		% 2022.07.05 [19:27:24] - Empty folder check update.

	import ciapkg.api.* % import CIAtah functions in ciapkg package API.

	%========================
	% Cell array of folders to add, particularly for GUI-less operations
	options.folderCellArray = {};
	% Str:
	options.inputMethod = '';
	% get options
	options = getOptions(options,varargin);
	% display(options)
	% unpack options into current workspace
	% fn=fieldnames(options);
	% for i=1:length(fn)
	% 	eval([fn{i} '=options.' fn{i} ';']);
	% end
	%========================
	try
		nExistingFolders = length(obj.inputFolders);
		if isempty(options.folderCellArray)
			sel = 0;
			if isempty(options.inputMethod)
				usrIdxChoiceStr = {...
					'manually enter folders to list',...
					'GUI select folders',...
					'Add CIAtah example folders.'};
				scnsize = get(0,'ScreenSize');
				try
					hFig = figure;
					uicontrol('Style','Text','String',['How to add folders to CIAtah?'],'Units','normalized','Position',[5 89 90 10]/100,'BackgroundColor','white','HorizontalAlignment','Left','FontWeight','bold');
					exitHandle = uicontrol('style','pushbutton','Units', 'normalized','position',[5 85 50 3]/100,'FontSize',9,'string','Click here to move to next screen','callback',@subfxnNextFig,'HorizontalAlignment','Left');
					hListbox = uicontrol(hFig, 'style','listbox','Units', 'normalized','position',[5,5,90,80]/100, 'string',usrIdxChoiceStr,'Value',1);
					set(hListbox,'Max',2,'Min',0);
					set(hListbox,'KeyPressFcn',@(src,evnt)onKeyPressRelease(evnt,'press',hFig))
					figure(hFig)
					uicontrol(hListbox)
					set(hFig, 'KeyPressFcn', @(source,eventdata) figure(hFig));
					uiwait(hFig)
				catch err
					disp(repmat('@',1,7))
					disp(getReport(err,'extended','hyperlinks','on'));
					disp(repmat('@',1,7))

					[sel, ok] = listdlg('ListString',usrIdxChoiceStr,'ListSize',[scnsize(3)*0.3 scnsize(4)*0.3],'Name','How to add folders to CIAtah?');
				end
				inputMethod = usrIdxChoiceStr{sel};
			else
				inputMethod = options.inputMethod;
			end

			switch inputMethod
				case 'Add CIAtah example folders.'
					disp('Adding example folders to path...')
					dataDir = ciapkg.getDirPkg('data');
					newFolderList = {...
						[dataDir filesep '2014_04_01_p203_m19_check01'],...
						[dataDir filesep 'batch' filesep '2014_08_05_p104_m19_PAV08'],...
						[dataDir filesep 'batch' filesep '2014_08_06_p104_m19_PAV09'],...
						[dataDir filesep 'batch' filesep '2014_08_07_p104_m19_PAV10'],...
						[dataDir filesep 'twoPhoton' filesep '2017_04_16_p485_m487_runningWheel02']...
					};
                    disp(newFolderList)
					nNewFolders = length(newFolderList);
					fileIdxArray = (nExistingFolders+1):(nExistingFolders+nNewFolders);
					nFolders = length(fileIdxArray);
					newFolderListCell = {};
					for thisFileNumIdx = 1:nFolders
						% strtrim(newFolderList(thisFileNumIdx,:))
						% class(strtrim(newFolderList(thisFileNumIdx,:)))
						% ciapkg.getDir() filesep
						% pathToAdd = [obj.defaultObjDir filesep newFolderList{thisFileNumIdx}];
						pathToAdd = [newFolderList{thisFileNumIdx}];
						newFolderListCell{thisFileNumIdx} = strtrim(pathToAdd);
					end
					disp(nFolders)
					disp(newFolderListCell)
				case 'manually enter folders to list'
					newFolderList = '';
					try

						hFig = figure;
						figure(hFig)

						% finishMethodHandle = uicontrol('style','pushbutton','Units', 'normalized','position',[1 94 38 2]/100,'FontSize',9,'string','Start selected method (or press enter)','BackgroundColor',[153 255 153]/255,'callback',@subfxnCloseFig);

						uicontrol('Style','Text','String',['Adding folders to CIAtah object.'],'Units','normalized','Position',[5 95 90 3]/100,'BackgroundColor','white','HorizontalAlignment','Left','FontWeight','bold');
						uicontrol('Style','Text','String',['One new line per folder path. Enter folder path WITHOUT any single/double quotation marks around the path.'],'Units','normalized','Position',[5 90 90 6]/100,'BackgroundColor','white','HorizontalAlignment','Left');
						exitHandle = uicontrol('style','pushbutton','Units', 'normalized','position',[5 85 50 3]/100,'FontSize',9,'string','Click here to finish','callback',@subfxnCloseFig,'HorizontalAlignment','Left');

						hListbox = uicontrol(hFig, 'style','edit','Units', 'normalized','position',[5,5,90,80]/100, 'string','','FontSize',obj.fontSizeGui,'HorizontalAlignment','left');
						set(hListbox,'Max',2,'Min',0);
						uiwait(hFig)
						% hListboxStruct.Value = hListbox.string;
					catch err
						disp(repmat('@',1,7))
						disp(getReport(err,'extended','hyperlinks','on'));
						disp(repmat('@',1,7))

						AddOpts.Resize='on';
						AddOpts.WindowStyle='normal';
						AddOpts.Interpreter='tex';
						% inputdlg
						newFolderList = inputdlgcol('One new line per folder path. Enter folder path WITHOUT any single/double quotation marks around the path.','Adding folders to CIAtah object.',[21 150],{''},AddOpts);
						if isempty(newFolderList)
							warning('No folders given. Please re-run modelAddNewFolders.')
							return
						end
						newFolderList = newFolderList{1,1};
					end
					% size(newFolderList)
					% class(newFolderList)

					% if obj.guiEnabled==1
					% 	scnsize = get(0,'ScreenSize');
					% 	[fileIdxArray, ok] = listdlg('ListString',obj.fileIDNameArray,'ListSize',[scnsize(3)*0.2 scnsize(4)*0.25],'Name','which folders to analyze?');
					% else
					% 	fileIdxArray = 1:length(obj.fileIDNameArray);
					% end
					nNewFolders = size(newFolderList,1);
					fileIdxArray = (nExistingFolders+1):(nExistingFolders+nNewFolders);
					nFolders = length(fileIdxArray);
					newFolderListCell = {};
					for thisFileNumIdx = 1:nFolders
						% strtrim(newFolderList(thisFileNumIdx,:))
						% class(strtrim(newFolderList(thisFileNumIdx,:)))
						newFolderListCell{thisFileNumIdx} = strtrim(newFolderList(thisFileNumIdx,:));
					end
				case 'GUI select folders'
					thisFileNumIdx = 1;
					newFolderListCell = {};
					pathToAdd = '';
					while ~isempty(pathToAdd)|length(newFolderListCell)<1
						try
							if ischar(pathToAdd)
								disp('Select a folder to add. Press cancel to stop adding folders.')
								pathToAdd = uigetdir(pathToAdd,'Select a folder to add. Press cancel to stop adding folders.');
								% If user cancels, do not add folder.
								if pathToAdd~=0
									fprintf('Adding folder #%d: %s.\n',thisFileNumIdx,pathToAdd);
									newFolderListCell{thisFileNumIdx} = pathToAdd;
									thisFileNumIdx = thisFileNumIdx+1;
								end
							else
								% Force exit
								pathToAdd = [];
							end
						catch err
							disp(repmat('@',1,7))
							disp(getReport(err,'extended','hyperlinks','on'));
							disp(repmat('@',1,7))
						end
					end
					nNewFolders = length(newFolderListCell);
					fileIdxArray = (nExistingFolders+1):(nExistingFolders+nNewFolders);
				otherwise
					% body
			end
		else
			newFolderListCell = options.folderCellArray;
			nNewFolders = length(newFolderListCell);
		end

		newFolderListCell
		disp('Checking and removing blank folders')
		emptyFolderIdx = cellfun(@isempty,newFolderListCell,'UniformOutput',false);
		emptyFolderIdx = cell2mat(emptyFolderIdx);
		if any(emptyFolderIdx)
			newFolderListCell = newFolderListCell(~emptyFolderIdx);
		end
		newFolderListCell
		nNewFolders = length(newFolderListCell);

		fileIdxArray = (nExistingFolders+1):(nExistingFolders+nNewFolders);
		% obj.foldersToAnalyze = fileIdxArray;
		nFolders = length(fileIdxArray);
		display(repmat('-',1,7))
		display('Existing folders:')
		cellfun(@display,obj.inputFolders)
		display(repmat('-',1,7))
		display(['Number new folders:' num2str(nFolders) ' | New folder indices: ' num2str(fileIdxArray)]);
		for thisFileNumIdx = 1:nFolders
			fileNum = fileIdxArray(thisFileNumIdx);
			obj.fileNum = fileNum;
			obj.inputFolders{obj.fileNum,1} = newFolderListCell{thisFileNumIdx};
			obj.dataPath{obj.fileNum,1} = newFolderListCell{thisFileNumIdx};
			% display(repmat('=',1,21))
			% display([num2str(fileNum) '/' num2str(nFolders) ': ' obj.fileIDNameArray{obj.fileNum}]);
		end
		display('Folders added:')
		cellfun(@display,obj.inputFolders((end-nFolders+1):end));
		display('Adding file info to class...')
		obj.modelGetFileInfo();
		% display('Getting model variables...')
		% obj.modelVarsFromFiles();
		% display('Running pipeline...')
		% obj.runPipeline();

		% Reset folders to analyze
		obj.foldersToAnalyze = [];
	catch err
		obj.foldersToAnalyze = [];
		display(repmat('@',1,7))
		disp(getReport(err,'extended','hyperlinks','on'));
		display(repmat('@',1,7))
	end
	function subfxnCloseFig(src,event)
		newFolderList = hListbox.String;
		close(hFig)
	end
	function subfxnNextFig(src,event)
		sel = hListbox.Value;
		close(hFig)
	end
	function onKeyPressRelease(evnt, pressRelease,hFig)
		% disp(evnt)
		% disp(pressRelease)
		if strcmp(evnt.Key,'return')
			sel = hListbox.Value;
			% hListboxStruct = struct(hListbox);
			close(hFig)
		end
		% If escape, close.
		if strcmp(evnt.Key,'escape')
			sel = 1;
			close(hFig)
		end
	end
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function [u,p,info] = stokessolver(A,B,f,g,u,p,elem,option)
%% STOKESSOLVER solvers for Stokes equations
%
% We compute the solution [u,p] of the Stokes problem:
%
%     [A B'][u] = [f]
%     [B 0 ][p] = [g]
%
% by various solvers:
%   * 'direct'  the built in direct solver \ (mldivide)
%   * 'uzawa'   inexact Uzawa's method
%   * 'mg'      multigrid based on LSC-DGS smoother
%   * 'diagprestokes' block diagonal preconditioned MINRES
%   * 'triprestokes'  block triangular preconditioned GMRES
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

ufreeDof = option.ufreeDof;
pDof = option.pDof;

switch option.solver
    case 'direct'
        t = cputime;
        Nu = length(f);
        Np = length(g);
        bigA = [A, B'; ...
                B, sparse(Np,Np)];
        bigF = [f; g];
        bigu = zeros(Nu+Np);
        bigFreeDof = 1:2*Nu+pDof;
        bigu(bigFreeDof) = bigA(bigFreeDof,bigFreeDof)\bigF(bigFreeDof);
        u = bigu(1:2*Nu);
        p = bigu(2*Nu+1:end);
        residual = norm(bigF - bigA*bigu);
        solverTime = cputime - t;
        info = struct('solverTime',solverTime,'itStep',0,'err',residual,'flag',2,'stopErr',residual);
    case 'mg'
        option.solver  = 'WCYCLE';
        [u,p,info] = mgstokes(A,B,f,g,u,p,elem,ufreeDof,option);
    case 'asmg'
        [u,p,info] = asmgstokes(A,B,f,g,u,p,node,elem,option.bdFlag,ufreeDof,option);
    case 'blkdiag'
end

repo: ifem


text:
%% Copyright (C) 2016-2022 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy is free software; you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published
%% by the Free Software Foundation; either version 3 of the License,
%% or (at your option) any later version.
%%
%% This software is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty
%% of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See
%% the GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public
%% License along with this software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod @@sym asec (@var{x})
%% Symbolic asec function.
%%
%% Example:
%% @example
%% @group
%% syms x
%% y = asec (x)
%%   @result{} y = (sym) asec(x)
%% @end group
%% @end example
%%
%% Note: this file is autogenerated: if you want to edit it, you might
%% want to make changes to 'generate_functions.py' instead.
%%
%% @end defmethod

function y = asec(x)
  if (nargin ~= 1)
    print_usage ();
  end
  y = elementwise_op ('asec', x);
end

%!error asec (sym(1), 2)
%!assert (isequaln (asec (sym(nan)), sym(nan)))

%!shared x, d
%! d = 2;
%! x = sym('2');

%!test
%! f1 = asec(x);
%! f2 = asec(d);
%! assert( abs(double(f1) - f2) < 1e-15 )

%!test
%! D = [d d; d d];
%! A = [x x; x x];
%! f1 = asec(A);
%! f2 = asec(D);
%! assert( all(all( abs(double(f1) - f2) < 1e-15 )))

%!test
%! % round trip
%! y = sym('y');
%! A = asec (d);
%! f = asec (y);
%! h = function_handle (f);
%! B = h (d);
%! assert (A, B, -eps)

repo: octsympy


text:
function laguerre_poly_coef_test ( )

%*****************************************************************************80
%
%% LAGUERRE_POLY_COEF_TEST tests LAGUERRE_POLY_COEF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 October 2008
%
%  Author:
%
%    John Burkardt
%
  n = 5;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'LAGUERRE_POLY_COEF_TEST\n' );
  fprintf ( 1, '  LAGUERRE_POLY_COEF determines the Laguerre \n' );
  fprintf ( 1, '  polynomial coefficients.\n' );

  c = laguerre_poly_coef ( n );

  for i = 0 : n
    fprintf ( 1, '\n' );
    fprintf ( 1, '  L(%d)\n', i );
    fprintf ( 1, '\n' );
    for j = i : -1 : 0
      if ( j == 0 )
        fprintf ( 1, '  %f\n', c(i+1,j+1) );
      elseif ( j == 1 )
        fprintf ( 1, '  %f * x\n', c(i+1,j+1) );
      else
        fprintf ( 1, '  %f * x^%d\n', c(i+1,j+1), j );
      end
    end
  end

  fact = 1.0;

  for i = 0 : n

    if ( 0 < i )
      fact = fact * i;
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Factorially scaled L(%d)\n', i );
    fprintf ( 1, '\n' );

    for j = i : -1 : 0
      if ( j == 0 )
        fprintf ( 1, '  %f\n', fact * c(i+1,j+1) );
      elseif ( j == 1 )
        fprintf ( 1, '  %f * x\n', fact * c(i+1,j+1) );
      else
        fprintf ( 1, '  %f * x^%d\n', fact * c(i+1,j+1), j );
      end
    end

  end

  return
end

repo: jburkardt-m