text:
function ui = mrViewInfoPanel(ui,dockFlag, vis);
%
% ui = mrViewInfoPanel(ui,[dockFlag], [vis]);
%
% Adds a panel to the mrViewer UI with
% a display of information on the mr objects
% currently loaded.
%
% dockFlag: if 1, will attach panel to the
% mrViewer figure. Otherwise, will make a
% separate figure for the panel.
%
% vis: if 0, will initialize the panel to be hidden (default when
% starting up mrViewer); otherwise, leave it visible (as when
% docking/undocking the panel).
%
% ras, 07/05/05.
% ras, 07/15/05 -- now places the panel in a separate figure
% by default.
if ~exist('ui', 'var') | isempty(ui), ui = mrViewGet; end
if ~exist('dockFlag', 'var') | isempty(dockFlag), dockFlag = 1; end
if ~exist('vis', 'var') | isempty(vis),  vis = 0; end

bgColor = [.9 .9 .9];

if dockFlag==1
	ui.panels.info = mrvPanel('right', 0.3, ui.fig, 'normalized');
else
	hfig = figure('Name', sprintf('Info [%s]',ui.tag), ...
		'Units', 'normalized', ...
		'Position',[.88 .5 .12 .5], ...
		'MenuBar', 'none', ...
		'NumberTitle', 'off', ...
		'UserData',ui.tag);
	ui.panels.info = uipanel('Parent',hfig,'Units', 'normalized', ...
		'Position',[0 0 1 1]);

	% close request function for this figure: just toggle visibility
    crf = 'tmp = findobj(''Parent'', gcf, ''Type'', ''uipanel''); ';
    crf = [crf 'mrvPanelToggle(tmp); clear tmp;'];
    set(hfig, 'CloseRequestFcn', crf);
end

set(ui.panels.info,'BackgroundColor',bgColor, ...
	'ShadowColor',bgColor,'Title', 'Info');

% popup for selecting mr object about which to view info
ui.controls.infoPopup = uicontrol('Parent',ui.panels.info, ...
	'Style', 'popup', 'String',{sprintf('Base: %s',ui.mr.name)}, ...
	'Units', 'normalized', 'Position', [.05 .84 .8 .06], ...
	'Callback', 'mrViewSet([],''infoPanel'');');
uicontrol('Parent', ui.panels.info, 'Style', 'text', 'String', 'MR Object:', ...
	'HorizontalAlignment', 'center', 'FontWeight', 'bold', ...
	'BackgroundColor', bgColor, 'ForegroundColor', 'k', ...
	'Units', 'normalized', 'Position', [.05 .9 .8 .06]);

% listbox for showing info
ui.controls.infoListbox = uicontrol('Parent',ui.panels.info, ...
	'Style', 'listbox', 'String', infoText(ui.mr), 'FontSize', 9, ...
	'FontName', 'Helvetica', ...
	'Units', 'normalized', 'Position',[.05 .1 .8 .7]);

%%%%%hide the panel if requested
if vis==0
	mrvPanelToggle(ui.panels.info, 'off');
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


function ddata=mydiff(data,h)

L=length(data);
ddata=NaN(1,L);
for i=1:L;
    if i==1;
        ddata(1)=(-3*data(1)+4*data(2)-data(3))/(2*h);
    else
        if i==L;
            ddata(L)=(-3*data(L-2)+4*data(L-1)-data(L))/(2*h);
        else
            ddata(i)=(-data(i-1)+data(i+1))/(2*h);
        end
    end
end


text:
% Support only 2 dim array
classdef qdouble
   properties (Constant)
      mex = qdouble.mexSelector();
   end

   methods (Static)
      function func = mexSelector()
         if isarm()
            func = @qdoubleArmMex;
         else
            func = @qdoubleMex;
         end
      end
   end

   properties
      x;
   end

   methods
      function o = qdouble(a)
         if iscell(a) || isa(a, 'uint8')
            o.x = a;
         else
            o.x = qdouble.toMex(a);
         end
      end

      %% Sizes
      function [m, n] = size(a, dim)
         a = a.x;
         if (iscell(a))
            s = double(size(a{1}));
         else
            s = [size(a, 2) size(a, 3)];
         end

         if (nargin == 1)
            if (nargout == 2)
               m = s(1);
               n = s(2);
            else
               m = s;
            end
         else
            m = s(dim);
         end
      end

      function r = length(a)
         s = size(a);
         if (min(s) == 0)
            r = 0;
         else
            r = max(s);
         end
      end

      function r = numel(a)
         r = prod(size(a)); %#ok<*PSIZE>
      end

      function r = isscalar(a)
         r = (numel(a) == 1);
      end

      function r = isvector(a)
         r = (size(a,1) == 1 || size(a,2) == 1);
      end

      function r = isempty(a)
         r = (numel(a) == 0);
      end

      function r = isrow(a)
         r = (size(a,1) == 1);
      end

      function r = iscolumn(a)
         r = (size(a,2) == 1);
      end

      function r = ismatrix(~)
         r = true;
      end

      function r = end(a, k, n)
         if n == 2
            r = size(a,k);
         else
            r = size(a,1) * size(a,2);
         end
      end

      %% Conversions
      function disp(a)
         disp(double(a))
      end

      function r = double(a)
         r = qdouble.UnaryOp('double', a);
      end

      function r = logical(a)
         r = qdouble.UnaryOp('logical', a);
      end

      function r = sparse(a,b,c,d,e,f)
         if nargin == 1
            if ~issparse(a)
               A = logical(sparse(double(a)));
               s = struct('type', '()');
               s.subs = {find(A)};
               s = qdouble(subsref(a, s));
               s.x = reshape(s.x, [size(s.x,1) numel(s.x)/size(s.x,1)]);
               r = qdouble({A;s.x});
            else
               r = qdouble(a);
            end
         else
            assert(nargin == 3 || nargin == 5 || nargin == 6, "unsupported number of arguments for sparse");
            if nargin == 3
               i = a; j = b; s = c; m = max(i); n = max(j); nzmax = length(s);
               if isempty(i)
                  m = 0; n = 0;
               end
            elseif nargin == 5
               i = a; j = b; s = c; m = d; n = e; nzmax = length(s);
            elseif nargin == 6
               i = a; j = b; s = c; m = d; n = e; nzmax = f;
            end
            A = sparse(i, j, ones(nzmax,1,'logical'),m,n,nzmax);
            assert(length(s) == nzmax);
            s = qdouble(s(:));
            r = qdouble({A;s.x});
         end
      end

      function r = full(a)
         if issparse(a)
            A = a.x{1};
            s = a.x{2};
            r = zeros(size(s,1), size(A,1), size(A,2), 'uint8');
            idx = find(a.x{1});
            r(:, idx) = s;
            r = qdouble(r);
         else
            r = qdouble(a);
         end
      end

      function r = issparse(a)
         r = iscell(a.x);
      end

      function r = issymmetric(a)
         if size(a,1) ~= size(a,2)
            r = false;
         else
            r = full(all(a == a', 'all'));
         end
      end

      function r = diag(a)
         assert(nargin == 1, "Support only one parameters")
         if isvector(a)
            r = qdouble(sparse(a));
            if size(a,1) == 1
               r = r';
            end
            r.x{1} = diag(double(a));
            if ~issparse(a)
               r = full(r);
            end
         else
            k = min(size(a));
            if k == 0
               r = qdouble([]);
            else
               idx = (1:k) + size(a,1) * (0:(k-1));
               r =  subsref(a, struct('type', '()', 'subs', {{idx}}));
            end
            if issparse(a)
               r = sparse(r);
            end
            if size(a,2) < size(a,1)
               r = r';
            end
         end
      end

      function r = reshape(a, s)
         if issparse(a)
            a.x{1} = reshape(a.x{1}, s);
         else
            a.x = reshape(a.x, [size(a.x,1) s]);
         end
         r = qdouble(a.x);
      end

      function r = eps(~)
         r = qdouble(qdouble.mex('eps'));
      end

      %% Comparisons
      function r = lt(a, b)
         r = qdouble.BinaryOp('lt', a, b);
         if ((issparse(a) || issparse(b)) && ~issparse(r)), r = sparse(r); end
      end

      function r = gt(a, b)
         r = qdouble.BinaryOp('gt', a, b);
         if ((issparse(a) || issparse(b)) && ~issparse(r)), r = sparse(r); end
      end

      function r = le(a, b)
         r = ~gt(a,b);
      end

      function r = ge(a, b)
         r = ~lt(a,b);
      end

      function r = ne(a, b)
         r = qdouble.BinaryOp('ne', a, b);
         if ((issparse(a) || issparse(b)) && ~issparse(r)), r = sparse(r); end
      end

      function r = eq(a, b)
         r = ~ne(a,b);
      end

      function r = and(a, b)
         r = qdouble.BinaryOp('and', a, b);
      end

      function r = or(a, b)
         r = qdouble.BinaryOp('or', a, b);
      end

      function r = not(a)
         r = ~logical(a);
      end

      %% Arithmetic
      function r = plus(a, b)
         r = qdouble.BinaryOp('plus', a, b);
      end

      function r = minus(a, b)
         r = qdouble.BinaryOp('minus', a, b);
      end

      function r = times(a, b)
         r = qdouble.BinaryOp('times', a, b);
      end

      function r = power(a, b)
         assert(isscalar(b) && b == 2, "support only power(a, 2)");
         r = a .* a;
      end

      function r = rdivide(a, b)
         if issparse(b)
            b = full(b);
         end
         r = qdouble.BinaryOp('rdivide', a, b);

         if issparse(a) % To fix the bug regarding x / 0
            r2 = double(a) ./ double(b);
            r(~isfinite(r2)) = r2(~isfinite(r2));
         end
      end

      function r = ldivide(a, b)
         r = b./a;
      end

      function r = uminus(a)
         r = qdouble.UnaryOp('uminus', a);
      end

      function r = uplus(a)
         r = a;
      end

      function r = abs(a)
         r = qdouble.UnaryOp('abs', a);
      end

      function r = sqrt(a)
         r = qdouble.UnaryOp('sqrt', a);
      end

      function r = sum(varargin)
         r = qdouble.ReductionOp('sum', varargin{:});
      end

      function r = prod(varargin)
         r = qdouble.ReductionOp('prod', varargin{:});
      end

      function r = all(a, dim)
         if nargin == 1
            r = all(logical(a));
         else
            r = all(logical(a), dim);
         end
      end

      function r = any(a, dim)
         if nargin == 1
            r = any(logical(a));
         else
            r = any(logical(a), dim);
         end
      end

      function r = max(a, b, c)
         if nargin == 1
            r = qdouble.ReductionOp('max', a);
         elseif nargin == 3 && size(b,1) == 0 && size(b,2) == 0
            r = qdouble.ReductionOp('max', a, c);
         else
            assert(nargin == 2);
            r = qdouble.BinaryOp('max2', a, b);
            if issparse(a) || issparse(b)
               r = sparse(r);
            end
         end
      end

      function r = min(a, b, c)
         if nargin == 1
            r = qdouble.ReductionOp('min', a);
         elseif nargin == 3 && size(b,1) == 0 && size(b,2) == 0
            r = qdouble.ReductionOp('min', a, c);
         else
            assert(nargin == 2);
            r = qdouble.BinaryOp('min2', a, b);
            if issparse(a) || issparse(b)
               r = sparse(r);
            end
         end
      end

      function r = norm(a)
         assert(isvector(a), "norm supports only vectors.");
         r = sqrt(full(sum(a.^2)));
      end

      function r = nnz(a)
         r = full(sum(logical(a),'all'));
      end

      %% Access
      function [i,j,v] = find(a)
         if nargout == 1
            i = find(logical(a));
         elseif nargout == 2
            [i,j] = find(logical(a));
         else
            a = sparse(qdouble(a));
            [i,j] = find(a.x{1});
            v = qdouble(a.x{2});
            if size(a,1) == 1
               v = v';
            end
            if numel(v) == 0
               v = qdouble([]);
            end
         end
      end

      %% Matrix operations
      function r = ctranspose(a)
         r = transpose(a);
      end

      function r = transpose(a)
         r = qdouble.UnaryOp('transpose', a);
      end

      function r = mtimes(a, b)
         if isscalar(a) || isscalar(b)
            r = a .* b;
         else
            r = qdouble.BinaryOp('mtimes', a, b);
         end
      end

      function r = mldivide(a, b)
         assert(size(a,1) == size(b,1), 'Incompatible Size.');
         s = [size(a,2), size(b,2)];
         if size(a,2) == 0 || size(a,1) == 0 || size(b,2) == 0
            r = qdouble(zeros(s));
         elseif isscalar(a)
            r = reshape(b ./ a, s);
         else
            b_full = full(b); % sparsity usually does not help. Do this for simplicity
            r = qdouble.BinaryOp('mldivide', a, b_full);
         end
         if (issparse(a) && issparse(b))
            r = sparse(r);
         end
      end

      function r = mrdivide(a, b)
         r = (b'\a')';
      end

      function r = horzcat(varargin)
         allDense = true;
         for k = 1:length(varargin)
            allDense = allDense & ~issparse(varargin{k});
         end

         if (allDense)
            cat_input = cell(length(varargin), 1);
            for k = 1:length(varargin)
               cat_input{k} = varargin{k}.x;
            end
            r = qdouble(cat(3, cat_input{:}));
         else
            A_bool = [];
            m = []; v = [];
            for k = 1:length(varargin)
               A = sparse(qdouble(varargin{k}));
               if k == 1
                  m = size(A,1);
                  A_bool = A.x{1};
                  v = A.x{2};
               else
                  assert(m == size(A,1), 'Incompatible size.');
                  A_bool = [A_bool A.x{1}];
                  v = [v A.x{2}];
               end
            end
            r = qdouble({A_bool;v});
         end
      end

      function r = vertcat(varargin)
         for i = 1:length(varargin)
            varargin{i} = varargin{i}';
         end
         r = horzcat(varargin{:})';
      end

      function r = subsref(a,s)
         if strcmp(s(1).type, '()')
            if issparse(a)
               idx = subsref(qdouble.ToIndex(a),s);
               A = logical(idx);
               [~,~,idx] = find(idx);
               v = a.x{2};
               v = v(:, idx);
               r = qdouble({A;v});
            else
               idx = subsref(qdouble.ToIndex(a),s);
               r = qdouble(reshape(a.x(:,idx), [size(a.x,1) size(idx)]));
            end
         else
            r = builtin('subsref',a,s);
         end
      end

      function r = subsasgn(a,s,b)
         if strcmp(s(1).type, '()')
            [a_idx, a_v] = qdouble.ToIndex(a);
            [b_idx, b_v] = qdouble.ToIndex(b);
            b_idx = max(a_idx,[],'all') + b_idx;
            v = [a_v b_v];

            r_idx = subsasgn(a_idx,s,b_idx);

            if issparse(r_idx)
               A = logical(r_idx);
               [~,~,idx] = find(r_idx);
               v = v(:, idx);
               r = {A;v};
            else
               r = reshape(v(:, r_idx), [size(v,1) size(r_idx)]);
            end
            r = qdouble(r);
         else
            r = builtin('subsasgn',a,s,b);
         end
      end

      %% cholesky related
      function r = dissect(a)
         r = dissect(double(a));
      end

      function r = amd(a)
         r = amd(double(a));
      end

      function r = colamd(a)
         r = colamd(double(a));
      end

      function r = symbfact(a)
         r = symbfact(double(a));
      end

      function [r1, r2] = etree(a)
         [r1, r2] = etree(double(a));
      end

      function r = chol(a)
         r = qdouble.UnaryOp('chol', a);
      end
   end

   methods (Static)
      function r = ones(varargin)
         r = qdouble(ones(varargin{:}));
      end

      function r = zeros(varargin)
         r = qdouble(zeros(varargin{:}));
      end

      function r = eye(varargin)
         r = qdouble(eye(varargin{:}));
      end

      function r = rand(varargin)
         r = qdouble(rand(varargin{:}));
      end

      function r = randn(varargin)
         r = qdouble(randn(varargin{:}));
      end

      function r = randi(varargin)
         r = qdouble(randi(varargin{:}));
      end

      function r = toMex(a)
         if isa(a, 'qdouble')
            r = a.x;
         elseif isobject(a)
            r = a.x;
            r = qdouble.mex('toMex', r);
         else
            r = qdouble.mex('toMex', double(a));
         end
      end

      function r = UnaryOp(cmd, a)
         a = qdouble.toMex(a);
         r = qdouble.mex(cmd, a);

         if iscell(r) || isa(r, 'uint8')
            r = qdouble(r);
         end
      end

      function r = BinaryOp(cmd, a, b)
         a = qdouble.toMex(a);
         b = qdouble.toMex(b);
         r = qdouble.mex(cmd, a, b);

         if iscell(r) || isa(r, 'uint8')
            r = qdouble(r);
         end
      end

      function [r, v] = ToIndex(a)
         if nargout == 2 && ~isa(a, 'qdouble')
            a = qdouble(a);
         end

         if issparse(a)
            [i,j] = find(double(a));
            r = sparse(i,j,1:length(i),size(a,1),size(a,2));

            if nargout == 2
               v = a.x{2};
            end
         else
            r = reshape(1:numel(a),[size(a,1) size(a,2)]);

            if nargout == 2
               v = reshape(a.x, [size(a.x,1) numel(a.x)/size(a.x,1)]);
            end
         end
      end

      function r = ReductionOp(cmd, a, dim_)
         a_size = size(a);
         a = qdouble(a);

         if nargin == 2
            if (size(a,1) > 1)
               dim = 1;
            elseif (size(a,2) > 1)
               dim = 2;
            else
               dim = 1;
               a = subsref(a, struct('type', '()', 'subs', {{':'}}));
            end
         else
            assert(strcmp(dim_, 'all') || isscalar(dim_), 'dim can either be all or a scalar');
            if strcmp(dim_, 'all')
               dim = 1;
               a = subsref(a, struct('type', '()', 'subs', {{':'}}));
            else
               dim = dim_;
            end
         end

         if size(a,1) == 0
            switch cmd
               case 'sum'
                  r = qdouble(0.0 * ones(1, size(a,2)));
               case 'prod'
                  r = qdouble(1.0 * ones(1, size(a,2)));
               case 'max'
                  r = qdouble(ones(0, a_size(2)));
               case 'min'
                  r = qdouble(ones(0, a_size(2)));
            end
            if dim == 2
               r = r';
            end
         elseif dim == 1
            r = qdouble.mex(cmd, a.x);
            r = qdouble(r);
         elseif dim == 2
            a = a';
            r = qdouble.mex(cmd, a.x);
            r = qdouble(r)';
         end

         if issparse(a)
            r = sparse(r);
         end
      end
   end
end

author: opencobra

repo: cobratoolbox


text:
function [matchedObjMaps euclideanStruct] = displayMatchingObjs(inputImages,globalIDs,varargin)
    % Displays information on matching cells across sessions.
    % Biafra Ahanonu
    % started: 2014.01.03 [19:13:01]
    % inputs
        % inputImages - cell array {1, N} N = each trial of [x y nFilters] matrices containing each set of filters, e.g. {imageSet1, imageSet2,...}
        % globalIDs, [M N] matrix with M = number of global IDs and N = each trial. Each m,n pair specifies the index of that global obj m in the data of trial n. If .globalIDs(m,n)==0, means no match was found.
    % outputs
        %
    % changelog
        % 2017.01.14 [20:06:04] - support switched from [nSignals x y] to [x y nSignals]
        % 2021.08.08 [19:30:20] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package.
        % 2021.12.31 [18:59:24] - Updated suptitle to ciapkg.overloaded.suptitle
    % TODO
        %

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

    %========================
    options.inputSignals = [];
    %
    options.permuteImages = 0;
    %
    options.globalIDCoords = [];
    %
    options.sortGlobalIDs = 1;
    %
    options.fps = 20;
    % 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
        guiRow = 2;
        guiCol = 3;

        matchedObjMaps = [];
        euclideanStruct = struct;
        if options.permuteImages==1
            for imageNo=1:length(inputImages)
                % inputImages{imageNo} = permute(inputImages{imageNo}, [3 1 2]);
            end
        end

        if ~isempty(options.inputSignals)
            [~, ~] = openFigure(292, '');
            [nrows, ncols] = cellfun(@size, options.inputSignals);
            maxCols = max(ncols);
            signalSnrAll = NaN([length(options.inputSignals) maxCols]);
            for j=1:length(options.inputSignals)
                [signalSnr ~] = computeSignalSnr(options.inputSignals{j});
                signalSnrAll(j,1:length(signalSnr)) = signalSnr;
            end
            plot(signalSnrAll');
            % legend(292);
        end

        reverseStr = '';
        nGlobals = size(globalIDs,1);
        nObjPerGlobal = sum(globalIDs>0,2);
        if options.sortGlobalIDs==1
            % sort global IDs by those with most
            [nObjPerGlobalSorted sortedIdx] = sort(nObjPerGlobal,'descend');
            globalIDs = globalIDs(sortedIdx,:);
        end

        for gID = 1:size(globalIDs,2)
            cumProb(gID) = sum(sum(~(globalIDs==0),2)==gID)/size(globalIDs,1);
        end
        [~, ~] = openFigure(1930, '');
        plot(cumProb);box off;
        % set(gca,'yscale','log');
        % set(gca,'xscale','log');
        % [~, ~] = openFigure(8282, '');
        % for i=1:nGlobals
        globalNo = 1;
        saveData = 0;

        [nSignalsSessions, nPointsSessions] = cellfun(@size, options.inputSignals);
        while saveData==0
            nMatchedIDs = sum(globalIDs(globalNo,:)~=0);
            directionOfNextChoice = 1;
            % if globalNo==10
            %     saveData = 1;
            % end

            [~, ~] = openFigure(8282, '');
            for j=1:length(inputImages)
                iIdx = globalIDs(globalNo,j);
                if iIdx==0
                    nullImage = zeros(size(squeeze(inputImages{1}(:,:,1))));
                    nullImage(1,1) = 1;
                    groupImages(:,:,j) = nullImage;
                else
                    % display([num2str(j) ',' num2str(iIdx)])
                    groupImages(:,:,j) = squeeze(inputImages{j}(:,:,iIdx));
                end
                subplot(1,length(inputImages),j)
                imagesc(squeeze(groupImages(:,:,j)));
            end
            % display('')
            % pause

            [~, ~] = openFigure(1929, '');

            if ~isempty(options.inputSignals)
                % compute shuffled signalFeatures pairwise
                % if nMatchedIDs>2
                    nShuffles = 2;
                    reverseStr = '';
                    tmpGlobalIDs = globalIDs;
                    matchIDList = globalIDs(globalNo,:);
                    matchIDIdx = matchIDList~=0;
                    maxGlobalIDs = max(globalIDs,[],1);

                    matchIDNames = matchIDList(matchIDIdx);
                    matchIDNames = cellfun(@num2str,(mat2cell(matchIDNames(:),ones([1 length(matchIDNames)]))),'UniformOutput',false);

                    for shuffleNo = 1:nShuffles
                        % obtain a random ID within the range for each trial...clever clever...jk
                        for maxGlobalShuffleNo = 1:length(maxGlobalIDs)
                            if maxGlobalIDs(maxGlobalShuffleNo)==0
                                shuffledIDs(maxGlobalShuffleNo) = 0;
                            else
                                shuffledIDs(maxGlobalShuffleNo) = randsample(min(nSignalsSessions),1,true);
                            end
                            % if nSignalsSessions
                        end
                        % shuffledIDs = arrayfun(@(x) randsample(x,1,true), maxGlobalIDs);
                        shuffledIDs(~matchIDIdx) = 0;
                        tmpGlobalIDs(globalNo,:) = shuffledIDs(:)';
                        % tmpGlobalIDs
                        [~, matchedSignals] = getGlobalData(inputImages,tmpGlobalIDs,options.inputSignals,globalNo);
                        matchIdx = sum(matchedSignals,2)~=0;
                        [~, ~, signalFeaturesEuclideanShuffle, ~] = getSignalFeatures(matchedSignals(matchIdx,:));
                        shuffledSignalFeaturesEuclidean(shuffleNo) = signalFeaturesEuclideanShuffle;
                        reverseStr = cmdWaitbar(shuffleNo,nShuffles,reverseStr,'inputStr','shuffling euclidean distances','waitbarOn',0,'displayEvery',1);
                    end
                    shuffledSignalFeaturesEuclideanMean(globalNo,1) = nanmean(shuffledSignalFeaturesEuclidean);
                    shuffledSignalFeaturesEuclideanStd(globalNo,1) = nanstd(shuffledSignalFeaturesEuclidean);

                    % matchIdx = sum(matchedSignals,2)~=0;
                    [~, matchedSignals] = getGlobalData(inputImages,globalIDs,options.inputSignals,globalNo);
                    [avgSpikeTrace signalFeatures signalFeaturesEuclideanTmp signalFeaturesNorm] = getSignalFeatures(matchedSignals);
                    signalFeaturesEuclidean(globalNo,1) = signalFeaturesEuclideanTmp;
                % else
                %     shuffledSignalFeaturesEuclideanMean(globalNo,1) = NaN;
                %     shuffledSignalFeaturesEuclideanStd(globalNo,1) = NaN;
                %     signalFeaturesEuclidean(globalNo,1) = NaN;
                % end

                subplot(guiRow,guiCol,1:2)

                    imagesc(signalFeaturesNorm')
                    % add labels for each feature
                    numYTicks = size(signalFeatures,2);
                    L = get(gca,'YLim');
                    set(gca,'YTick',0:1:numYTicks*2)
                    set(gca,'YTickLabel',{'','signalSnr','numOfPeaks','slopeRatio','avgFwhm','avgPeakAmplitude'})

                    % display txt in the boxes
                    sizeX = size(signalFeaturesNorm);
                    idxX = 1:(sizeX(1)*sizeX(2));
                    [I,J] = ind2sub(sizeX,idxX);
                    % round values
                    txtVals = round(signalFeatures .* 10^3) ./ 10^3;
                    text(I,J,num2str(txtVals(:)),'HorizontalAlignment','center')
                    title([num2str(globalNo) '/' num2str(nGlobals) ' globalID list: ' num2str(globalIDs(globalNo,:))])
                    box off;
                    colormap(gca,[customColormap([])]);

                subplot(guiRow,guiCol,4)
                    % title('example trace')
                    peakROI = [-40:40];
                    % set(gca,'ColorOrder',copper(nMatchedIDs));
                    plot(peakROI, avgSpikeTrace');
                    box off;
                    % axis off;
                    % add in zero line
                    hold on;
                    xval = 0;
                    x=[xval,xval];
                    y=[min(avgSpikeTrace(:)) max(avgSpikeTrace(:))];
                    h = plot(x,y,'k'); box off;
                    % uistack(h,'bottom');

                    legend(matchIDNames,'Location','northeast')

                    % xval = 0;
                    % x=[xval,xval];
                    % y=[min(avgSpikeTrace(:)) max(avgSpikeTrace(:))];
                    % h = plot(x,y,'r'); box off;
                    % uistack(h,'bottom');

                    hold off;
                    title(['euclidean distance: ' num2str(signalFeaturesEuclidean(globalNo)) ' shuffled: ' num2str(shuffledSignalFeaturesEuclideanMean(globalNo))])
                % imagesc([0:0.1:1]); colorbar
                % title('reference for main plot')

                subplot(guiRow,guiCol,[5 6])
                    plotSignalsGraph(matchedSignals,'newAxisColorOrder','');
                    legend(matchIDNames)
            end

            % setup the next subplot
            subplot(guiRow,guiCol,3)
            [groupImages matchedSignals] = getGlobalData(inputImages,globalIDs,options.inputSignals,globalNo);
            % groupVector = 1:length(inputImages);
            % [groupedImages] = groupImagesByColor(thresholdImages(groupImages,'binary',1),[]);
            % matchedObjMaps(:,:,globalNo) = createObjMap(groupedImages);
            % imagesc(matchedObjMaps(:,:,globalNo));
            imagesc(squeeze(nansum(thresholdImages(groupImages,'binary',1,'threshold',0.4),3))/nMatchedIDs*100);
            colormap(gca,[0 0 0;customColormap([])]);
            title('heatmap of percent overlap object maps')
            colorbar
            axis equal tight;
            zoom on;
            box off

            ciapkg.overloaded.suptitle(['globalID ' num2str(globalNo) '/' num2str(nGlobals) '   f:finish    left/right: forward/back'])

            set(gcf,'currentch','3');
            % keyIn = get(gcf,'CurrentCharacter');
            keyIn = '3';
            while strcmp(keyIn,'3')
                keyIn = get(gcf,'CurrentCharacter');
                % if frameNo==frameNoMax
                %     frameNo = 1;
                % end
                pause(1/options.fps);
                % frameNo = frameNo + 1;
                % writeVideo(writerObj,getframe(mainFig));
            end

            reply = double(keyIn);
            figure(1929)
            set(gcf,'currentch','3');

            % [x y reply] = ginput(1);
            % reply = 29;
            if isequal(reply, 102)
                % return;
                saveData = 1;
            elseif isequal(reply, 28)
                % go back, left
                directionOfNextChoice=-1;
            elseif isequal(reply, 29)
                % go forward, right
                directionOfNextChoice=1;
            end
            % loop if user gets to either end
            globalNo=globalNo+directionOfNextChoice;
            if globalNo<=0 globalNo = nGlobals; end
            % if globalNo>nGlobals globalNo = 1; end
            if globalNo>nGlobals saveData = 1; end

            clear groupImages
            % if mod(i,7)==0|i==nGlobals
            %     reverseStr = cmdWaitbar(i,nGlobals,reverseStr,'inputStr','creating cell maps');
            % end

        end
        % playMovie(cellmap);
        euclideanStruct.shuffledSignalFeaturesEuclideanMean = shuffledSignalFeaturesEuclideanMean;
        euclideanStruct.shuffledSignalFeaturesEuclideanStd = shuffledSignalFeaturesEuclideanStd;
        euclideanStruct.signalFeaturesEuclidean = signalFeaturesEuclidean;

    catch err
        display(repmat('@',1,7))
        disp(getReport(err,'extended','hyperlinks','on'));
        display(repmat('@',1,7))
    end
end

function [groupImages matchedSignals] = getGlobalData(inputImages,globalIDs,inputSignals,globalNo)
    matchIDList = globalIDs(globalNo,:);
    matchIDIdx = matchIDList~=0;
    nMatchGlobalIDs = sum(matchIDIdx);
    if ~isempty(inputSignals)
        % get max length
        [nrows, ncols] = cellfun(@size, inputSignals);
        maxCols = max(ncols);
        matchedSignals = zeros(length(inputSignals),maxCols);
    end

    idxNo = 1;
    for j=1:length(inputImages)
        iIdx = globalIDs(globalNo,j);
        if iIdx==0
            nullImage = NaN(size(squeeze(inputImages{1}(:,:,1))));
            nullImage(1,1) = 1;
            groupImages(:,:,j) = nullImage;
        else
            % size(inputImages{j})
            % iIdx
            try
                groupImages(:,:,j) = squeeze(inputImages{j}(:,:,iIdx));
            catch
                display([num2str(j) ',' num2str(iIdx)])
            end
            if ~isempty(inputSignals)
                iSignal = inputSignals{j}(iIdx,:);
                matchedSignals(j,1:length(iSignal)) = iSignal;
            end
            idxNo = idxNo + 1;
        end
    end
end

function [avgSpikeTrace signalFeatures signalFeaturesEuclidean signalFeaturesNorm] = getSignalFeatures(matchedSignals)

    % get the SNR for traces
    [signalSnr ~] = computeSignalSnr(matchedSignals,'waitbarOn',0);
    % get the peak statistics
    [peakOutputStat] = computePeakStatistics(matchedSignals,'waitbarOn',0);
    slopeRatio = peakOutputStat.slopeRatio;
    avgFwhm = peakOutputStat.avgFwhm;
    avgPeakAmplitude = peakOutputStat.avgPeakAmplitude;
    avgSpikeTrace = peakOutputStat.avgSpikeTrace;
    % get the number of spikes
    [signalPeaks, signalPeakIdx] = computeSignalPeaks(matchedSignals,'makeSummaryPlots',0,'waitbarOn',0);
    numOfPeaks = sum(signalPeaks,2);
    signalFeatures = horzcat(signalSnr(:),...
        numOfPeaks(:),...
        slopeRatio(:),...
        avgFwhm(:),...
        avgPeakAmplitude(:));

    % signalFeatures(find(isnan(signalFeatures))) = 0;
    % normalize feature matrix
    X = signalFeatures;
    [m n] = size(signalFeatures);
    xMax = repmat(nanmax(X),m,1);
    xMin = repmat(nanmin(X),m,1);
    signalFeaturesNorm = (X - xMin)./(xMax - xMin);

    % compute pairwise distance
    % pdist(signalFeaturesNorm')
    % signalFeaturesEuclidean = nanmean(pdist(signalFeatures));
    signalFeaturesEuclidean = nanmean(pdist(signalFeaturesNorm));
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
classdef SymmetricFourthOrder3DTensor < FourthOrder3DTensor

    methods (Access = public)

        function obj = SymmetricFourthOrder3DTensor()
        end

        function createRandomTensor(obj)
            obj.createRandomTensor@FourthOrder3DTensor();
            obj.MakeMajorAndMinorSymmetrization();
        end

        function MakeMajorAndMinorSymmetrization(obj)
            Symmetrizer = FourthOrderSymmetrizer;
            obj.tensorValue = Symmetrizer.symmetrize(obj.tensorValue);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function s = addnaka(s)
%ADDNAKA Add the Nakagami distribution.

%   Copyright 1993-2004 The MathWorks, Inc.
%   $Revision: 1.1.6.7 $  $Date: 2003/12/11 03:50:49 $

j = length(s) + 1;
s(j).name = 'Nakagami';
s(j).code = 'nakagami';
s(j).pnames = {'mu' 'omega'};
s(j).pdescription = {'shape' 'scale'};
s(j).prequired = [false false];
s(j).fitfunc = @nakafit;
s(j).likefunc = @nakalike;
s(j).cdffunc = @nakacdf;
s(j).pdffunc = @nakapdf;
s(j).invfunc = @nakainv;
s(j).statfunc = @nakastat;
s(j).loginvfunc = [];
s(j).logcdffunc = [];
s(j).hasconfbounds = false;
s(j).censoring = true;
s(j).paramvec = true;
s(j).support = [0 Inf];
s(j).closedbound = [false false];
s(j).iscontinuous = true;
s(j).islocscale = false;
s(j).uselogpp = false;

% ==== Nakagami distribution functions ====

% these distribution functions do not yet handle arrays of parameters

function y = nakapdf(x, mu, omega)
%NAKAPDF Nakagami probability density function (pdf).
mu(mu <= 0) = NaN;
omega(omega <= 0) = NaN;

x(x<0) = 0;
% equivalent to y = 2.*x .* gampdf(x.^2, mu, omega./mu), but the version here
% puts all the x terms into gampdf, so Inf*0, etc. is handled there.
y = 2.*sqrt(omega./mu).*exp(gammaln(mu+.5) - gammaln(mu)) .* gampdf(x.^2, mu+.5, omega./mu);

function p = nakacdf(x, mu, omega)
%NAKACDF Nakagami cumulative distribution function (cdf).
mu(mu <= 0) = NaN;
omega(omega <= 0) = NaN;

x(x<0) = 0;
p = gamcdf(x.^2, mu, omega./mu);

function x = nakainv(p, mu, omega)
%NAKAINV Inverse of the Nakagami cumulative distribution function (cdf).
mu(mu <= 0) = NaN;
omega(omega <= 0) = NaN;

x = sqrt(gaminv(p,mu,omega./mu));

function r = nakarnd(mu, omega, varargin)
%NAKARND Random arrays from the Nakagami distribution.
mu(mu <= 0) = NaN;
omega(omega <= 0) = NaN;

[err, sizeOut] = statsizechk(2,mu,omega,varargin{:});
if err > 0
    error('stats:nakarnd:InconsistentSizes','Size information is inconsistent.');
end

r = sqrt(gamrnd(mu,omega./mu,sizeOut));

function [m,v] = nakastat(mu, omega)
%NAKASTAT Mean and variance for the Nakagami distribution.
mu(mu <= 0) = NaN;
omega(omega <= 0) = NaN;

gamratio = exp(gammaln(mu+.5) - gammaln(mu));
m = gamratio .* sqrt(omega./mu);
v = omega .* (1 - gamratio.^2 ./ mu);

function [nlogL,acov] = nakalike(params,data,cens,freq)
%NAKALIKE Negative log-likelihood for the Nakagami distribution.
if nargin < 4 || isempty(freq), freq = ones(size(data)); end
if nargin < 3 || isempty(cens), cens = zeros(size(data)); end

nlogL = naka_nloglf(params, data, cens, freq);
if nargout > 1
    acov = mlecov(params, data, 'nloglf',@naka_nloglf, 'cens',cens, 'freq',freq);
end

% ==== Nakagami fitting functions ====

function [phat,pci] = nakafit(x,alpha,cens,freq,opts)
%NAKAFIT Parameter estimates and confidence intervals for Nakagami data.

if nargin < 2 || isempty(alpha), alpha = .05; end
if nargin < 3 || isempty(cens), cens = zeros(size(x)); end
if nargin < 4 || isempty(freq), freq = ones(size(x)); end
if nargin < 5, opts = []; end

if any(x <= 0)
    error('stats:nakafit:BadData','The data in X must be positive');
end

phat = gamfit(x.^2,alpha,cens,freq,opts);
phat(2) = phat(1).*phat(2); % (a,b) -> (mu,omega)
if nargout > 1
    acov = mlecov(phat, x, 'nloglf',@naka_nloglf, 'cens',cens, 'freq',freq);
    probs = [alpha/2; 1-alpha/2];
    se = sqrt(diag(acov))';
    pci = norminv(repmat(probs,1,numel(phat)), [phat; phat], [se; se]);
    % CI on the log scale for omega?
end

function [nll,ngrad] = naka_nloglf(parms, x, cens, freq)
%NAKA_NLOGLF Objective function for Nakagami maximum likelihood.

% do all the calculations in terms of the gamma dist'n
a = parms(1);
b = parms(2)./parms(1); % (mu,omega) -> (a,b)
loggama = gammaln(a);
logb = log(b);

xsq = x.^2;
z = xsq ./ b;
logz = log(z);
L = (a-1).*logz - z - loggama - logb + log(2.*x);
ncen = sum(freq.*cens);
if ncen > 0
    cen = (cens == 1);
    zcen = z(cen);
    if nargout == 1
        Scen = gammainc(zcen,a,'upper');
    else
        [dScen,Scen] = dgammainc(zcen,a,'upper');
    end
    L(cen) = log(Scen);
end
nll = -sum(freq .* L);

if nargout > 1
    dL1 = logz - psi(a);
    dL2 = (z - a)./b;
    if ncen > 0
        dL1(cen) = dScen ./ Scen;
        dL2(cen) = exp(a.*logz(cen) - logb - zcen - loggama) ./ Scen;
    end
    ngrad = -[sum(freq .* dL1) sum(freq .* dL2)];

    % transform back to Nakagami parameters
    ngrad = ngrad * [1 0; -b./a 1./a]; % (a,b) -> (mu,omega)
end

author: zouchuhang

repo: LayoutNet


text:
function varargout = plotgrid(varargin)
% VL_PLOTGRID  Plot a 2-D grid
%  VL_PLOTGRID(X,Y) plots a grid with vertices (X,Y). X and Y are MxN
%  matrices, with one entry per vertex.
%
%  H=VL_PLOTGRID(...) returns the handle to the grid object.
%
%  See also: VL_HELP().
[varargout{1:nargout}] = vl_plotgrid(varargin{:});

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
function mpc = case14
%CASE14    Power flow data for IEEE 14 bus test case.
%   Please see CASEFORMAT for details on the case file format.
%   This data was converted from IEEE Common Data Format
%   (ieee14cdf.txt) on 15-Oct-2014 by cdf2matp, rev. 2393
%   See end of file for warnings generated during conversion.
%
%   Converted from IEEE CDF file from:
%       https://labs.ece.uw.edu/pstca/
%
%  08/19/93 UW ARCHIVE           100.0  1962 W IEEE 14 Bus Test Case

%   MATPOWER

%% MATPOWER Case Format : Version 2
mpc.version = '2';

%%-----  Power Flow Data  -----%%
%% system MVA base
mpc.baseMVA = 100;

%% bus data
%	bus_i	type	Pd	Qd	Gs	Bs	area	Vm	Va	baseKV	zone	Vmax	Vmin
mpc.bus = [
	1	3	0	0	0	0	1	1.06	0	0	1	1.06	0.94;
	2	2	21.7	12.7	0	0	1	1.045	-4.98	0	1	1.06	0.94;
	3	2	94.2	19	0	0	1	1.01	-12.72	0	1	1.06	0.94;
	4	1	47.8	-3.9	0	0	1	1.019	-10.33	0	1	1.06	0.94;
	5	1	7.6	1.6	0	0	1	1.02	-8.78	0	1	1.06	0.94;
	6	2	11.2	7.5	0	0	1	1.07	-14.22	0	1	1.06	0.94;
	7	1	0	0	0	0	1	1.062	-13.37	0	1	1.06	0.94;
	8	2	0	0	0	0	1	1.09	-13.36	0	1	1.06	0.94;
	9	1	29.5	16.6	0	19	1	1.056	-14.94	0	1	1.06	0.94;
	10	1	9	5.8	0	0	1	1.051	-15.1	0	1	1.06	0.94;
	11	1	3.5	1.8	0	0	1	1.057	-14.79	0	1	1.06	0.94;
	12	1	6.1	1.6	0	0	1	1.055	-15.07	0	1	1.06	0.94;
	13	1	13.5	5.8	0	0	1	1.05	-15.16	0	1	1.06	0.94;
	14	1	14.9	5	0	0	1	1.036	-16.04	0	1	1.06	0.94;
];

%% generator data
%	bus	Pg	Qg	Qmax	Qmin	Vg	mBase	status	Pmax	Pmin	Pc1	Pc2	Qc1min	Qc1max	Qc2min	Qc2max	ramp_agc	ramp_10	ramp_30	ramp_q	apf
mpc.gen = [
	1	232.4	-16.9	10	0	1.06	100	1	332.4	0	0	0	0	0	0	0	0	0	0	0	0;
	2	40	42.4	50	-40	1.045	100	1	140	0	0	0	0	0	0	0	0	0	0	0	0;
	3	0	23.4	40	0	1.01	100	1	100	0	0	0	0	0	0	0	0	0	0	0	0;
	6	0	12.2	24	-6	1.07	100	1	100	0	0	0	0	0	0	0	0	0	0	0	0;
	8	0	17.4	24	-6	1.09	100	1	100	0	0	0	0	0	0	0	0	0	0	0	0;
];

%% branch data
%	fbus	tbus	r	x	b	rateA	rateB	rateC	ratio	angle	status	angmin	angmax
mpc.branch = [
	1	2	0.01938	0.05917	0.0528	0	0	0	0	0	1	-360	360;
	1	5	0.05403	0.22304	0.0492	0	0	0	0	0	1	-360	360;
	2	3	0.04699	0.19797	0.0438	0	0	0	0	0	1	-360	360;
	2	4	0.05811	0.17632	0.034	0	0	0	0	0	1	-360	360;
	2	5	0.05695	0.17388	0.0346	0	0	0	0	0	1	-360	360;
	3	4	0.06701	0.17103	0.0128	0	0	0	0	0	1	-360	360;
	4	5	0.01335	0.04211	0	0	0	0	0	0	1	-360	360;
	4	7	0	0.20912	0	0	0	0	0.978	0	1	-360	360;
	4	9	0	0.55618	0	0	0	0	0.969	0	1	-360	360;
	5	6	0	0.25202	0	0	0	0	0.932	0	1	-360	360;
	6	11	0.09498	0.1989	0	0	0	0	0	0	1	-360	360;
	6	12	0.12291	0.25581	0	0	0	0	0	0	1	-360	360;
	6	13	0.06615	0.13027	0	0	0	0	0	0	1	-360	360;
	7	8	0	0.17615	0	0	0	0	0	0	1	-360	360;
	7	9	0	0.11001	0	0	0	0	0	0	1	-360	360;
	9	10	0.03181	0.0845	0	0	0	0	0	0	1	-360	360;
	9	14	0.12711	0.27038	0	0	0	0	0	0	1	-360	360;
	10	11	0.08205	0.19207	0	0	0	0	0	0	1	-360	360;
	12	13	0.22092	0.19988	0	0	0	0	0	0	1	-360	360;
	13	14	0.17093	0.34802	0	0	0	0	0	0	1	-360	360;
];

%%-----  OPF Data  -----%%
%% generator cost data
%	1	startup	shutdown	n	x1	y1	...	xn	yn
%	2	startup	shutdown	n	c(n-1)	...	c0
mpc.gencost = [
	2	0	0	3	0.0430292599	20	0;
	2	0	0	3	0.25	20	0;
	2	0	0	3	0.01	40	0;
	2	0	0	3	0.01	40	0;
	2	0	0	3	0.01	40	0;
];

%% bus names
mpc.bus_name = {
	'Bus 1     HV';
	'Bus 2     HV';
	'Bus 3     HV';
	'Bus 4     HV';
	'Bus 5     HV';
	'Bus 6     LV';
	'Bus 7     ZV';
	'Bus 8     TV';
	'Bus 9     LV';
	'Bus 10    LV';
	'Bus 11    LV';
	'Bus 12    LV';
	'Bus 13    LV';
	'Bus 14    LV';
};

% Warnings from cdf2matp conversion:
%
% ***** check the title format in the first line of the cdf file.
% ***** Qmax = Qmin at generator at bus    1 (Qmax set to Qmin + 10)
% ***** MVA limit of branch 1 - 2 not given, set to 0
% ***** MVA limit of branch 1 - 5 not given, set to 0
% ***** MVA limit of branch 2 - 3 not given, set to 0
% ***** MVA limit of branch 2 - 4 not given, set to 0
% ***** MVA limit of branch 2 - 5 not given, set to 0
% ***** MVA limit of branch 3 - 4 not given, set to 0
% ***** MVA limit of branch 4 - 5 not given, set to 0
% ***** MVA limit of branch 4 - 7 not given, set to 0
% ***** MVA limit of branch 4 - 9 not given, set to 0
% ***** MVA limit of branch 5 - 6 not given, set to 0
% ***** MVA limit of branch 6 - 11 not given, set to 0
% ***** MVA limit of branch 6 - 12 not given, set to 0
% ***** MVA limit of branch 6 - 13 not given, set to 0
% ***** MVA limit of branch 7 - 8 not given, set to 0
% ***** MVA limit of branch 7 - 9 not given, set to 0
% ***** MVA limit of branch 9 - 10 not given, set to 0
% ***** MVA limit of branch 9 - 14 not given, set to 0
% ***** MVA limit of branch 10 - 11 not given, set to 0
% ***** MVA limit of branch 12 - 13 not given, set to 0
% ***** MVA limit of branch 13 - 14 not given, set to 0

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function varargout = dblquad(varargin)
%DBLQUAD   Complete definite integral of SPHEREFUN.
%   I = DBLQUAD(F, a, b, c, d), returns the definite integral of a SPHEREFUN over
%   the region [a, b, c, d].
%
%   This function is a wrapper for quad2d.
%
% See also SPHEREFUN/QUAD2D, SPHEREFUN/INTEGRAL2, SPHEREFUN/SUM2.

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

[varargout{1:nargout}] = dblquad@separableApprox(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function ni = niftiSet(ni,param,val,varargin)
% Set data for various nifti data structures
%
%     ni = niftiSet(ni,param,val,varargin)
%
% Parameter values
%
%  'pixdim'    - Voxel dimensions in millimeters (though settable)
%  'dim'       - Number of pixels in each dimension (row,col, depth, time/orientation)
%  'data'      - The data
%  'voxelsize' - What is this?  Why isn't it pixdim?
%
%  'checkqto'
%
% Example:
%
%   ni = niftiCreate;
%   ni = niftiSet(ni,'voxel size',[0.375, 0.375, 2]);
%
%
% EKA/BW Vistasoft Team, 2015

%% Check the parameters
if notDefined('ni'), error('Nifti data structure variable required'); end
if notDefined('param'), error('Parameter field required.'); end
if ~exist('val','var'), error('Val required'); end

%% Squeeze spaces out, force lower case
param = mrvParamFormat(param);

%TODO: Add a nifti paramaterMapField
switch param
    case 'checkqto'

        sz = size(niftiGet(ni,'data'));
        dim = niftiGet(ni,'dim');

        if(any(dim(1:3)~=sz(1:3)))
            warning('[%s] NIFTI volume dim wrong- setting it to the actual data size.\n',mfilename);
            dim(1:3) = sz(1:3);
            ni = niftiSet(ni,'dim',dim);
        end

        if(ni.qform_code==0 && ni.sform_code~=0)
            warning('[%s] ni.qform_code is zero and sform_code ~=0. Setting ni.qto_* from ni.sto_*...\n',mfilename);
            %ni = niftiSetQto(ni, ni.sto_xyz);
            ni = niftiSet(ni,'qto',niftiGet(ni,'sto_xyz'));
        end

        dim = niftiGet(ni,'dim');
        qto_ijk = niftiGet(ni,'qto_ijk');
        origin = [qto_ijk(1:3,:)*[0 0 0 1]']';
        if(any(origin<2)||any(origin>dim(1:3)-2))
            [~,r,s,k] = affineDecompose(niftiGet(ni, 'qto_ijk'));
            t = ni.dim/2;
            warning_string = [sprintf('[%s] Qto matrix defines an origin very far away from the isocenter.\n',mfilename),...
                              sprintf('This implies that the Qto matrix may be bad - please check qto_ijk. An automatic fix will be attempted.\n'),...
                              sprintf('Origin to the image center is at [%2.3f,%2.3f,%2.3f] pix.\n',t(1),t(2),t(3))];
            warning(warning_string);
            %ni = niftiSetQto(ni, inv(affineBuild(t,r,s,k)));
            ni = niftiSet(ni,'qto',inv(affineBuild(t,r,s,k)));

        end

    case 'data'
        ni.data = val;

    case 'dim'
        ni.dim = val;

    case 'filepath'
        ni.fname = val;

    case 'freqdim'
        ni.freq_dim = val;

    case 'nifti'
        ni = val; %This means that we are passing in an entire Nifti!

    case 'phasedim'
        ni.phase_dim = val;

    case 'pixdim'
        ni.pixdim = val;

    case 'qfac'
        ni.qfac= val;

    case 'qform_code'
        ni.qform_code = val;

    case 'qoffset_x'
        ni.qoffset_x = val;

    case 'qoffset_y'
        ni.qoffset_y = val;

    case 'qoffset_z'
        ni.qoffset_z = val;

    case 'qto'
        xformXyz = val;
        q = matToQuat(xformXyz);
        ni = niftiSet(ni,'qform_code',2);
        ni = niftiSet(ni,'qto_xyz',xformXyz);
        ni = niftiSet(ni,'qto_ijk',inv(xformXyz));
        ni = niftiSet(ni,'quatern_b',q.quatern_b);
        ni = niftiSet(ni,'quatern_c ',q.quatern_c);
        ni = niftiSet(ni,'quatern_d',q.quatern_d);
        ni = niftiSet(ni,'qoffset_x',q.quatern_x);
        ni = niftiSet(ni,'qoffset_y',q.quatern_y);
        ni = niftiSet(ni,'qoffset_z',q.quatern_z);
        ni = niftiSet(ni,'qfac',q.qfac);

        if length(varargin) < 1
            setStoToo = '';
        else
            setStoToo = varargin{1};
        end

        if(~isempty(setStoToo)&&setStoToo)
            ni = niftiSet(ni,'sto_xyz',niftiGet(ni,'qto_xyz'));
            ni = niftiSet(ni,'sto_ijk',niftiGet(ni,'qto_ijk'));
        end

    case 'qto_ijk'
        ni.qto_ijk = val;

    case 'qto_xyz'
        ni.qto_xyz = val;

    case 'quatern_b'
        ni.quatern_b = val;

    case 'quatern_c'
        ni.quatern_c = val;

    case 'quatern_d'
        ni.quatern_d = val;

    case 'slicedim'
        ni.slice_dim = val;

    case 'sto_ijk'
        ni.sto_ijk = val;

    case 'sto_xyz'
        ni.sto_xyz = val;

    case 'voxelsize'
        ni.voxelSize = val;

    otherwise
        error('Unknown parameter %s\n',param);

end %switch

return


function [ save_totalNum,  visCells ] = deleteZeroPoints( weightFile, visFilesDir, compressResultsDir, saveReIndexFile, saveNewQFile )
%DELETEZEROPOINTS Summary of this function goes here
%   Detailed explanation goes here

    % get the union result
    disp('Prepare the next optimization');

    fileExt = '*.txt';
    compress_files = dir(fullfile(compressResultsDir,fileExt));
    all_ID = [];
    for i = 0:length(compress_files)-1

        fileCnt = num2str(i);
        fileName = [compressResultsDir, fileCnt, '.txt'];

        file_t = fopen(fileName);
        point_ID = fscanf(file_t, '%d');
        fclose(file_t);

        all_ID = union(all_ID, point_ID);

    end
    save_totalNum = length(all_ID);

    %%
    % get the one2one indexes relations and save
    % the map point index, in C++ from 0, in Matlab from 1
    % indexList id from ZERO
    indexList = [];
    indexList(:,1) = all_ID-1;
    for i= 1:length(all_ID)
        indexList(i,2) = i-1;  % map point re-index from zero
    end
    dlmwrite(saveReIndexFile, indexList, 'precision', '%d');

    % update the visMatrix and save

    fileExt = '*.txt';
    vis_files = dir(fullfile(visFilesDir,fileExt));
    for i=0:length(vis_files)-1
%         tic
%         disp(i);

        fileCnt = num2str(i);
        old_fileName = [visFilesDir, fileCnt, '.txt'];
        file_old = fopen(old_fileName);
        vis_point_ID_old = fscanf(file_old, '%d');

        % thought 3
        isInIndex = ismember(indexList(:,1), vis_point_ID_old);
        isSaved = find(isInIndex==1);
        newnew = indexList(:,2);
        vis_point_ID_new = newnew(isSaved);

        visCells{i+1} = vis_point_ID_new;

%         toc
    end

    % update the q-matrix-vector
    % weight is from ONE
    file_q = fopen(weightFile);
    weights = fscanf(file_q, '%d');
    new_weights = [];
    for i = 1:length(indexList)
        new_weights(i,:) = weights(indexList(i,1)+1);
    end
    dlmwrite(saveNewQFile, new_weights, 'precision', '%d');

    disp('Finished');

end

%         % visMatrix is from
%         vis_point_ID_new = [];
%         remain_cnt = 1;
        % thought 1
%         for j =1:length(vis_point_ID_old)
%             row_indexList = find(indexList(:,1) == (vis_point_ID_old(j)));
%             if size(row_indexList,1) == 1
%                 vis_point_ID_new(remain_cnt,:) = indexList(row_indexList,2);  % already from zero before this part
%                 remain_cnt = remain_cnt + 1;
%             else  % zero matrix
%             end
%         end
%         % thought 2
%         vis_point_ID_old_inter = intersect(vis_point_ID_old, indexList(:,1));
%         for j = 1:length(vis_point_ID_old_inter)
% %             rowLine = find(indexList(:,1) == vis_point_ID_old_inter(j));
%             vis_point_ID_old_inter(j,:) = indexList(find(indexList(:,1) == vis_point_ID_old_inter(j)), 2);
%         end
%


text:
function varargout = spm_changepath(Sf, oldp, newp, verbose)
% Recursively replace all occurences of a text pattern in a MATLAB variable.
% FORMAT S = spm_changepath(Sf, oldp, newp)
%
% Sf       - MATLAB variable to fix, or char array of MAT filenames,
%            or directory name (all found MAT files will be analysed)
% oldp     - old string to replace
% newp     - new string replacing oldp
%
% S        - updated MATLAB variable (only if Sf is one)
%
% If the pattern is found in a string, any occurence of an invalid file
% separator is replaced to match that of the current system.
%
% If MAT filenames are specified, they will be overwritten with the new
% version. A backup of the initial version is made with a ".old" suffix.
%__________________________________________________________________________
% Copyright (C) 2009-2016 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin
% $Id: spm_changepath.m 6948 2016-11-24 10:34:01Z guillaume $

%-Input arguments
%--------------------------------------------------------------------------
if ~nargin
    Sf = spm_select(Inf,'mat','Select MAT files to fix');
end

if nargin <= 1
    oldp = spm_input('Old pattern','+1','s');
end

if nargin <= 2
    newp = spm_input('New pattern','+1','s');
end

if nargin <= 3
    verbose = true;
end

%-Replace pattern in given MAT-files
%--------------------------------------------------------------------------
if ischar(Sf)
    if nargout
        error('Output argument only valid for MATLAB variable input');
    end
    if exist(Sf,'dir')
        Sf = spm_select('FPList',Sf,'^.*\.mat$'); % FPListRec for recursive
    end
    if isempty(Sf), Sf = {}; else Sf = cellstr(Sf); end
    for i=1:numel(Sf)
        f = Sf{i};
        try
            S = load(f);
        catch
            error('Cannot load "%s".',f);
        end
        tmp = changepath(S,oldp,newp,verbose);
        if ~isequalwithequalnans(tmp,S)
            if verbose, fprintf('=> Fixing %s\n',f); end
            [sts, msg] = movefile(f,[f '.old']);
            if ~sts, error(msg); end
            save(f ,'-struct','tmp', spm_get_defaults('mat.format'));
        end
    end
else
    varargout = { changepath(Sf,oldp,newp,verbose) };
end

%==========================================================================
function S = changepath(S,oldp,newp,verbose)

check = false;

switch class(S)
    case {'double','single','logical','int8','uint8','int16','uint16',...
            'int32','uint32','int64','uint64','function_handle'}

    case 'cell'
        for i=1:numel(S)
            S{i} = changepath(S{i},oldp,newp,verbose);
        end

    case 'struct'
        for i=1:numel(S)
            fn = fieldnames(S);
            for j=1:length(fn)
                S(i).(fn{j}) = changepath(S(i).(fn{j}),oldp,newp,verbose);
            end
        end

    case 'char'
        f = {'\' '/'};
        if ispc, f = fliplr(f); end
        tmp = cellstr(S);
        for i=1:numel(tmp)
            t = strrep(tmp{i},oldp,newp);
            if ~isequal(tmp{i},t)
                t = strrep(t,f{1},f{2});
                if check
                    if spm_existfile(t)
                        sts = ' (OK) ';
                    else
                        sts = ' (NOT FOUND) ';
                    end
                else
                    sts = '';
                end
                if verbose, fprintf('%s%s\n',t,sts); end
            end
            tmp{i} = t;
        end
        S = char(tmp);

    case 'nifti'
        for i=1:numel(S)
            S(i).dat = changepath(S(i).dat,oldp,newp,verbose);
        end

    case 'file_array'
        S.fname = changepath(S.fname,oldp,newp,verbose);

    case 'gifti'
        if isfield(S,'cdata') && isa(S.cdata,'file_array')
            fa = struct(S.cdata);
            fa.fname = changepath(fa.fname,oldp,newp,verbose);
            S.cdata = file_array(fa);
        end

    otherwise
        warning('Unknown class "%s".',class(S));
end

author: spm


text:
%DEMEV3	Demonstrate Bayesian regression for the RBF.
%
%	Description
%	The problem consists an input variable X which sampled from a
%	Gaussian distribution, and a target variable T generated by computing
%	SIN(2*PI*X) and adding Gaussian noise. An RBF network with linear
%	outputs is trained by minimizing a sum-of-squares error function with
%	isotropic Gaussian regularizer, using the scaled conjugate gradient
%	optimizer. The hyperparameters ALPHA and BETA are re-estimated using
%	the function EVIDENCE. A graph  is plotted of the original function,
%	the training data, the trained network function, and the error bars.
%
%	See also
%	DEMEV1, EVIDENCE, RBF, SCG, NETEVFWD
%

%	Copyright (c) Ian T Nabney (1996-2001)

clc;
disp('This demonstration illustrates the application of Bayesian')
disp('re-estimation to determine the hyperparameters in a simple regression')
disp('problem using an RBF netowk. It is based on a the fact that the')
disp('posterior distribution for the output weights of an RBF is Gaussian')
disp('and uses the evidence maximization framework of MacKay.')
disp(' ')
disp('First, we generate a synthetic data set consisting of a single input')
disp('variable x sampled from a Gaussian distribution, and a target variable')
disp('t obtained by evaluating sin(2*pi*x) and adding Gaussian noise.')
disp(' ')
disp('Press any key to see a plot of the data together with the sine function.')
pause;

% Generate the matrix of inputs x and targets t.

ndata = 16;			% Number of data points.
noise = 0.1;			% Standard deviation of noise distribution.
randn('state', 0);
rand('state', 0);
x = 0.25 + 0.07*randn(ndata, 1);
t = sin(2*pi*x) + noise*randn(size(x));

% Plot the data and the original sine function.
h = figure;
nplot = 200;
plotvals = linspace(0, 1, nplot)';
plot(x, t, 'ok')
xlabel('Input')
ylabel('Target')
hold on
axis([0 1 -1.5 1.5])
fplot('sin(2*pi*x)', [0 1], '-g')
legend('data', 'function');

disp(' ')
disp('Press any key to continue')
pause; clc;

disp('Next we create a two-layer MLP network having 3 hidden units and one')
disp('linear output. The model assumes Gaussian target noise governed by an')
disp('inverse variance hyperparmeter beta, and uses a simple Gaussian prior')
disp('distribution governed by an inverse variance hyperparameter alpha.')
disp(' ');
disp('The network weights and the hyperparameters are initialised and then')
disp('the output layer weights are optimized with the scaled conjugate gradient')
disp('algorithm using the SCG function, with the hyperparameters kept')
disp('fixed. After a maximum of 50 iterations, the hyperparameters are')
disp('re-estimated using the EVIDENCE function. The process of optimizing')
disp('the weights with fixed hyperparameters and then re-estimating the')
disp('hyperparameters is repeated for a total of 3 cycles.')
disp(' ')
disp('Press any key to train the network and determine the hyperparameters.')
pause;

% Set up network parameters.
nin = 1;		% Number of inputs.
nhidden = 3;		% Number of hidden units.
nout = 1;		% Number of outputs.
alpha = 0.01;		% Initial prior hyperparameter.
beta_init = 50.0;	% Initial noise hyperparameter.

% Create and initialize network weight vector.
net = rbf(nin, nhidden, nout, 'tps', 'linear', alpha, beta_init);
[net.mask, prior] = rbfprior('tps', nin, nhidden, nout, alpha, alpha);
net = netinit(net, prior);

options = foptions;
options(14) = 5;  % At most 5 EM iterations for basis functions
options(1) = -1;  % Turn off all messages
net = rbfsetbf(net, options, x);  % Initialise the basis functions

% Now train the network
nouter = 5;
ninner = 2;
options = foptions;
options(1) = 1;
options(2) = 1.0e-5;		% Absolute precision for weights.
options(3) = 1.0e-5;		% Precision for objective function.
options(14) = 50;		% Number of training cycles in inner loop.

% Train using scaled conjugate gradients, re-estimating alpha and beta.
for k = 1:nouter
  net = netopt(net, options, x, t, 'scg');
  [net, gamma] = evidence(net, x, t, ninner);
  fprintf(1, '\nRe-estimation cycle %d:\n', k);
  fprintf(1, '  alpha =  %8.5f\n', net.alpha);
  fprintf(1, '  beta  =  %8.5f\n', net.beta);
  fprintf(1, '  gamma =  %8.5f\n\n', gamma);
  disp(' ')
  disp('Press any key to continue.')
  pause;
end

fprintf(1, 'true beta: %f\n', 1/(noise*noise));

disp(' ')
disp('Network training and hyperparameter re-estimation are now complete.')
disp('Compare the final value for the hyperparameter beta with the true')
disp('value.')
disp(' ')
disp('Notice that the final error value is close to the number of data')
disp(['points (', num2str(ndata),') divided by two.'])
disp(' ')
disp('Press any key to continue.')
pause; clc;
disp('We can now plot the function represented by the trained network. This')
disp('corresponds to the mean of the predictive distribution. We can also')
disp('plot ''error bars'' representing one standard deviation of the')
disp('predictive distribution around the mean.')
disp(' ')
disp('Press any key to add the network function and error bars to the plot.')
pause;

% Evaluate error bars.
[y, sig2] = netevfwd(netpak(net), net, x, t, plotvals);
sig = sqrt(sig2);

% Plot the data, the original function, and the trained network function.
[y, z] = rbffwd(net, plotvals);
figure(h); hold on;
plot(plotvals, y, '-r')
xlabel('Input')
ylabel('Target')
plot(plotvals, y + sig, '-b');
plot(plotvals, y - sig, '-b');
legend('data', 'function', 'network', 'error bars');

disp(' ')
disp('Notice how the confidence interval spanned by the ''error bars'' is')
disp('smaller in the region of input space where the data density is high,')
disp('and becomes larger in regions away from the data.')
disp(' ')
disp('Press any key to end.')
pause; clc; close(h);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
%% DeglitchStack
% removes images from stack
% S = complete stack struct
% badEnergy=

function Snew=DeglitchStack(S,badEnergy)

numbad=length(badEnergy);
resolution=0.1;            % Half of best energy resolution in scan definition

Snew=S;

for counter=1:numbad            % Find bad energy points

    rmindex=find(Snew.eVenergy < badEnergy(counter)+resolution & ...
        Snew.eVenergy > badEnergy(counter)-resolution);

     if isempty(rmindex)==0

        Snew.eVenergy(rmindex)=[];
        Snew.spectr(:,:,rmindex)=[];

    else

        displayeV=num2str(badEnergy(counter));
        display(strcat('Error: no energy match for input = ',displayeV,' eV! Image not removed!!!'))
        clear displayeV;

    end

    clear rmindex

end

return

author: Sable

repo: mcbench-benchmarks


text:
function centroids = kMeansInitCentroids(X, K)
%KMEANSINITCENTROIDS This function initializes K centroids that are to be
%used in K-Means on the dataset X
%   centroids = KMEANSINITCENTROIDS(X, K) returns K initial centroids to be
%   used with the K-Means on the dataset X
%

% You should return this values correctly
centroids = zeros(K, size(X, 2));

% ====================== YOUR CODE HERE ======================
% Instructions: You should set centroids to randomly chosen examples from
%               the dataset X
%

for i = 1:K
    centroids(i, :) = X(randperm(length(X) ,1), :);
end

% =============================================================

end

repo: machine-learning


text:
function toms443_test ( )

%*****************************************************************************80
%
%% TOMS443_TEST tests the TOMS443 library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 June 2014
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TOMS443_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the TOMS443 library.\n' );

  toms443_test01 ( );
  toms443_test02 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TOMS433_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [x,y] = sample_lds(F, H, Q, R, init_state, T, models, G, u)
% SAMPLE_LDS Simulate a run of a (switching) stochastic linear dynamical system.
% [x,y] = switching_lds_draw(F, H, Q, R, init_state, models, G, u)
%
%   x(t+1) = F*x(t) + G*u(t) + w(t),  w ~ N(0, Q),  x(0) = init_state
%   y(t) =   H*x(t) + v(t),  v ~ N(0, R)
%
% Input:
% F(:,:,i) - the transition matrix for the i'th model
% H(:,:,i) - the observation matrix for the i'th model
% Q(:,:,i) - the transition covariance for the i'th model
% R(:,:,i) - the observation covariance for the i'th model
% init_state(:,i) - the initial mean for the i'th model
% T - the num. time steps to run for
%
% Optional inputs:
% models(t) - which model to use at time t. Default = ones(1,T)
% G(:,:,i) - the input matrix for the i'th model. Default = 0.
% u(:,t)   - the input vector at time t. Default = zeros(1,T)
%
% Output:
% x(:,t)    - the hidden state vector at time t.
% y(:,t)    - the observation vector at time t.

if ~iscell(F)
  F = num2cell(F, [1 2]);
  H = num2cell(H, [1 2]);
  Q = num2cell(Q, [1 2]);
  R = num2cell(R, [1 2]);
end

M = length(F);
%T = length(models);

if nargin < 7,
  models = ones(1,T);
end
if nargin < 8,
  G = num2cell(repmat(0, [1 1 M]));
  u = zeros(1,T);
end

[os ss] = size(H{1});
state_noise_samples = cell(1,M);
obs_noise_samples = cell(1,M);
for i=1:M
  state_noise_samples{i} = sample_gaussian(zeros(length(Q{i}),1), Q{i}, T)';
  obs_noise_samples{i} = sample_gaussian(zeros(length(R{i}),1), R{i}, T)';
end

x = zeros(ss, T);
y = zeros(os, T);

m = models(1);
x(:,1) = init_state(:,m);
y(:,1) = H{m}*x(:,1) + obs_noise_samples{m}(:,1);

for t=2:T
  m = models(t);
  x(:,t) = F{m}*x(:,t-1) + G{m}*u(:,t-1) + state_noise_samples{m}(:,t);
  y(:,t) = H{m}*x(:,t)  + obs_noise_samples{m}(:,t);
end


text:
function sparse_interp_nd_test ( )

%*****************************************************************************80
%
%% SPARSE_INTERP_ND_TEST tests the SPARSE_INTERP_ND library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 October 2012
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../r8lib' )

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPARSE_INTERP_ND_TEST\n' );
  fprintf ( 1, '  MATLAB version.\n' );
  fprintf ( 1, '  Test the SPARSE_INTERP_ND library.\n' );
  fprintf ( 1, '  The R8LIB library is also required.\n' );

  m = 1;
  sparse_max = 9;
  sparse_interp_nd_test01 ( m, sparse_max );

  m = 2;
  sparse_max = 9;
  sparse_interp_nd_test01 ( m, sparse_max );

  m = 3;
  sparse_max = 9;
  sparse_interp_nd_test01 ( m, sparse_max );

  m = 4;
  sparse_max = 7;
  sparse_interp_nd_test01 ( m, sparse_max );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPARSE_INTERP_ND_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  rmpath ( '../r8lib' )

  return
end


text:
function c4_uniform_01_test ( )

%*****************************************************************************80
%
%% C4_UNIFORM_01_TEST tests C4_UNIFORM_01.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 June 2006
%
%  Author:
%
%    John Burkardt
%
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'C4_UNIFORM_01_TEST\n' );
  fprintf ( 1, '  C4_UNIFORM_01 computes pseudorandom complex values\n' );
  fprintf ( 1, '  in the unit circle.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The initial seed is %d\n', seed );
  fprintf ( 1, '\n' );

  for i = 1 : 10
    [ x, seed ] = c4_uniform_01 ( seed );
    fprintf ( 1, '  %6d  ( %f, %f )\n', i, real ( x ), imag ( x ) );
  end

  return
end


text:
function pass = test_emptyObjects( pref )
% For empty diskfun objects, does each command deal with them
% appropriately?

% Check things work for empty diskfuns.
f = diskfun();
try
    f + f;
    2*f;
    f*2;
    f.^2;
    2.^f;
    f.^f;
    sqrt(f);
    sum(f);
    integral2(f);
    norm(f);
    squeeze(f);
    diff(f);
    cos(f);
    sin(f);
    sinh(f);
    f.^f + f;
    diag(f);
    trace(f);
    mean(f);
    minandmax2(f);
    median(f);
    flipud(f);
    flipdim(f,1);
    pass = 1;
catch
    pass = 0;
end
end


text:
function [ywarp,dwarp_dt,dwarp_dtheta,d2warp_dthetadt] = outwarp_negscaledpow(hyp,y,invflag)
%GPLITE_NOISEFUN Noise function for lite Gaussian Process regression.
%   SN2 = GPLITE_NOISEFUN(HYP,X,NOISEFUN) computes the GP noise function
%   NOISEFUN, that is the variance of observation noise evaluated at test
%   points X. HYP is a single column vector of noise function
%   hyperparameters. NOISEFUN is a numeric array whose elements specify
%   features of the noise function, as follows:
%
%   See also GPLITE_COVFUN, GPLITE_MEANFUN.

if nargin < 2; y = []; end
if nargin < 3 || isempty(invflag); invflag = false; else; invflag = true; end

if invflag && nargout > 1
    error('outwarp_fun:InverseOnly', ...
        ['When calling for the inverse output warping function, only one function output is expected.']);
end

%--------------------------------------------------------------------------
% CUSTOM: Number of hyperparameters
Noutwarp = 3;       % # hyperparameters of the output warping function
%--------------------------------------------------------------------------

N = size(y,1);      % Number of training points

% Return number of output warping function hyperparameters and additional info
if ischar(hyp)
    ywarp = Noutwarp;
    if nargout > 1

        % Initialize bounds for all hyperparameters
        outwarp_info.LB = -Inf(1,Noutwarp);
        outwarp_info.UB = Inf(1,Noutwarp);
        outwarp_info.PLB = -Inf(1,Noutwarp);
        outwarp_info.PUB = Inf(1,Noutwarp);
        outwarp_info.x0 = NaN(1,Noutwarp);

        %------------------------------------------------------------------
        % CUSTOM: Initialize hyperparameter bounds and other details

        % Threshold parameter
        outwarp_info.LB(1) = min(y);
        outwarp_info.UB(1) = max(y);
        outwarp_info.PLB(1) = min(y);
        outwarp_info.PUB(1) = max(y);
        outwarp_info.x0(1) = NaN;

        % Scaling parameter a (log space)
        outwarp_info.LB(2) = -Inf;
        outwarp_info.UB(2) = Inf;
        outwarp_info.PLB(2) = -2;
        outwarp_info.PUB(2) = 2;
        outwarp_info.x0(2) = 0;

        % Power exponent k (log space)
        outwarp_info.LB(3) = -Inf;
        outwarp_info.UB(3) = Inf;
        outwarp_info.PLB(3) = -3;
        outwarp_info.PUB(3) = 3;
        outwarp_info.x0(3) = 0;

        %------------------------------------------------------------------

        % Assign handle of current output warping function
        outwarp_info.outwarpfun = str2func(mfilename);

        % Plausible starting point
        idx_nan = isnan(outwarp_info.x0);
        outwarp_info.x0(idx_nan) = 0.5*(outwarp_info.PLB(idx_nan) + outwarp_info.PUB(idx_nan));

        dwarp_dt = outwarp_info;

    end

    return;
end

[Nhyp,Ns] = size(hyp);      % Hyperparameters and samples

if Nhyp ~= Noutwarp
    error('outwarp_fun:WrongLikHyp', ...
        ['Expected ' num2str(Noutwarp) ' output warping function hyperparameters, ' num2str(Nhyp) ' passed instead.']);
end
if Ns > 1
    error('outwarp_fun:nosampling', ...
        'Output warping function output is available only for one-sample hyperparameter inputs.');
end

%--------------------------------------------------------------------------
% CUSTOM: Compute output warping function and gradients

% Read hyperparameters
y0 = hyp(1);
a = exp(hyp(2));
k = exp(hyp(3));

% Compute output warping or inverse warping
ywarp = y;
idx = y < y0;
if invflag      % Inverse output warping
    ywarp(idx) = y0 - ((y0 - y(idx)).^(1/k))/a;
else            % Direct output warping
    adelta = a*(y0 - y(idx));
    adeltak = adelta.^k;
    ywarp(idx) = y0 - adeltak;
end

if nargout > 1
    % First-order derivative of output warping function in output space
    dwarp_dt = ones(size(y));
    adeltakm1 = adelta.^(k-1);

    dwarp_dt(idx) = a*k*adeltakm1;

    if nargout > 2
        % Gradient of output warping function wrt hyperparameters
        dwarp_dtheta = zeros(N,Noutwarp);

        dwarp_dtheta(idx,1) = 1 - a*k*adeltakm1;            % y0
        dwarp_dtheta(idx,2) = -k*adeltak;                   % log(a)
        dwarp_dtheta(idx,3) = -k*adeltak.*log(adelta);      % log(k)

        if nargout > 3
            % Gradient of derivative of output warping function
            d2warp_dthetadt = zeros(N,Noutwarp);

            d2warp_dthetadt(idx,1) = a^2*k*(k-1)*adelta.^(k-2);                % y0
            d2warp_dthetadt(idx,2) = a*k^2*adeltakm1;                           % log(a)
            d2warp_dthetadt(idx,3) = a*k*adeltakm1 + a*k^2*adeltakm1.*log(adelta);  % log(k)

        end

    end
end

end

author: acerbilab

repo: vbmc

save_path: github-repos/MATLAB/acerbilab-vbmc


text:
function inout_flow_png ( )

%*****************************************************************************80
%
%% INOUT_FLOW_PNG creates PNG image files from flow data.
%
%  Discussion:
%
%    This MATLAB script file reads the In/Out flow data:
%
%      geometry (XY values at nodes, assumed to be in 'xy.txt')
%      flow (UV values at nodes, in a sequence of files starting with 'up001.txt')
%
%    and plots the velocity vectors (U,V)(X,Y), and saves each plot as a
%    PNG file, presumably so the PNG files can be gathered into an
%    animation.
%
%    The file plots either the velocity vector field, or the velocity
%    direction field, depending on the value of the internal logical
%    parameter "normalized".
%
%    The MATLAB routine QUIVER will internally scale the vectors, but this
%    can be adjusted by setting SCALE to a value that is not 1.
%
%    For the unnormalized case, the velocity field is computed to a fixed
%    scale determined by finding the largest velocity vector over time.
%
%    This routine requires some auxiliary routines in order to "increment"
%    the name of the current velocity data file to get the name of the next
%    one, with the main such routine being called FILE_NAME_INC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 July 2003
%
%  Author:
%
%    John Burkardt
%
  FALSE = 0;
  TRUE = 1;
  scale = 1.0;
  normalized = TRUE;

  nframes = 500;
%
%  Get the XY coordinates of the nodes.
%
  load xy.txt;

  x = xy(:,1);
  y = xy(:,2);
%
%  Thin the data.
%
  thin_factor = 2;
  thin_dex = thin_index ( x, y, thin_factor );
  thin_num = length ( thin_dex );
  x = x(thin_dex);
  y = y(thin_dex);
%
%  For unnormalized plots, you need to make sure that a fixed scale is preserved from
%  step to step.  The only way I can see to do this requires that I determine the
%  maximum velocity over all time, and then append one extract node and velocity
%  to the data structure.
%
  if ( normalized == FALSE )

    upfile = 'up000.txt';
    vnorm_max = 0.0;

    for ( i = 1 : nframes )
      fprintf ( 1, '  Checking file %d.\n', i );
      upfile = file_name_inc ( upfile );
      uv = load ( upfile );
      u = uv(thin_dex,1);
      v = uv(thin_dex,2);
      norm = sqrt ( u.^2 + v.^2 );
      vnorm_max = max ( vnorm_max, max ( norm ) );
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Maximum visible velocity magnitude is %f\n', vnorm_max );

  end
%
%  Set the coordinates of the boundary lines, and an extra set of
%  lines that will be invisible, but which make some space for the
%  file name to be displayed within the plot area.
%
  bx1 = [ 0.00, 1.00, 1.00, 0.99, 0.99, 0.00, 0.00 ];
  by1 = [ 0.00, 0.00, 0.80, 0.80, 0.01, 0.01, 0.00 ];

  bx2 = [ 0.00, 0.00, 1.00, 1.00, 0.01, 0.01, 0.00 ];
  by2 = [ 0.20, 1.00, 1.00, 0.99, 0.99, 0.20, 0.20 ];

  bx3 = [ -0.10,  1.10, 1.10, -0.10, -0.10 ];
  by3 = [ -0.10, -0.10, 1.10,  1.10, -0.10 ];
%
%  Set the name of the 0th (nonexistent) velocity file.
%  All the velocity files with have this format, with the
%  numeric part of the name incremented to get the next one.
%  In particular, the first velocity file is called UP001.TXT.
%
  upfile = 'up000.txt';
  pngfile = 'up000.png';

  for ( i = 1 : nframes )

    fprintf ( 1, '  Converting file %d.\n', i );

    upfile = file_name_inc ( upfile );
    pngfile = file_name_inc ( pngfile );
    uv = load ( upfile );
    u = uv(thin_dex,1);
    v = uv(thin_dex,2);

    if ( normalized == TRUE )
      norm = sqrt ( u.^2 + v.^2 );
      nonzero = find ( norm ~= 0.0 );
      u(nonzero) = u(nonzero) ./ norm(nonzero);
      v(nonzero) = v(nonzero) ./ norm(nonzero);
      vnorm_max = 1.0;
    end

    x(thin_num+1) = 0.00;
    y(thin_num+1) = 1.05;
    u(thin_num+1) = vnorm_max;
    v(thin_num+1) = 0.0;

    quiver ( x,  y, u, v, scale )
%
%  Draw the boundary lines, and the invisible bounding line.
%
    line ( bx1, by1, 'color', 'r' )
    line ( bx2, by2, 'color', 'r' )
    line ( bx3, by3, 'color', 'w' )

    axis equal
    if ( normalized == TRUE )
      title ( 'In/Out Direction Field' )
    else
      title ( 'In/Out Flow Field' )
    end
    text ( 0.420, 1.03, upfile )
%
%  Here's where we take a snapshot of the current image, and save it to
%  a PNG file.  Doing this was surprisingly complicated, and the help
%  system was not very clear.  Surely there's a more straightforward way!
%
    F = getframe;
    [X,map] = frame2im ( F );
    imwrite ( X, pngfile, 'PNG' );

  end

  return
end
function file_name = file_name_inc ( file_name )

%*****************************************************************************80
%
%% FILE_NAME_INC generates the next filename in a series.
%
%  Discussion:
%
%    It is assumed that the digits in the name, whether scattered or
%    connected, represent a number that is to be increased by 1 on
%    each call.  If this number is all 9's on input, the output number
%    is all 0's.  Non-numeric letters of the name are unaffected..
%
%    If the name is empty, then the routine stops.
%
%    If the name contains no digits, the empty string is returned.
%
%  Example:
%
%      Input            Output
%      -----            ------
%      'a7to11.txt'     'a7to12.txt'  (typical case.  Last digit incremented)
%      'a7to99.txt'     'a8to00.txt'  (last digit incremented, with carry.)
%      'a9to99.txt'     'a0to00.txt'  (wrap around)
%      'cat.txt'        ' '           (no digits in input name.)
%      ' '              STOP!         (error.)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 September 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character FILE_NAME(*), the string to be incremented.
%
%    Output, character FILE_NAME_NEW(*), the incremented string.
%
  lens = length ( file_name );

  if ( lens <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_NAME_INC - Fatal error!\n' );
    fprintf ( 1, '  The input filename is empty.\n' );
    error ( 'FILE_NAME_INC - Fatal error!' );
  end

  change = 0;

  for i = lens : -1 : 1

    c = file_name(i);

    if ( '0' <= c & c <= '8' )

      change = change + 1;

      c = c + 1;

      file_name(i) = c;

      return

    elseif ( c == '9' )

      change = change + 1;

      c = '0';

      file_name(i) = c;

    end

  end

  if ( change == 0 )
    file_name = ' ';
  end

  return
end
function thin_dex = thin_index ( x, y, thin_factor )

%*****************************************************************************80
%
%%  THIN_INDEX determines thinning indices for a X, Y data.
%
%  Discussion:
%
%    A set of X, Y data is given, that is presumably, not too far off
%    from being on a rectangular grid.
%
%    The input value of THIN_FACTOR indicates by how much the data should
%    be thinned.
%
%    The X and Y ranges are computed, and only those data points are
%    retained for which both X and Y lie in an appropriate subrange.
%
%    For instance, a THIN_FACTOR of 2 would essentially save data
%    that lay in the black squares of a checkerboard.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 June
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X(NODE_NUM), Y(NODE_NUM), the X and Y coordinates
%    of the nodes.
%
%    Input, integer THIN_FACTOR, the thinning factor.
%
%    Output, integer THIN_DEX(NODE_NUM), contains in (1:THIN_NUM) the
%    indices into X and Y of the vectors to be retained after thinning.
%
  TRUE = 1;
  FALSE = 0;

  node_num = length ( x );

  x_unique_num = 0;

  for ( i = 1 : node_num )

    unique = TRUE;

    for ( j = 1 : x_unique_num )
      if ( x(i) == x_unique(j) )
        unique = FALSE;
        break;
      end
    end

    if ( unique )
      x_unique_num = x_unique_num + 1;
      x_unique(x_unique_num) = x(i);
    end

  end

  sort ( x_unique );

  y_unique_num = 0;

  for ( i = 1 : node_num )

    unique = TRUE;

    for ( j = 1 : y_unique_num )
      if ( y(i) == y_unique(j) )
        unique = FALSE;
        break;
      end
    end

    if ( unique )
      y_unique_num = y_unique_num + 1;
      y_unique(y_unique_num) = y(i);
    end

  end

  sort ( y_unique );

  thin_num = 0;

  for ( i = 1 : node_num )

    for ( j = 1 : x_unique_num-1 )
      if ( x_unique(j) <= x(i) & x(i) <= x_unique(j+1) )
        x_bin = j;
        break;
      end
    end

    for ( j = 1 : y_unique_num-1 )
      if ( y_unique(j) <= y(i) & y(i) <= y_unique(j+1) )
        y_bin = j;
        break;
      end
    end

    if ( mod ( y_bin, thin_factor ) == thin_factor / 2 & ...
         mod ( x_bin, thin_factor ) == thin_factor / 2 )

      thin_num = thin_num + 1;
      thin_dex(thin_num) = i;

    end

  end

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end


text:
function [ fxx, fxy, fyy ] = f02_f2 ( n, x, y )

%*****************************************************************************80
%
%% F02_F2 returns second derivatives of function 2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 January 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real X(N,1), Y(N,1), the evalution points.
%
%    Output, real FXX(N,1), FXY(N,1), FYY(N,1), second derivatives.
%
  t1(1:n,1) = 18.0 * ( y(1:n,1) - x(1:n,1) );

  fxx(1:n,1) = 18.0 * tanh ( 0.5 * t1(1:n,1) ) ...
    .* ( tanh ( 9.0 * ( y(1:n,1) - x(1:n,1) ) ) + 1.0 ) / 9.0
  fxy(1:n,1) = - fxx(1:n,1);
  fyy(1:n,1) = fxx(1:n,1);

  return
end


text:
function fx = p01_fun ( n, x )

%*****************************************************************************80
%
%% P01_FUN evaluates the integrand for problem 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 January 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Gwynne Evans,
%    Practical Numerical Integration,
%    Wiley, 1993,
%    ISBN: 047193898X,
%    LC: QA299.3E93.
%
%  Parameters:
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real X(2,N), the evaluation points.
%
%    Output, real FX(N,1), the integrand values.
%
  fx(1:n,1) = 1.0 ./ ( 1.0 - x(1,1:n) .* x(2,1:n) );

  return
end


text:
function [vertices, faces] = freesurfer_read_surf(fname)

% freesurfer_read_surf - FreeSurfer I/O function to read a surface file
%
% [vertices, faces] = freesurfer_read_surf(fname)
%
% Reads the vertex coordinates (mm) and face lists from a surface file.
%
% Surface files are stored as either triangulations or quadrangulations.
% That is, for a triangulation, each face is defined by 3 vertices.  For a
% quadrangulation, each face is defined by 4 vertices.  The rows of 'faces'
% contain indices into the rows of 'vertices', the latter holds the XYZ
% coordinates of each vertex.
%
% The freesurfer faces index the vertices in counter-clockwise order (when
% viewed from the outside of the surface).  This is consistent with a
% right-hand rule.  If we have vertices
%
% C           B
%
%
%       A
%
% Then we can calculate an edge vector from A to B (ie, AB = B - A) and
% another edge vector from A to C (ie, AC = C - A).  If you form a "gun"
% with your thumb and forefinger of the right hand, then align your thumb
% with the AB vector and your forefinger with the AC vector, your palm is
% facing out of the screen and extending your middle finger in the
% orthogonal direction to the plane of the screen will give the outward
% surface normal of the triangle ABC.  (If you lookup "triangle" on
% Wolfram's mathworld, you can see that AB is referred to as c and AC is
% referred to as b.)
%
% However, if this surface is read into matlab, it will give INWARD surface
% normals in the matlab patch command.  For some reason, matlab is not
% following the right hand rule.  To get OUTWARD normals with the matlab
% patch command, use faces(:,[1 3 2]) (see below).
%
% The vertex coordinates are in mm.  The FreeSurfer coordinate
% system for surfaces is quite simple, but relating to their MRI
% cor-??? files is too confusing to explain here; see the FreeSurfer
% homepage or google the documentation by Graham Wideman.  For the
% surfaces, at least, the origin is somewhere in the center of the
% head, and the vertex XYZ coordinates are oriented such that +X is
% right, +Y is anterior and +Z is superior (this is the
% FreeSurfer RAS coordinate system).
%
% Note that reading the faces of a quad file can take a long
% time due to their compact storage format.  In this case, the return of
% vertices can be faster if the face output variable is not specified; in
% this case, the faces are not read.
%
% Try this to visualize the surface:
% Hp = patch('vertices',vertices,'faces',faces(:,[1 3 2]),...
%       'facecolor',[.5 .5 .5],'edgecolor','none')
% camlight('headlight','infinite')
% vertnormals = get(Hp,'vertexnormals');
%
% See also freesurfer_write_surf, freesurfer_read_curv,
%          freesurfer_read_wfile
%

% $Revision: 1.2 $ $Date: 2011/02/07 21:47:40 $

% Copyright (C) 2000  Darren L. Weber

% History:  08/2000, Darren.Weber_at_radiology.ucsf.edu
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ver = '$Revision: 1.2 $ $Date: 2011/02/07 21:47:40 $';
fprintf('FREESURFER_READ_SURF [v %s]\n',ver(11:15));

if(nargin < 1)
    help freesurfer_read_surf;
    return;
end

%QUAD_FILE_MAGIC_NUMBER =  (-1 & 0x00ffffff) ;
%NEW_QUAD_FILE_MAGIC_NUMBER =  (-3 & 0x00ffffff) ;

TRIANGLE_FILE_MAGIC_NUMBER  =  16777214;
QUAD_FILE_MAGIC_NUMBER      =  16777215;

% open it as a big-endian file
fid = fopen(fname, 'rb', 'b');
if (fid < 0),
    str = sprintf('could not open surface file %s.', fname);
    error(str);
end

fprintf('...reading surface file: %s\n', fname);
tic;

magic = freesurfer_fread3(fid);

if (magic == QUAD_FILE_MAGIC_NUMBER),
    Nvertices = freesurfer_fread3(fid);
    Nfaces = freesurfer_fread3(fid);
    fprintf('...reading %d quad file vertices\n',Nvertices);
    vertices = fread(fid, Nvertices*3, 'int16') ./ 100 ;
    if (nargout > 1),
        fprintf('...reading %d quad file faces (please wait)\n',Nfaces);
        faces = zeros(Nfaces,4);
        for iface = 1:Nfaces,
            for n=1:4,
                faces(iface,n) = freesurfer_fread3(fid) ;
            end
            if(~rem(iface, 10000)), fprintf(' %7.0f',iface); end
            if(~rem(iface,100000)), fprintf('\n'); end
        end
    end
elseif (magic == TRIANGLE_FILE_MAGIC_NUMBER),
    fprintf('...reading triangle file\n');
    tline = fgets(fid); % read creation date text line
    tline = fgets(fid); % read info text line

    Nvertices = fread(fid, 1, 'int32'); % number of vertices
    Nfaces = fread(fid, 1, 'int32'); % number of faces

    % vertices are read in column format and reshaped below
    vertices = fread(fid, Nvertices*3, 'float32');

    % faces are read in column format and reshaped
    faces = fread(fid, Nfaces*3, 'int32');
    faces = reshape(faces, 3, Nfaces)';
else
    str = sprintf('unknown magic number in surface file %s.', fname);
    error(str);
end

vertices = reshape(vertices, 3, Nvertices)';
fclose(fid);

fprintf('...adding 1 to face indices for matlab compatibility.\n');
faces = faces + 1;

t=toc; fprintf('...done (%6.2f sec)\n\n',t);

return


text:
function x = L4Pinv(cf,y)
%L4PINV The inverse of the 4 parameters logistic equation.
% The Four Parameters Logistic Regression or 4PL nonlinear regression model
% is commonly used for curve-fitting analysis in bioassays or immunoassays
% such as ELISAs or dose-response curves.
%
% Syntax: x=L4Pinv(cf,y)
%
% Inputs:
%           cf is the object containing the 4 parameters of logistic
%           equation computed by L4P function. Alternatively, it can be a
%           1x4 array.
%
%           y is the array of the response that you want to iterpolate.
%
% Outputs:
%           x is the vector of interpolated data.
%
% Example:
%
% xs=[0 4.5 10.6 19.7 40 84 210]; ys=[0.0089 0.0419 0.0873 0.2599 0.7074 1.528 2.7739];
%
% Calling on MatLab the function: [cf G]=L4P(x,y);
%
% you will find the 5 parameters of this curve.
%
% Calling on MatLab the function L4Pinv(cf,1.782315);
%
%           Answer is:
% ans =
%
%  99.9982
%
% Alternatively, you can do:
%
% P=[0.0010 1.5153 108.0035 3.7841]; L4Pinv(P,1.782315);
%
% with the same result.
%
%           Created by Giuseppe Cardillo
%           giuseppe.cardillo-edta@poste.it
%
% See also L4P, L5P, L5Pinv, L3P, L3Pinv
%
% To cite this file, this would be an appropriate format:
% Cardillo G. (2012) Four parameters logistic regression - There and back again
%

%--------------------Input errors handling section-------------------------
if nargin < 2
    error('Almost two inputs are required')
end
if isobject(cf)
    p=coeffvalues(cf);
else
    p=cf;
    if ~isvector(p) || length(p)~=4 || ~all(isfinite(p))
        error('cf must be a fit object or a 1x4 vector of real and finite numbers')
    end
end

%-------------------------Interpolate--------------------------------------
x=NaN(size(y)); ok=isfinite(y);
x(ok)=p(3).*(((p(1)-p(4))./(y(ok)-p(4)))-1).^(1/p(2));


text:
function [J, grad] = costFunctionReg(theta, X, y, lambda)
%COSTFUNCTIONREG Compute cost and gradient for logistic regression with regularization
%   J = COSTFUNCTIONREG(theta, X, y, lambda) computes the cost of using
%   theta as the parameter for regularized logistic regression and the
%   gradient of the cost w.r.t. to the parameters.

% Initialize some useful values
m = length(y); % number of training examples

% You need to return the following variables correctly
J = 0;
grad = zeros(size(theta));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost of a particular choice of theta.
%               You should set J to the cost.
%               Compute the partial derivatives and set grad to the partial
%               derivatives of the cost w.r.t. each parameter in theta

hx = sigmoid(X * theta); %hypothesis, m * 1

J = 1 / m * sum(-y' * log(hx) - (1 - y)' * log(1 - hx)) + lambda / (2 * m) * theta(2:end)' * theta(2:end);

gradf = (1 / m) * (X(:, 1)' * (hx - y));
gradb = (1 / m) * (X(:, 2:end)' * (hx - y)) + lambda * theta(2:end) / m;

grad = [gradf;gradb];

% =============================================================

end

author: scruel

repo: Notes-ML-AndrewNg


text:
function varargout = pm_ff_unwrap(varargin)
% Performs phase-unwrapping along the shore of a flood-fill
% progressing from some seed in a low-variance area towards
% higher variance areas.
% FORMAT [pm,wm] = pm_ff_unwrap(pm,vm,wm,mask,thres)
%
% Input:
% pm     : 2 or 3D phasemap where some voxels have been unwrapped
%          and some not.
% vm     : Variance map, indicating the variance of the phase estimate.
% wm     : Wrap-map, where a non-zero value indicates corresponding
%          phase-value in pm has been unwrapped.
% mask   : Mask that indicates which voxels are worth
%          bothering with and which are not.
% thres  : Array of thresholds such that unwrapping proceeds by
%          serial dilates from some seed point in wm until all
%          connected voxels with vm<thres(1) have been unwrapped,
%          after which the same happens with vm<thres(2) etc etc
%
% Output:
% pm     : Same as pm but with unwrapping now completed as far as
%          indicated by wm on output.
% wm     : Indicates how far unwrapping has progressed.
%
% Note that the following are equivalent in terms of the final
% result
%
% [pm,wm] = pm_ff_unwrap(pm,vm,wm,mask,thres);
%
% for i=1:length(thres)
%   [pm,vm] = pm_ff_unwrap(pm,vm,wm,maskm,thres(i));
% end
%
% where the later would be slightly slower but would allow
% reporting of progress to some GUI.
%_________________________________________________________________
% Jesper Andersson 30/9-03
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Jesper Andersson
% $Id: pm_ff_unwrap.m 1317 2008-04-08 16:16:38Z chloe $

error('mex-function pm_ff_unwrap.c not compiled');


text:
function [regions, thresholds, maxSig] = imMultiOtsuThreshold(img, nClasses, varargin)
% Multilevel Thresholding using Otsu Method.
%
%   CLASSES = imMultiOtsuThreshold(IMG, NC)
%   Computes a segmentation of the input grayscale image IMG into NC
%   classes. The number of classes must be comprised between 2 and 5.
%
%   [CLASSES, THRESHOLDS] = imMultiOtsuThreshold(IMG, NC)
%   Also returns the list of threshold. The number of threshold values is
%   the number of classes minus one.
%
%   Example
%     % segment cameraman image into three classes
%     img = imread('cameraman.tif');
%     [classes, threshs] = imMultiOtsuThreshold(img, 3);
%     rgb = label2rgb(classes, 'jet', [1 0 0], 'shuffle');
%     figure; imshow(rgb);
%     % also displays histogram with threshold levels
%     figure; imHistogram(img);
%     hold on;
%     for i = 1:2
%        plot(threshs([i i]), [0 1800], 'color', 'r', 'linewidth', 2);
%     end
%
%   References
%   * Ping-Sung Liao, Tse-Sheng Chen, Pau-Choo Chung (2001). "A Fast
%   Algorithm for Multilevel Thresholding". Journal of Information Science
%   and Engineering, Vol. 17 No. 5, pp. 713-727.
%   https://jise.iis.sinica.edu.tw/JISESearch/pages/View/PaperView.jsf?keyId=86_1302#
%   * https://imagej.net/Multi_Otsu_Threshold
%
%   See also
%     imOtsuThreshold, imHistogram, imMaxEntropyThreshold
%

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% INRAE - BIA Research Unit - BIBS Platform (Nantes)
% Created: 2021-01-29,    using Matlab 9.8.0.1323502 (R2020a)
% Copyright 2021 INRAE.

%% Input arguments

if ~isa(img, 'uint8')
    error('Requires a grayscale uint8 image as first input');
end

if nClasses < 2 || nClasses > 5
    error('The number of classes must be comprised between 2 and 5.');
end

%% Initialisations

% compute histogram, and convert into probability density
[h, levels] = imHistogram(img, varargin{:});
h = h / sum(h);

% number of gray levels
nLevels = length(levels);

Puv = zeros(nLevels, nLevels);
Suv = zeros(nLevels, nLevels);
Huv = zeros(nLevels, nLevels);

% initialize diagonal terms
for i = 1:nLevels
    Puv(i, i) = h(i);
    Suv(i, i) = h(i) * i;
end

% initialize first rows
for i = 2:nLevels
    Puv(1, i) = Puv(1, i-1) + h(i);
    Suv(1, i) = Suv(1, i-1) + i*h(i);
end
% initialize the rest of the matrices
for u = 2:nLevels
    for v = u+1:nLevels
        Puv(u, v) = Puv(1, v) - Puv(1, u-1);
        Suv(u, v) = Suv(1, v) - Suv(1, u-1);
    end
end

% now calculate Huv
for u = 1:nLevels
    for v = u+1:nLevels
        if Puv(u,v) > 0
            Huv(u,v) = Suv(u,v) * Suv(u,v) / Puv(u,v);
        end
    end
end

%% Main processing

% create array for threshold values (first one is zeros, last one is inf)
thresholds = zeros(nClasses+1, 1);
thresholds(end) = inf;

% iterate over all possible combinations
maxSig = 0.0;
switch nClasses
    case 2
        % two classes -> only need to find second threshold value
        for i = 1:nLevels-1
            Sq = Huv(1,i) + Huv(i+1, end);
            if Sq >= maxSig
                thresholds(2) = i;
                maxSig = Sq;
            end
        end

    case 3
        % three classes
        for i = 1:nLevels-2
            for j = i+1:nLevels-1
                Sq = Huv(1,i) + Huv(i+1, j) + Huv(j+1, end);
                if Sq >= maxSig
                    thresholds(2) = i;
                    thresholds(3) = j;
                    maxSig = Sq;
                end
            end
        end

    case 4
        % four classes
        for i = 1:nLevels-3
            for j = i+1:nLevels-2
                for k = j+1:nLevels-1
                    Sq = Huv(1,i) + Huv(i+1, j) + Huv(j+1, k) + Huv(k+1, end);
                    if Sq >= maxSig
                        thresholds(2) = i;
                        thresholds(3) = j;
                        thresholds(4) = k;
                        maxSig = Sq;
                    end
                end
            end
        end

    case 5
        % five classes
        for i = 1:nLevels-3
            for j = i+1:nLevels-2
                for k = j+1:nLevels-1
                    for m = k+1:nLevels-1
                        Sq = Huv(1,i) + Huv(i+1, j) + Huv(j+1, k) + Huv(k+1, m) + Huv(m+1, end);
                        if Sq >= maxSig
                            thresholds(2) = i;
                            thresholds(3) = j;
                            thresholds(4) = k;
                            thresholds(5) = m;
                            maxSig = Sq;
                        end
                    end
                end
            end
        end

    otherwise
        error('Can not manage %d number of classes', nClasses);
end

%% Create region image

regions = zeros(size(img), 'uint8');
for i = 1:nClasses
    inds = img >= thresholds(i) & img < thresholds(i+1);
    regions(inds) = i;
end

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function pwl_interp_1d_test ( )

%*****************************************************************************80
%
%% PWL_INTERP_1D_TEST tests the PWL_INTERP_1D library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 October 2012
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../test_interp' );
  addpath ( '../r8lib' );

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PWL_INTERP_1D_TEST:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the PWL_INTERP_1D library.\n' );
  fprintf ( 1, '  The R8LIB library is needed.\n' );
  fprintf ( 1, '  The test needs the TEST_INTERP library.\n' );

  prob_num = p00_prob_num ( );
  for prob = 1 : prob_num
    pwl_interp_1d_test01 ( prob );
  end

  prob = 7;
  pwl_interp_1d_test02 ( prob );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PWL_INTERP_1D_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  rmpath ( '../r8lib' );
  rmpath ( '../test_interp' );

  return
end