text:
classdef FileSelector < wt.abstract.BaseWidget &...
        wt.mixin.Enableable & wt.mixin.FontStyled & wt.mixin.Tooltipable &...
        wt.mixin.FieldColorable & wt.mixin.ButtonColorable
    % A file/folder selection control with browse button

    % Copyright 2020-2021 The MathWorks Inc.

    %% Public properties
    properties (AbortSet)

        % The current value shown. If a RootDirectory is used, this will be
        % a relative path to the root.
        Value (1,1) string

    end %properties

    properties (Dependent, SetAccess = protected)

        % Absolute path to the file. If RootDirectory is used, this will
        % show the full path combining the root and the Value property.
        FullPath (1,1) string

        % Indicates the current value is a valid file path
        ValueIsValidPath (1,1) logical

    end %properties

    properties (AbortSet)

        % Selection type: file or folder
        SelectionType (1,1) wt.enum.FileFolderState = wt.enum.FileFolderState.file

        % Optional root directory. If unspecified, Value uses an absolute
        % path (default). If specified, Value will show a relative path to
        % the root directory and Value must be beneath RootDirectory.
        RootDirectory (1,1) string

        % Indicates whether to show a dropdown of recent folders
        ShowHistory (1,1) matlab.lang.OnOffSwitchState = false

        % List of recently selected folders to display in dropdown
        History (:,1) string

    end %properties

    % These properties do not trigger the update method
    properties (AbortSet, UsedInUpdate = false)

        % Optional default directory to start in when Value does not exist
        % and user clicks the browse button.
        DefaultDirectory (1,1) string

    end %properties

    %% Events
    events (HasCallbackProperty, NotifyAccess = protected)

        % Triggered on value changed, has companion callback
        ValueChanged

    end %events

    %% Internal Properties
    properties ( Transient, NonCopyable, ...
            Access = {?wt.abstract.BaseWidget, ?wt.test.BaseWidgetTest} )

        % Button
        ButtonControl (1,1) matlab.ui.control.Button

        % Edit control or dropdown
        EditControl (1,1) matlab.ui.control.EditField
        DropdownControl (1,1) matlab.ui.control.DropDown

    end %properties

    %% Protected methods
    methods (Access = protected)

        function setup(obj)

            % Adjust default size
            obj.Position(3:4) = [200 25];

            % Call superclass setup first to establish the grid
            obj.setup@wt.abstract.BaseWidget();

            % Configure Grid
            obj.Grid.ColumnWidth = {'1x',25};
            obj.Grid.RowHeight = {'1x'};

            % Create the standard edit control
            obj.EditControl = matlab.ui.control.EditField(...
                "Parent",obj.Grid,...
                "ValueChangedFcn",@(h,e)obj.onTextChanged(e));

            % Create the optional dropdown control
            obj.DropdownControl = matlab.ui.control.DropDown(...
                "Parent",[],...
                "Editable",true,...
                "Value","",...
                "ValueChangedFcn",@(h,e)obj.onTextChanged(e));
            %"Visible",false,...
            %obj.DropdownControl.Layout.Column = 1;

            % Create Button
            obj.ButtonControl = matlab.ui.control.Button(...
                "Parent",obj.Grid,...
                "Text","",...
                "ButtonPushedFcn",@(h,e)obj.onButtonPushed(e));
            obj.updateButtonIcon();

            % Update the internal component lists
            obj.FontStyledComponents = [obj.EditControl, obj.DropdownControl];
            obj.FieldColorableComponents = [obj.EditControl, obj.DropdownControl];
            obj.EnableableComponents = [obj.EditControl, obj.DropdownControl, obj.ButtonControl];
            obj.ButtonColorableComponents = obj.ButtonControl;
            obj.TooltipableComponents = [obj.EditControl, obj.DropdownControl, obj.ButtonControl];

        end %function

        function update(obj)

            % Is history being shown?
            if obj.ShowHistory
                % YES - Using dropdown control

                % Get the dropdown (history) items
                histItems = obj.History;

                % Current value must be in the list
                if ~isempty(obj.Value) && ~any(histItems == obj.Value)
                    histItems = vertcat(obj.Value,histItems);
                end

                % Update the items and value
                obj.DropdownControl.Items = histItems;
                obj.DropdownControl.Value = obj.Value;

            else
                % NO - Using edit control

                % Update the edit control text
                obj.EditControl.Value = obj.Value;

            end %if obj.ShowHistory

        end %function

        function updateButtonIcon(obj)

            % Update the button icon
            if obj.SelectionType == "file"
                obj.ButtonControl.Icon = "folder_file_24.png";
            else
                obj.ButtonControl.Icon = "folder_24.png";
            end

        end %function

        function updateControlType(obj)

            % Is history being shown? If so, update history and items
            if obj.ShowHistory
                % Using dropdown control

                if isempty(obj.DropdownControl.Parent)
                    obj.EditControl.Parent = [];
                    obj.DropdownControl.Parent = obj.Grid;
                    obj.DropdownControl.Layout.Column = 1;
                    obj.DropdownControl.Layout.Row = 1;
                end
            else
                % Using edit control

                if isempty(obj.EditControl.Parent)
                    obj.DropdownControl.Parent = [];
                    obj.EditControl.Parent = obj.Grid;
                    obj.EditControl.Layout.Column = 1;
                    obj.EditControl.Layout.Row = 1;
                end
            end %if obj.ShowHistory

        end %function

        function onButtonPushed(obj,~)
            % Triggered on button pushed

            % What folder should the prompt start at?
            if obj.ValueIsValidPath
                initialPath = obj.FullPath;
            elseif exist(obj.DefaultDirectory,"dir")
                initialPath = obj.DefaultDirectory;
            else
                initialPath = pwd;
            end

            % Prompt user for the path
            if obj.SelectionType == "file"
                [fileName,pathName] = uigetfile(initialPath,"Select a file");
            else
                pathName = uigetdir(initialPath,"Select a folder");
                fileName = "";
            end

            % Proceed if user didn't cancel
            if ~isequal(pathName,0)
                obj.setValueFromFullPath( fullfile(pathName,fileName) );
            end %if ~isequal(pathName,0)

        end %function

        function onTextChanged(obj,evt)
            % Triggered on text interaction

            % Prepare event data
            evtOut = wt.eventdata.PropertyChangedData('Value', evt.Value, obj.Value);

            % Store new result
            obj.Value = evt.Value;

            % Trigger event
            notify(obj,"ValueChanged",evtOut);

        end %function

        function addToHistory(obj,value)
            % Add new item to history and clean up

            % Add this to the top of the history
            histFiles = obj.History;
            histFiles = vertcat(value,histFiles);

            % Clean up the history
            histFiles = wt.utility.cleanPath(histFiles);
            histFiles = unique(histFiles,"stable");

            % Filter to valid paths
            if obj.SelectionType == "file"
                fcn = @(x)gt(exist(fullfile(obj.RootDirectory, x),"file"), 0);
            else
                fcn = @(x)eq(exist(fullfile(obj.RootDirectory, x),"dir"), 7);
            end
            isValidPath = arrayfun(fcn,histFiles);
            histFiles(~isValidPath) = [];

            % Limit history length
            histFiles(11:end) = [];

            % Store the result
            obj.History = histFiles;

        end %function

        function setValueFromFullPath(obj,fullPath)

            arguments
                obj (1,1)
                fullPath (1,1) string
            end %arguments

            fullPath = wt.utility.cleanPath(fullPath);

            if strlength(obj.RootDirectory) && strlength(fullPath)
                % Calculate the relative path for the value

                % Split the paths apart
                rootParts = strsplit(obj.RootDirectory,filesep);
                fullParts = strsplit(fullPath,filesep);

                % Find where the paths diverge
                idx = 1;
                smallestPath = min(numel(rootParts), numel(fullParts));
                while idx<=smallestPath && strcmpi(rootParts(idx),fullParts(idx))
                    idx = idx+1;
                end

                % Is the specified path outside of the root directory?
                numAbove = max(numel(rootParts) - idx + 1, 0);
                if numAbove>0

                    obj.throwError("Path '%s' is not within the root directory '%s'.",...
                        fullPath, obj.RootDirectory);

                else

                    % In case full path is above the RootPath, add ".." paths
                    parentPaths = string(repmat(['..' filesep],1,numAbove));

                    % Form the relative path
                    relPath = filesep + fullfile(parentPaths, fullParts{idx:end});

                    % What if paths are still the same?
                    if isempty(relPath)
                        relPath = "." + filesep;
                    end

                    % Prepare event data
                    evtOut = wt.eventdata.PropertyChangedData('Value',relPath, obj.Value);

                    % Store new result
                    obj.Value = relPath;

                    % Trigger event
                    notify(obj,"ValueChanged",evtOut);

                end

            else

                    % Prepare event data
                    evtOut = wt.eventdata.PropertyChangedData('Value',fullPath, obj.Value);

                    % Store new result
                    obj.Value = fullPath;

                    % Trigger event
                    notify(obj,"ValueChanged",evtOut);

            end %if

        end %function

    end % methods

    %% Accessors
    methods

        function set.Value(obj,value)
            value = wt.utility.cleanPath(value);
            obj.Value = value;
            obj.addToHistory(value)
        end

        function set.RootDirectory(obj,value)
            value = wt.utility.cleanPath(value);
            obj.RootDirectory = value;
            obj.addToHistory(value)
        end

        function set.SelectionType(obj,value)
            obj.SelectionType = value;
            obj.updateButtonIcon()
        end

        function set.ShowHistory(obj,value)
            obj.ShowHistory = value;
            obj.updateControlType()
        end

        function value = get.FullPath(obj)
            value = fullfile(obj.RootDirectory, obj.Value);
        end

        function value = get.ValueIsValidPath(obj)
            filePath = fullfile(obj.RootDirectory, obj.Value);
            value = ( obj.SelectionType == "file" && exist(filePath,"file") ) || ...
                ( obj.SelectionType == "folder" && exist(filePath,"dir") );
        end

    end % methods

end % classdef

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [nll,g,H,T] = LogisticLoss(w,X,y)
% w(feature,1)
% X(instance,feature)
% y(instance,1)

[n,p] = size(X);

Xw = X*w;
yXw = y.*Xw;

nll = sum(mylogsumexp([zeros(n,1) -yXw]));

if nargout > 1
    if nargout > 2
        sig = 1./(1+exp(-yXw));
        g = -X.'*(y.*(1-sig));
    else
        %g = -X.'*(y./(1+exp(yXw)));
        g = -(X.'*(y./(1+exp(yXw))));
    end
end

if nargout > 2
    H = X.'*diag(sparse(sig.*(1-sig)))*X;
end

if nargout > 3
    T = zeros(p,p,p);
    for j1 = 1:p
        for j2 = 1:p
            for j3 = 1:p
                T(j1,j2,j3) = sum(y(:).^3.*X(:,j1).*X(:,j2).*X(:,j3).*sig.*(1-sig).*(1-2*sig));
            end
        end
    end
end

author: zouchuhang

repo: LayoutNet


text:
function [hx, hy, hz] = ft_plot_ortho(dat, varargin)

% FT_PLOT_ORTHO plots 3 orthographic slices through a 3-D volume and interpolates if needed
%
% Use as
%   ft_plot_ortho(dat, ...)
% or
%   ft_plot_ortho(dat, mask, ...)
% where dat and mask are equal-sized 3-D arrays.
%
% Additional options should be specified in key-value pairs and can be
%   'style'        = string, 'subplot' or 'intersect' (default = 'subplot')
%   'orientation'  = 3x3 matrix specifying the directions orthogonal through the planes which will be plotted
%   'parents'      = (optional) 3-element vector containing the handles of the axes for the subplots (when style = 'subplot')
%   'surfhandle'   = (optional) 3-element vector containing the handles of the surfaces for each of the sublots (when style = 'subplot'). Parents and surfhandle are mutually exclusive
%   'update'       = (optional) 3-element boolean vector with the axes that should be updated (default = [true true true])
%
% The following options are supported and passed on to FT_PLOT_SLICE
%   'clim'                = [min max], lower and upper color limits
%   'transform'           = 4x4 homogeneous transformation matrix specifying the mapping from voxel space to the coordinate system in which the data are plotted
%   'location'            = 1x3 vector specifying the intersection point at which the three slices will be plotted. The coordinates should be expressed in the coordinate system of the data.
%   'datmask'             = 3D-matrix with the same size as the matrix dat, serving as opacitymap if the second input argument to the function contains a matrix, this will be used as the mask
%   'maskstyle'           = string, 'opacity' or 'colormix', defines the rendering
%   'background'          = needed when maskstyle is 'colormix', 3D-matrix with the same size as the data matrix, serving as grayscale image that provides the background
%   'interpmethod'        = string specifying the method for the interpolation, see INTERPN (default = 'nearest')
%   'colormap'            = string, see COLORMAP
%   'unit'                = string, can be 'm', 'cm' or 'mm (default is automatic)
%   'intersectmesh'       = triangulated mesh, see FT_PREPARE_MESH
%   'intersectcolor'      = string, color specification
%   'intersectlinestyle'  = string, line specification
%   'intersectlinewidth'  = number
%
% See also FT_PLOT_SLICE, FT_PLOT_MONTAGE, FT_SOURCEPLOT

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

% parse first input argument(s). it is either
% (dat, varargin)
% (dat, msk, varargin)
% (dat, [], varargin)
% this is done in ft_plot_slice

sellist = 1:numel(varargin);
if ~isempty(sellist)
  if isempty(varargin{1}) || isnumeric(varargin{1})
    sellist(1) = [];
  end
end

% get the optional input arguments
% other options such as location and transform are passed along to ft_plot_slice
style     = ft_getopt(varargin(sellist), 'style', 'subplot');
ori       = ft_getopt(varargin(sellist), 'orientation', eye(3));

if strcmp(style, 'subplot')
  parents    = ft_getopt(varargin(sellist), 'parents');
  surfhandle = ft_getopt(varargin(sellist), 'surfhandle');
  update     = ft_getopt(varargin(sellist), 'update', [true true true]);
  if ~isempty(surfhandle) && ~isempty(parents)
    ft_error('if specifying handles, you should either specify handles to the axes or to the surface objects, not both');
  end
end

if ~isa(dat, 'double')
  dat = cast(dat, 'double');
end

% determine the orientation key-value pair
keys = varargin(sellist(1:2:end));

sel  = find(strcmp('orientation', keys));
if isempty(sel)
  % add orientation key-value pair if it does not exist
  sel             = numel(keys)+1;
  varargin{2*sel-1} = 'orientation';
  varargin{2*sel} = [];
end

switch style
  case 'subplot'
    if isempty(parents) && isempty(surfhandle)
      Hx = subplot(2,2,1);
      Hy = subplot(2,2,2);
      Hz = subplot(2,2,4);
    elseif ~isempty(parents) && isempty(surfhandle)
      Hx = parents(1);
      Hy = parents(2);
      Hz = parents(3);
    elseif isempty(parents) && ~isempty(surfhandle)
      % determine the parents from the surface handle and use the
      % surfhandle for efficient visualization (overwriting existing data)
      if update(1), Hx = get(surfhandle(1), 'parent'); else Hx = []; end
      if update(2), Hy = get(surfhandle(2), 'parent'); else Hy = []; end
      if update(3), Hz = get(surfhandle(3), 'parent'); else Hz = []; end
    end

    if ~isempty(Hx)
      if ~isempty(surfhandle) && update(1)
        varargin(sellist) = ft_setopt(varargin(sellist), 'surfhandle', surfhandle(1));
      end
      % swap the first 2 dimensions because of meshgrid vs ndgrid issues
      varargin{2*sel} = ori(2,:);
      set(gcf,'currentaxes',Hx);
      hx = ft_plot_slice(dat, varargin{:});
      set(Hx, 'view', [0 0]); %, 'xlim', [0.5 size(dat,1)-0.5], 'zlim', [0.5 size(dat,3)-0.5]);
      if isempty(parents)
        % only change axis behavior if no parents are specified
        axis off
      end
    end

    if ~isempty(Hy)
      if ~isempty(surfhandle) && update(2)
        varargin(sellist) = ft_setopt(varargin(sellist), 'surfhandle', surfhandle(2));
      end
      varargin{2*sel} = ori(1,:);
      set(gcf,'currentaxes',Hy);
      hy = ft_plot_slice(dat, varargin{:});
      set(Hy, 'view', [90 0]); %, 'ylim', [0.5 size(dat,2)-0.5], 'zlim', [0.5 size(dat,3)-0.5]);
      if isempty(parents)
        % only change axis behavior if no parents are specified
        axis off
      end
    end

    if ~isempty(Hz)
      if ~isempty(surfhandle) && update(3)
        varargin(sellist) = ft_setopt(varargin(sellist), 'surfhandle', surfhandle(3));
      end
      varargin{2*sel} = ori(3,:);
      set(gcf,'currentaxes',Hz);
      hz = ft_plot_slice(dat, varargin{:});
      set(Hz, 'view', [0 90]); %, 'xlim', [0.5 size(dat,1)-0.5], 'ylim', [0.5 size(dat,2)-0.5]);
      if isempty(parents)
        % only change axis behavior if no parents are specified
        axis off
      end
    end

  case 'intersect'
    holdflag = ishold;
    if ~holdflag
      hold on
    end

    varargin{2*sel} = ori(1,:);
    hx = ft_plot_slice(dat, varargin{:});

    varargin{2*sel} = ori(2,:);
    hy = ft_plot_slice(dat, varargin{:});

    varargin{2*sel} = ori(3,:);
    hz = ft_plot_slice(dat, varargin{:});
    axis equal; axis tight; axis off;axis vis3d
    view(3);

    if ~holdflag
      hold off
    end

  otherwise
    ft_error('unsupported style %s', style);

end % switch style

% if strcmp(interactive, 'yes')
%   flag = 1;
%   while flag
%
%   end
% end

author: spm


text:
function [c,cd]=index0(s)
% INDEX :
% alternative expression of cluster structure
% of binary patterns based on various possible encodings
% named 'Cluster Indices'
%
% -----------------------------------------------------
% RETURNS :
% c : cluster vector, cd : cluster dim.
%
% Theophanes E. Raptis, DAT-NCSRD 2010
% http://cag.dat.demokritos.gr
% rtheo@dat.demokritos.gr
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    n = length(s);
    % identify cluster structure
    sflag = 1 - s(1);
    c = zeros(1,n);k = 0;
    for j=1:n
        if s(j)==0
           if sflag == 1 k = k+1; end
           sflag = 0;
           c(k) = c(k) - 1;
        else
           if sflag == 0 k = k+1; end
           sflag = 1;
           c(k) = c(k) + 1;
        end
    end
    % compute cluster functions
    m = length(find(c==0));
    cd = n-m; c = c(1:cd);
end

author: Sable

repo: mcbench-benchmarks


text:
function drawplat3d(nodez,lc,lm)
import rvctools.*
hold on
nodenum = size(nodez,1);
for j=1:nodenum
    if size(nodez,2)==12
        xc = [nodez(j,1:2),nodez(j,12),0,0,nodez(j,3)];
        xr = [nodez(j,1:2),0,0,0,nodez(j,3)]; a = nodez(j,4:7); th = nodez(j,8:11);
    else
        xc = [nodez(j,1:2),nodez(j,4),0,0,0];
        xr = [nodez(j,1:2),0,0,0,0]; a = nodez(j,3)*ones(1,4); th = zeros(1,4);
    end
    cub = pkgMechanics.RigidCuboid(1,xc,lc);
    patch('Vertices',cub.verticesStates.position','Faces',cub.faces,'FaceColor','b','FaceAlpha',0.5);

    xm = xr2m(xr,a,th,lc,lm);
    % build mobile platform cuboid
    for i=1:4
        plat(i) = pkgMechanics.RigidCuboid(1,xm(i,:),lm);
        patch('Vertices',plat(i).verticesStates.position','Faces',plat(i).faces,'FaceColor','y','FaceAlpha',0.5);
    end
end
hold off
end

repo: paper-simulation


text:
function value = r8_fall ( x, n )

%*****************************************************************************80
%
%% R8_FALL computes the falling factorial function [X]_N.
%
%  Discussion:
%
%    Note that the number of "injections" or 1-to-1 mappings from
%    a set of N elements to a set of M elements is [M]_N.
%
%    The number of permutations of N objects out of M is [M]_N.
%
%    Moreover, the Stirling numbers of the first kind can be used
%    to convert a falling factorial into a polynomial, as follows:
%
%      [X]_N = S^0_N + S^1_N * X + S^2_N * X^2 + ... + S^N_N X^N.
%
%  Formula:
%
%    [X]_N = X * ( X - 1 ) * ( X - 2 ) * ... * ( X - N + 1 ).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 June 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the argument of the falling factorial function.
%
%    Input, integer N, the order of the falling factorial function.
%    If N = 0, FALL = 1, if N = 1, FALL = X.  Note that if N is
%    negative, a "rising" factorial will be computed.
%
%    Output, real VALUE, the value of the falling factorial function.
%
  value = 1.0;

  arg = x;

  if ( 0 < n )

    for i = 1 : n
      value = value * arg;
      arg = arg - 1.0;
    end

  elseif ( n < 0 )

    for i = -1 : -1 : n
      value = value * arg;
      arg = arg + 1.0;
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ekg_RSlinB_am_ELF, ekg_RSlinB_bw_ELF, ekg_RSlinB_fm_ELF] = ELF_is(ekg_RSlinB_am, ekg_RSlinB_bw, ekg_RSlinB_fm, up)
%ELF eliminates very low frequencies from resampled respiratory

for rel_var_name_no = 1 : length(up.al.options.FMe)
%     temp = strfind(rel_var_names{rel_var_name_no},'_');
%     start_el = temp(1); clear temp
%     eval(['save_name = ''' curr_sig, up.paths.filenames.elim_vlf2, rel_var_names{rel_var_name_no}(start_el:end) ''';']); clear start_el
%     exist_log = check_exists(savepath, save_name);
%     if exist_log
%         continue
%     end

    %% Load relevant data
    rel_name = ['ekg_RSlinB_' up.al.options.FMe{rel_var_name_no} ];
    %loadpath = [up.paths.data_save_folder, num2str(subj), up.paths.filenames.int_respSigs];
    %load(loadpath, rel_name);
    eval(['old_data = ' rel_name ';']);

    %% Eliminate VLFs
    data.hpfilt.t = old_data.t;
    try
        data.hpfilt.v = elim_vlfs(old_data, up);
    catch
        % if there aren't enough points to use the filter, simply carry forward the previous data
        data.hpfilt.v = old_data.v;
    end
    data.hpfilt.fs = old_data.fs;
    %eval(['save_name = ' rel_name '_ELF']);
    save_name = [rel_name '_ELF'];
    eval([save_name ' = data.hpfilt;']);

    %% Save processed data
    %save_or_append_data
end

end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function [ns_RESULT, TimeStamp, Data, DataSize] = ns_GetEventData(hFile, EntityID, Index);

%ns_GetEventData   Retrieves event data by index
%
%   Usage:
%      [ns_RESULT, TimeStamp, Data, DataSize] =
%                                   ns_GetEventData(hFile, EntityID, Index)
%
%   Description:
%       Returns the data values from the file referenced by hFile and the
%       Event Entity EntityID.  The Event data entry specified by Index
%       is written to Data and the timestamp of the entry is returned to
%       TimeStamp.  Upon return of the function, the value at DataSize
%       contains the number of bytes actually written to Data.
%
%   Parameters:
%       hFile	        Handle/Indentification number to an open file.
%       EntityID	    Identification number of the entity in the data
%                       file.
%       Index	        The index number of the requested Event data item.
%
%   Return Values:
%       TimeStamp	    Variable that receives the timestamp of the Event
%                       data item.
%       Data	        Variable that receives the data for the Event entry.
%                       The format of Event data is specified by the member
%                       EventType in ns_EVENTINFO.
%       DataSize	    Variable that receives the actual number of bytes
%                       of data retrieved in the data buffer.
%       ns_RESULT   This function returns ns_OK if the file is successfully
%                   opened. Otherwise one of the following error codes is
%                   generated:
%
%                       ns_BADFILE	    Invalid file handle passed to
%                                       function
%                       ns_BADENTITY	Invalid or inappropriate entity
%                                       identifier specified
%                       ns_BADINDEX	    Invalid entity index specified
%                       ns_FILEERROR	File access or read error
%
%   Copyright (C) 2003 Neuroshare Project
%   Author: Almut Branner
%   Last modification: 8/11/2003

[ns_RESULT, TimeStamp, Data, DataSize] = mexprog(6, hFile, EntityID - 1, Index - 1);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [Z] = zca_whiten(X, Wsep, Wmix, mX)

Z = bsxfun(@minus, X, mX) * Wsep * Wmix;

author: kyunghyuncho

repo: deepmat


text:
function L = log_lik_complete(bnet, cases, clamped)
% LOG_LIK_COMPLETE Compute sum_m sum_i log P(x(i,m)| x(pi_i,m), theta_i) for a completely observed data set
% L = log_lik_complete(bnet, cases, clamped)
%
% If there is a missing data, you must use an inference engine.
% cases(i,m) is the value assigned to node i in case m.
% (If there are vector-valued nodes, cases should be a cell array.)
% clamped(i,m) = 1 if node i was set by intervention in case m (default: clamped = zeros)
% Clamped nodes contribute a factor of 1.0 to the likelihood.

if iscell(cases), usecell = 1; else usecell = 0; end

n = length(bnet.dag);
ncases = size(cases, 2);
if n ~= size(cases, 1)
  error('data should be of size nnodes * ncases');
end

if nargin < 3, clamped = zeros(n,ncases); end

L = 0;
for i=1:n
  ps = parents(bnet.dag, i);
  e = bnet.equiv_class(i);
  u = find(clamped(i,:)==0);
  L = L + log_prob_node(bnet.CPD{e}, cases(i,u), cases(ps,u));
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function DCM = spm_dcm_csd_data(DCM)
% gets cross-spectral density data-features using a VAR model
% FORMAT DCM = spm_dcm_csd_data(DCM)
% DCM    -  DCM structure
% requires
%
%    DCM.xY.Dfile        - name of data file
%    DCM.M.U             - channel subspace
%    DCM.options.trials  - trial to evaluate
%    DCM.options.Tdcm    - time limits
%    DCM.options.Fdcm    - frequency limits
%    DCM.options.D       - Down-sampling
%
% sets
%
%    DCM.xY.pst     - Peristimulus Time [ms] sampled
%    DCM.xY.dt      - sampling in seconds [s] (down-sampled)
%    DCM.xY.U       - channel subspace
%    DCM.xY.y       - cross spectral density over sources
%    DCM.xY.csd     - cross spectral density over sources
%    DCM.xY.It      - Indices of time bins
%    DCM.xY.Ic      - Indices of good channels
%    DCM.xY.Hz      - Frequency bins
%    DCM.xY.code    - trial codes evaluated
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_dcm_csd_data.m 6481 2015-06-16 17:01:47Z karl $

% Set defaults and Get D filename
%-------------------------------------------------------------------------
try
    Dfile = DCM.xY.Dfile;
catch
    errordlg('Please specify data and trials');
    error('')
end

% ensure spatial modes have been computed (see spm_dcm_csd)
%-------------------------------------------------------------------------
try
    DCM.M.U;
catch
    Nm      = DCM.options.Nmodes;
    DCM.M.U = spm_dcm_eeg_channelmodes(DCM.M.dipfit,Nm);
end

% load D
%--------------------------------------------------------------------------
try
    D = spm_eeg_load(Dfile);
catch
    try
        [p,f]        = fileparts(Dfile);
        D            = spm_eeg_load(f);
        DCM.xY.Dfile = fullfile(pwd,f);
    catch
        try
            [f,p]        = uigetfile('*.mat','please select data file');
            name         = fullfile(p,f);
            D            = spm_eeg_load(name);
            DCM.xY.Dfile = fullfile(name);
        catch
            warndlg([Dfile ' could not be found'])
            return
        end
    end
end

% indices of EEG channel (excluding bad channels)
%--------------------------------------------------------------------------
if ~isfield(DCM.xY, 'modality')
    [mod, list] = modality(D, 0, 1);

    if isequal(mod, 'Multimodal')
        qstr = 'Only one modality can be modelled at a time. Please select.';
        if numel(list) < 4

            % Nice looking dialog
            %--------------------------------------------------------------
            options = [];
            options.Default = list{1};
            options.Interpreter = 'none';
            DCM.xY.modality = questdlg(qstr, 'Select modality', list{:}, options);
        else

            % accomodate more buttons
            %--------------------------------------------------------------
            ind = menu(qstr, list);
            DCM.xY.modality = list{ind};
        end
    else
        DCM.xY.modality = mod;
    end
end

if ~isfield(DCM.xY, 'Ic')
    DCM.xY.Ic  = D.indchantype(DCM.xY.modality,'GOOD');
end

Ic        = DCM.xY.Ic;
Nm        = size(DCM.M.U,2);
DCM.xY.Ic = Ic;

% options
%--------------------------------------------------------------------------
try
    DT    = DCM.options.D;
catch
    DT    = 1;
end
try
    trial = DCM.options.trials;
catch
    trial = 1:D.nconditions;
end

% check data are not oversampled (< 4ms)
%--------------------------------------------------------------------------
if DT/D.fsample < 0.004
    DT            = ceil(0.004*D.fsample);
    DCM.options.D = DT;
end

% get peristimulus times
%--------------------------------------------------------------------------
try

    % time window and bins for modelling
    %----------------------------------------------------------------------
    DCM.xY.Time = time(D, [], 'ms');
    T1          = DCM.options.Tdcm(1);
    T2          = DCM.options.Tdcm(2);
    [i, T1]     = min(abs(DCM.xY.Time - T1));
    [i, T2]     = min(abs(DCM.xY.Time - T2));

    % Time [ms] of down-sampled data
    %----------------------------------------------------------------------
    It          = [T1:DT:T2]';               % indices - bins
    DCM.xY.pst  = DCM.xY.Time(It);           % PST
    DCM.xY.It   = It;                        % Indices of time bins
    DCM.xY.dt   = DT/D.fsample;              % sampling in seconds
    Nb          = length(It);                % number of bins

catch
    errordlg('Please specify time window');
    error('')
end

% get frequency range
%--------------------------------------------------------------------------
try
    Hz1     = DCM.options.Fdcm(1);          % lower frequency
    Hz2     = DCM.options.Fdcm(2);          % upper frequency
catch
    pst     = DCM.xY.pst(end) - DCM.xY.pst(1);
    Hz1     = max(ceil(2*1000/pst),4);
    if Hz1 < 8;
        Hz2 = 48;
    else
        Hz2 = 128;
    end
end

% Frequencies
%--------------------------------------------------------------------------
DCM.xY.Hz  = fix(Hz1:Hz2);             % Frequencies
Nf         = length(DCM.xY.Hz);        % number of frequencies
Ne         = length(trial);            % number of trial types

% Cross spectral density for each trial type
%==========================================================================
condlabels = D.condlist;               % condition or trial type labels
DCM.xY.csd = cell(1,Ne);               % CSD for each condition

w     = min(fix(2/DCM.xY.dt),Nb);      % window length (bins)
m     = 1;                             % retain principal mode
for i = 1:Ne;

    % trial indices
    %----------------------------------------------------------------------
    c = D.indtrial(condlabels(trial(i)), 'GOOD');
    fprintf('\nevaluating CSD for condition %i\n',i)

    % use only the first 512 trial
    %----------------------------------------------------------------------
    try c = c(1:512); end
    Nt    = length(c);

    % Get data
    %----------------------------------------------------------------------
    Nw    = max(8*(fix(Nb/w) - 1),1);
    K     = zeros(Nf*Nm*Nm,Nw);
    for k = 1:Nw
        P     = zeros(Nf,Nm,Nm);
        for j = 1:Nt
            Iw  = It((1:w) + fix((k - 1)*w/8));
            Y   = full(double(D(Ic,Iw,c(j))'*DCM.M.U));
            mar = spm_mar(Y,8);
            mar = spm_mar_spectra(mar,DCM.xY.Hz,1/DCM.xY.dt);
            P   = P + mar.P;
        end

        % store
        %------------------------------------------------------------------
        K(:,k) = spm_vec(P/Nt);
    end

    % retain principal eigenmode
    %----------------------------------------------------------------------
    [u s v]       = spm_svd(K,1);
    P             = u(:,m)*s(m,m)*mean(v(:,m));
    DCM.xY.csd{i} = spm_unvec(P,mar.P);

end

% place cross-spectral density in xY.y
%==========================================================================
DCM.xY.y    = DCM.xY.csd;
DCM.xY.U    = DCM.M.U;
DCM.xY.code = condlabels(trial);


text:
%kkbitrshift 'Bitwise RIGHT SHIFT of Input 1 by (Input 2 or Constant) Bits'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros kbitrshift.pane file
%
% Parameters:
% InputFile: i1 'Input 1', required: 'First input data object'
% OutputFile: o 'Output', required: 'Resulting output data object'
% InputFile: i2 'Input 2', optional: 'Second input data object'
% Integer: real 'Constant', default: 0: 'Constant value'
%
% Example: o = kkbitrshift({i1, i2}, {'i1','';'o','';'i2','';'real',0})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% kbitrshift - Bitwise RIGHT SHIFT of Input 1 by (Input 2 or Constant) Bits
%
%  DESCRIPTION
% The \fBRight Shift" operator performs a bitwise right shift on each point in
% \fBInput 1" by the number of bits defined by either the corresponding
% point in \fBInput 2" or the \fBConstant\fP value, which
% ever is specified by the user.
% Zero's are shifted into the
% most significant bits, and, in effect, each data value is divided by two and the
% remainder is discarded.
%
%  "Data Type - Single Input" 5
% .cI $DATAMANIP/repos/shared/man/sections/value_type_1input
%
%  "Data Type - Two Input Objects" 5
% .cI $DATAMANIP/repos/shared/man/sections/value_type_2input
% The bitwise operators will not operate on float or complex data types.
%
%  "Map Data" 5
% The bitwise operations have not been written to be fully polymorphic yet.
% They does not
% check for map data, and will therefore always operate on the value data,
% even if a map exists.  This will most likely corrupt indexing into the map.
% In the case of a single input, it is recommended to use the "Copy to
% Value" (kcptoval) segment operator to temporarily move the map data into
% the value segment, run the bitwise operation on the data, then move it
% back to the map with the "Copy from Value" (kcpfromval) operator.
% When operating with two input objects, where at least one of the objects
% contains map data, data should be mapped prior to running the bitwise
% operator.  The "Map Data" operator (kmapdata) can be used to perform
% the mapping operation.
%
%  "Validity Mask - Single Input" 5
% .cI $DATAMANIP/repos/shared/man/sections/mask_1input
%
%  "Validity Mask - Two Input Objects" 5
% Masking has not been implemented yet for the bitwise operators.  Therefore,
% only the mask from the first input object will be transferred to the output.
%
%  "Input Objects of Different Sizes" 5
% .cI $DATAMANIP/repos/shared/man/sections/resize_2input
% The value used to pad the data when the input files are not the same size is
% zero.
%
%  "Explicit Location and Time Data - Single Input" 5
% .cI $DATAMANIP/repos/shared/man/sections/loc_and_time_1input
%
%  "Explicit Location and Time Data - Two Input Objects" 5
% The bitwise operations have not been extended to understand location and
% time data.  Therefore, only location and time data present in the first input
% object will be transferred to the output.
%
% Executing \fBRight Shift" runs the program "kbitwise\fP with the -rs flag.
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
% DATAMANIP::kbitwise
%
%  RESTRICTIONS
%
% THIS ROUTINE is still under development, and has not been modified to fully
% support the polymorphic data model.  See paragraphs above for discussions
% concerning the polymorphic data segments.
%
% All processing is currently performed as UNSIGNED long, so the operations
% may not work properly on negative values.
%
% The bitwise operators will not operate on float or complex data types.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kkbitrshift(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kkbitrshift(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i1', '__input';'o', '__output';'i2', '__input';'real', 0};
maxval={0,0,1,0};
minval={0,0,1,0};
istoggle=[0,0,1,1];
was_set=istoggle * 0;
paramtype={'InputFile','OutputFile','InputFile','Integer'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'kbitwise"  -rs'],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function im = IdealLowPass(im0,fc)
% fc is the circular cutoff frequency which is normalized to [0 1], that is,
% the highest radian frequency \pi of digital signals is mapped to 1.

[ir,ic,iz] = size(im0);
hr = (ir-1)/2;
hc = (ic-1)/2;
[x, y] = meshgrid(-hc:hc, -hr:hr);

mg = sqrt((x/hc).^2 + (y/hr).^2);
lp = double(mg <= fc);

IM = fftshift(fft2(double(im0)));
IP = zeros(size(IM));
for z = 1:iz
    IP(:,:,z) = IM(:,:,z) .* lp;
end
im = abs(ifft2(ifftshift(IP)));


text:
function affineTransformationMatrix = get_affine_matrix(this)
% Transforms geometry parameters into affine 4x4 matrix (T*R*Z*S)
% with T = Translation, R = Rotation, Z = Zooming (Scaling with Resolution)
%      S = Shearing
% i.e. performing operations in the order from right to left, in particular
% rotation before translation
%
%   Y = MrImageGeometry()
%   Y.get_affine_matrix(inputs)
%
% This is a method of class MrImageGeometry.
%
% IN
%
% OUT
%
% EXAMPLE
%   get_affine_matrix
%
%   See also MrImageGeometry tapas_uniqc_spm_matrix

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

P(1:3) = this.offcenter_mm;
P(4:6) = this.rotation_deg*pi/180;
P(7:9) = this.resolution_mm;
P(10:12) = this.shear;

affineTransformationMatrix = tapas_uniqc_spm_matrix(P);

author: translationalneuromodeling

repo: tapas


text:
 function x = read_zubal_attn(varargin)
%function x = read_zubal_attn(options)
%| read in zubal attenuation phantom from data directory,
%| and assign it attenuation coefficients in inverse mm units.
%| options
%|	'nx'		desired size
%|	'ny'
%|	'ddir'		data directory

if nargin==1 && streq(varargin{1}, 'test'), read_zubal_attn_test, return, end

if nargout == 0, ir_usage, end

arg.ddir = '';
arg.file = 'zubal,attn.raw';
arg.nx = 128;
arg.ny = [];
arg = vararg_pair(arg, varargin);
if isempty(arg.ny), arg.ny = arg.nx; end

x = read_zubal_attn_do(arg.ddir, arg.file, arg.nx, arg.ny);

function x = read_zubal_attn_do(ddir, file, nx, ny)

% guess .../data directory by looking parallel to '.../transmission' directory
if ~isvar('ddir') || isempty(ddir)
	t = path_find_dir([filesep 'transmission']);
	ddir = strrep(t, 'transmission', 'data');
end

if ~exist(ddir, 'dir')
	warning(sprintf('cannot find data directory %s', ddir))
	error(sprintf('edit path in %s.m for your installation', mfilename))
end

file = [ddir filesep file];
if ~exist(file, 'file')
%	os_run(sprintf([ddir filesep 'do,attn,zubal %s'], file)) % here is how to create!
	error 'cannot find zubal phantom raw data'
end

fp = fopen(file, 'rb');
if (fp == -1), error 'open file', end
x = fread(fp, [128 128], 'uint8');
if fclose(fp), error 'close file', end

x = x(:,[(end-10):end 1:(end-11)]); % center it nicely

% assign attenuation coefficients in inverse mm units
mulist = [0.002 0.0096 0.0120];
for ii=1:length(mulist)
	x(x == ii) = mulist(ii);
end

x = ir_phantom_resize(x, nx, ny);

function read_zubal_attn_test
x = read_zubal_attn('nx', 120, 'ny', 90);
im(x, 'zubal phantom: transmission'), cbar

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [r,t,p]=spear(x,y)
%Syntax: [r,t,p]=spear(x,y)
%__________________________
%
% Spearman's rank correalation coefficient.
%
% r is the Spearman's rank correlation coefficient.
% t is the t-ratio of r.
% p is the corresponding p-value.
% x is the first data series (column).
% y is the second data series, a matrix which may contain one or multiple
%     columns.
%
%
% Reference:
% Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P.(1996):
% Numerical Recipes in C, Cambridge University Press. Page 641.
%
%
% Example:
% x = [1 2 3 3 3]';
% y = [1 2 2 4 3; rand(1,5)]';
% [r,t,p] = spear(x,y)
%
%
% Products Required:
% Statistics Toolbox
%
% Alexandros Leontitsis
% Department of Education
% University of Ioannina
% 45110- Dourouti
% Ioannina
% Greece
%
% University e-mail: me00743@cc.uoi.gr
% Lifetime e-mail: leoaleq@yahoo.com
% Homepage: http://www.geocities.com/CapeCanaveral/Lab/1421
%
% 3 Feb 2002.

% x and y must have equal number of rows
if size(x,1)~=size(y,1)
    error('x and y must have equal number of rows.');
end

% Find the data length
N = length(x);

% Get the ranks of x
R = crank(x)';

for i=1:size(y,2)

    % Get the ranks of y
    S = crank(y(:,i))';

    % Calculate the correlation coefficient
    r(i) = 1-6*sum((R-S).^2)/N/(N^2-1);

end

% Calculate the t statistic
if r == 1 | r == -1
    t = r*inf;
else
    t=r.*sqrt((N-2)./(1-r.^2));
end

% Calculate the p-values
p=2*(1-tcdf(abs(t),N-2));

function r=crank(x)
%Syntax: r=crank(x)
%__________________
%
% Assigns ranks on a data series x.
%
% r is the vector of the ranks
% x is the data series. It must be sorted.
%
%
% Reference:
% Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P.(1996):
% Numerical Recipes in C, Cambridge University Press. Page 642.
%
%
% Alexandros Leontitsis
% Department of Education
% University of Ioannina
% 45110- Dourouti
% Ioannina
% Greece
%
% University e-mail: me00743@cc.uoi.gr
% Lifetime e-mail: leoaleq@yahoo.com
% Homepage: http://www.geocities.com/CapeCanaveral/Lab/1421
%
% 3 Feb 2002.

u = unique(x);
[xs,z1] = sort(x);
[z1,z2] = sort(z1);
r = (1:length(x))';
r=r(z2);

for i=1:length(u)

    s=find(u(i)==x);

    r(s,1) = mean(r(s));

end


text:
% check whether the sequences are of the same length.
% The sequences are aligned from begining. If some sequence is shorter, it
% will be padded with a big negative number, i.e. -1e10, in the first dimension.
%
function [mask, variableLength] = CheckTrajectoryLength(data)

mask = permute(data(1,:,:), [2 3 1]) == -1e10;

variableLength = sum(mask(:));

end

author: singaxiong

repo: SignalGraph


text:
function title = p07_title ( )

%*****************************************************************************80
%
%% P07_TITLE returns a title for problem 7.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 February 2002
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, a title for the problem.
%
  title = 'The dying snake, ( x + sin(x) ) * e^(-x^2).';

  return
end


text:
function [del,obj] = setProj(obj,proj,projtype,insert,bbox)
    % [del,obj] = setProj(obj,proj,projtype,insert)
    % kjr generic function to parse projected space options
    % returns del flag to delete overlapping elements when plotting
    % global models.
    % if insert = 1, then automatically insert the global m_proj variables
    % into the msh obj. insert is 0 by default.

    if nargin < 4
        insert = 0;
    end

    if nargin < 5
        bbox = [];
    end

    % process bounds of mesh (or supply your own)
    if isempty(bbox)
        lon_mi = min(obj.p(:,1)); lon_ma = max(obj.p(:,1));
        lat_mi = min(obj.p(:,2)); lat_ma = max(obj.p(:,2));
    else
        lon_mi = bbox(1,1); lon_ma = bbox(1,2);
        lat_mi = bbox(2,1); lat_ma = bbox(2,2);
    end
    lat_mea = mean(obj.p(:,2)); lon_mea = mean(obj.p(:,1));

    % some defaults
    rad = 100; rot = 15;
    del = 0 ;

    % process the projtype as a varargin type
    I = find(strcmp(projtype,'long'));
    if ~isempty(I)
        if length(projtype{I+1}) == 1
            lon_mea = projtype{I+1};
        else
            lon_mi = projtype{I+1}(1); lon_ma = projtype{I+1}(2);
        end
    end
    I = find(strcmp(projtype,'lat'));
    if ~isempty(I)
        if length(projtype{I+1}) == 1
            lat_mea = projtype{I+1};
        else
            lat_mi = projtype{I+1}(1); lat_ma = projtype{I+1}(2);
        end
    end
    I = find(strcmp(projtype,'rad'));
    if ~isempty(I)
        rad = projtype{I+1};
    end
    I = find(strcmp(projtype,'rot'));
    if ~isempty(I)
        rot = projtype{I+1};
    end
    if proj == 0
        % normal geographic coordinates
        m_proj('equi','lat',[lat_mi lat_ma],'long',[lon_mi lon_ma]) ;
        del = 1;
    else
        if ~ischar(projtype)
            projtype = projtype{1};
        end
        projtype = lower(projtype);
        if ~isempty(regexp(projtype,'ste'))
            % Special treatment of Stereographic projection
            if lat_ma < 0
                % center Antarctica
                m_proj(projtype,'lat',-90,...
                      'long',0.5*(lon_mi+lon_ma),...
                      'radius',min(lat_ma+90,180),'rot',rot);
            else
                % center Arctic
                m_proj(projtype,'lat',90,...
                      'long',0.5*(lon_mi+lon_ma),...
                      'radius',min(90-lat_mi,180),'rot',rot);
            end
            m_proj('get') ;
        elseif  ~isempty(regexp(projtype,'ort')) || ...
                ~isempty(regexp(projtype,'gno')) || ...
                ~isempty(regexp(projtype,'azi')) || ...
                ~isempty(regexp(projtype,'sat'))
            m_proj(projtype,'lat',lat_mea,'long',lon_mea,...
                   'radius',rad,'rot',rot);
            m_proj('get') ;
            del = 1;
        elseif ~isempty(regexp(projtype,'obl'))
            % Oblique Mercator projection
            asp = (lon_ma-lon_mi)/(lat_ma - lat_mi);
            dir = 'hor';
            if asp > 1
                asp = 1./asp; dir = 'ver';
            end
            m_proj(projtype,'lon',[lon_mi lon_ma],...
                            'lat',[lat_mi lat_ma],...
                            'aspect',asp,'dir',dir) ;
            m_proj('get') ;
            del = 1;
        else
            % Cylindrical, Conic or Global type projections
            del = 1;
            m_proj(projtype,'lon',[lon_mi lon_ma],...
                            'lat',[lat_mi lat_ma]) ;
            m_proj('get') ;
        end
    end
    if insert
        global MAP_PROJECTION MAP_COORDS MAP_VAR_LIST
        obj.proj   = MAP_PROJECTION ;
        obj.coord  = MAP_COORDS ;
        obj.mapvar = MAP_VAR_LIST ;
    end
end

author: CHLNDDEV


text:
function [eventLog] = ma_executeCoast_nBody_goto_dt(dt, initialState, eventNum, forceModel, considerSoITransitions, soiSkipIds, massLoss, maxPropTime, events, contOnSoITrans, celBodyData)
%ma_executeCoast_nBody_goto_dt Summary of this function goes here
%   Detailed explanation goes here
    if(dt < 0)
        errorStr = ['Cannot coast ', num2str(dt), ' seconds.  Reverse time error, skipping event (any revs prior to this event may still be executed).'];
        addToExecutionErrors(errorStr, eventNum, initialState(8), celBodyData);

        eventLog = initialState;
        eventLog(:,13) = eventNum;
        return;
    elseif(dt == 0)
        eventLog = initialState;
        eventLog(:,13) = eventNum;
        return;
    end

    if(dt > maxPropTime)
        dt = maxPropTime;
    end

    ut = initialState(1);
	rVect = initialState(2:4);
    vVect = initialState(5:7);
    dryMass = initialState(9);
    fuelOxMass = initialState(10);
    monoMass = initialState(11);
    xenonMass = initialState(12);

	bodyID = initialState(8);
    bodyInfo = getBodyInfoByNumber(bodyID, celBodyData);
    topLevelBodyInfo = getTopLevelCentralBody(celBodyData);
    bodyInfoSun = topLevelBodyInfo;
    forceModel.gravModelBodyIds(end+1) = bodyInfoSun.id;
    forceModel.gravModelBodyIds = unique(forceModel.gravModelBodyIds);

    [rVect,vVect] = getAbsPositBetweenSpacecraftAndBody(ut, rVect', bodyInfo, bodyInfoSun, celBodyData, vVect');

    rhs = @(t,x) propRHS(t, x, bodyInfoSun, forceModel, massLoss, celBodyData);
    tspan = [ut, ut + dt]; %linspace(ut, ut + dt, numPts); %[ut, ut + dt];
    x0 = [-rVect', -vVect', dryMass, fuelOxMass, monoMass, xenonMass];
    eventFun = @(T,Y) ma_runNBodyEvents(T,Y, events);
    options = odeset('RelTol',1E-8, 'AbsTol',1E-8, 'Events',eventFun);
    [T,Y,TE,YE,IE] = ode113(rhs,tspan,x0,options);

    deltaT = T(end) - T(1);

    [rVect,vVect] = getAbsPositBetweenSpacecraftAndBody(T', Y(:,1:3)', bodyInfoSun, bodyInfo, celBodyData, Y(:,4:6)');

    eventLog = zeros(length(T), length(initialState));
    eventLog(:,1) = T;
    eventLog(:,2:4) = -rVect';
    eventLog(:,5:7) = -vVect';
    eventLog(:,8) = bodyID;
    eventLog(:,9) = Y(:,7);   %Dry Mass
    eventLog(:,10) = Y(:,8);  %Fuel/Ox Mass
    eventLog(:,11) = Y(:,9);  %Monoprop Mass
    eventLog(:,12) = Y(:,10); %Xenon Mass
    eventLog(:,13) = eventNum;

    if(~isempty(IE))
        [~,~,~,eventDescs] = eventFun(TE, YE');
        eDesc = eventDescs{IE};

        if(contains(eDesc, 'SoI Transition'))
            ut =  TE;
            rVect = YE(1:3);
            vVect = YE(4:6);

            if(contains(eDesc,'Up'))
                parentBodyInfo = bodyInfo.getParBodyInfo(celBodyData);
                [rVectNew, vVectNew] = convertRVVectOnUpwardsSoITransition(bodyInfo, celBodyData, ut, rVect, vVect);
                newBodyId = parentBodyInfo.id;
            elseif(contains(eDesc,'Down'))
                tokens = regexpi(eDesc,'- (\d+)\)','tokens');
                bodyId = str2double(tokens{1});
                childBodyInfo = getBodyInfoByNumber(bodyId, celBodyData);

                [rVectNew, vVectNew] = convertRVVectOnDownwardsSoITransition(childBodyInfo, celBodyData, ut, rVect, vVect);
                newBodyId = childBodyInfo.id;
            end

            soiTransInitState = [ut, rVectNew', vVectNew', newBodyId, eventLog(end,9:13)];
            eventLog(end+1,:) = soiTransInitState;

            if(contOnSoITrans == true)
                events = ma_updateEventFcnWithNewBodyInfo(events, getBodyInfoByNumber(newBodyId, celBodyData), celBodyData);

                soiTranseventLog = ma_executeCoast_nBody_goto_dt(dt-deltaT, soiTransInitState, eventNum, forceModel, considerSoITransitions, soiSkipIds, massLoss, maxPropTime-deltaT, events, contOnSoITrans, celBodyData);
                eventLog = [eventLog;soiTranseventLog];
            end
        end
    end
end

function xdot = propRHS(ut, x, cBodyInfo, forceModel, massLoss, celBodyData)
    gBodyIds = forceModel.gravModelBodyIds;

    rVectSC = x(1:3);
    vVectSC = x(4:6);

    gAccel = [0;0;0];
    for(i=1:length(gBodyIds)) %#ok<*NO4LP>
        gBodyId = gBodyIds(i);
        bodyOther = getBodyInfoByNumber(gBodyId, celBodyData);

        [rVect,~] = getAbsPositBetweenSpacecraftAndBody(ut, rVectSC, cBodyInfo, bodyOther, celBodyData);
        r = norm(rVect);

        gAccel = gAccel + (bodyOther.gm/r^3)*rVect; %the negative sign that is normally here is taken care of in the rVect already by virtue of how getAbsPositBetweenSpacecraftAndBody() works
    end

	xdot(1:3) = vVectSC;
    xdot(4:6) = gAccel;
    xdot(7) = 0;
    xdot(8) = 0;
    xdot(9) = 0;
    xdot(10) = 0;

	if(massLoss.use == 1 && ~isempty(massLoss.lossConvert))
        resRates = ma_getResRates(massLoss);

        xdot(8) = resRates(1);
        xdot(9) = resRates(2);
        xdot(10) = resRates(3);
	end

    xdot = xdot';
end


text:
%     NeuroSLAM System Copyright (C) 2018-2019
%     NeuroSLAM: A Brain inspired SLAM System for 3D Environments
%
%     Fangwen Yu (www.yufangwen.com), Jianga Shang, Youjian Hu, Michael Milford(www.michaelmilford.com)
%
%     The NeuroSLAM V1.0 (MATLAB) was developed based on the OpenRatSLAM (David et al. 2013).
%     The RatSLAM V0.3 (MATLAB) developed by David Ball, Michael Milford and Gordon Wyeth in 2008.
%
%     Reference:
%     Ball, David, Scott Heath, Janet Wiles, Gordon Wyeth, Peter Corke, and Michael Milford.
%     "OpenRatSLAM: an open source brain-based SLAM system." Autonomous Robots 34, no. 3 (2013): 149-176.
%
%     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/>.

% Jan 3, 2019
% groundTruthFile = 'C:\Dataset\3dratslam_groud_truth\01_ground_truth_data\03_3d_ground_truth_same_view_perspective_0918.txt';
odoMapFile = 'C:\NeuroSLAM_Datasets\03_NeuroSLAM_Experiments_Results\QUTCarparData\02_odo_map_ml.txt';
%     plot3(odomap_y*1*(-1), odomap_x*1+3, odomap_z*(1), '.b');
%     plot3((gt_x - gt_x(1))*17,(gt_y - gt_y(1))*23, (gt_z - gt_z(1))*(10), '.r');

plot_3d_odomap_qut(groundTruthFile, odoMapFile, 0.8, -0.8, 0.1, 0, 0.5, 0, 20,20, 20);

author: cognav

repo: NeuroSLAM

save_path: github-repos/MATLAB/cognav-NeuroSLAM


text:
function varargout = drawRect2(varargin)
%DRAWRECT2 Draw centered rectangle on the current axis
%
%   r = drawRect2(x, y, w, h) draw rectangle with width W and height H,
%   whose center is located at (x, y);
%
%   The four corners of rectangle are then :
%   (X-W/2, Y-H/2), (X+W/2, Y-H/2), (X+W/2, Y+H/2), (X-W/2, Y+H/2).
%
%   r = drawRect2(x, y, w, h, theta) also specifies orientation for
%   rectangle. Theta is given in radians.
%
%   r = drawRect2(coord) is the same as DRAWRECT2(X,Y,W,H), but all
%   parameters are packed into one array, whose dimensions is 4*1 or 5*1.
%
%   deprecated: use 'drawOrientedBox' instead
%
%   See Also :
%   drawRect
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 04/05/2004.
%

warning('polygons2d:deprecated', ...
    'This function is deprecated, use "drawOrientedBox" instead');

% default values
theta = 0;

% get entered values
if length(varargin)>3
    x = varargin{1};
    y = varargin{2};
    w = varargin{3};
    h = varargin{4};
    if length(varargin)>4
        theta = varargin{5};
    end
else
    coord = varargin{1};
    x = coord(:, 1);
    y = coord(:, 2);
    w = coord(:, 3);
    h = coord(:, 4);

    if length(coord)>4
        theta = coord(:, 5);
    else
        theta = zeros(size(x));
    end
end

% use only the half length of each rectanhle
w = w/2;
h = h/2;

hr = zeros(length(x), 1);
for i=1:length(x)
    tx = zeros(5, 1);
    ty = zeros(5, 1);

    tx(1) = x(i) - w(i)*cos(theta(i)) + h(i)*sin(theta(i));
    ty(1) = y(i) - w(i)*sin(theta(i)) - h(i)*cos(theta(i));

    tx(2) = x(i) + w(i)*cos(theta(i)) + h(i)*sin(theta(i));
    ty(2) = y(i) + w(i)*sin(theta(i)) - h(i)*cos(theta(i));

    tx(3) = x(i) + w(i)*cos(theta(i)) - h(i)*sin(theta(i));
    ty(3) = y(i) + w(i)*sin(theta(i)) + h(i)*cos(theta(i));

    tx(4) = x(i) - w(i)*cos(theta(i)) - h(i)*sin(theta(i));
    ty(4) = y(i) - w(i)*sin(theta(i)) + h(i)*cos(theta(i));

    tx(5) = x(i) - w(i)*cos(theta(i)) + h(i)*sin(theta(i));
    ty(5) = y(i) - w(i)*sin(theta(i)) - h(i)*cos(theta(i));

    hr(i) = line(tx, ty);
end

if nargout > 0
    varargout{1} = hr;
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips