text:
function [val,vec,info] = syev(jobz,uplo,a)
%SYEV   Computing eigen information.
%   [VAL,VEC,INFO] = SYEV(JOBZ,UPLO,A) computes eigenvalues and,
%	optionally, eigenvectors of a symmetric matrix A.
%
%	JOBZ allows the user to choose if eigenvectors are to be computed.
%	JOBZ
%          = 'N':  Compute eigenvalues only;
%          = 'V':  Compute eigenvalues and eigenvectors.
%
%	UPLO allows the user to choose which part of the matrix will be referenced.
%	UPLO
%          = 'U':  Upper triangle of A is stored;
%          = 'L':  Lower triangle of A is stored.
%
%   If only eigenvalues are desired, the QR algorithm is used.
%
%	No input parameter are optional. If you want an easy to use interface,
%	see SYEV_DRIVER.

    vec = [];
	val = [];
	info = [];

	if (nargin~=3),
		disp('Wrong number of input parameters.');
		return;
	end;

	    %Do all the boring checking
	if (~isnumeric(a)),
		disp('The matrix is not composed of numeric values.');
		return;
	elseif (isinteger(a)),
		a=double(a);
	end;

    if (size(a,1)~=size(a,2)),
        disp('The matrix must be square.');
        return;
	end;
	n=int32(size(a,1));
	lda=int32(n);

	if (isComplex(a)),
		if (isDouble(a))

			w=zeros(double(n),1);
			work=complex(zeros(1));
			lrwork=int32(max(1,3*n-2));
			rwork=zeros(double(lrwork),1);
			lwork=int32(-1);
			info=int32(-1);
			[jobz, uplo, n, vec, lda, val, work, lwork, rwork,info]=lapack_zheev(jobz, uplo, n, a, lda, w, work, lwork, rwork,info);
			lwork=int32(work(1));
			work=complex(zeros(double(lwork),1));
			[jobz, uplo, n, vec, lda, val, work, lwork, rwork,info]=lapack_zheev(jobz, uplo, n, a, lda, w, work, lwork, rwork,info);
		else
			w=single(zeros(double(n),1));
			work=single(complex(zeros(1,1)));
			work(1)=sqrt(-1)+1;
			lrwork=int32(single(max(1,3*n-2)));
			rwork=single(zeros(double(lrwork),1));
			lwork=int32(single(-1));
			info=int32(single(-1));
			[jobz, uplo, n, vec, lda, val, work, lwork, rwork,info]=lapack_cheev(jobz, uplo, n, a, lda, w, work, lwork, rwork,info);
			lwork=int32(work(1));
			work=single(complex(zeros(double(lwork),1)));
			[jobz, uplo, n, vec, lda, val, work, lwork, rwork,info]=lapack_cheev(jobz, uplo, n, a, lda, w, work, lwork, rwork,info);

		end;
	else
		if (isDouble(a)),
			w=zeros(double(n),1);
			work=zeros(1);
			lwork=int32(-1);
			info=int32(-1);
			[jobz, uplo, n, vec, lda, val, work, lwork, info]=lapack_dsyev(jobz, uplo, n, a, lda, w, work, lwork, info);
			lwork=int32(work(1));
			work=zeros(double(lwork),1);
			[jobz, uplo, n, vec, lda, val, work, lwork, info]=lapack_dsyev(jobz, uplo, n, a, lda, w, work, lwork, info);
		else
			w=single(zeros(double(n),1));
			work=single(zeros(1));
			lwork=int32(-1);
			info=int32(-1);
			[jobz, uplo, n, vec, lda, val, work, lwork, info]=lapack_ssyev(jobz, uplo, n, a, lda, w, work, lwork, info);
			lwork=int32(work(1));
			work=single(zeros(double(lwork),1));
			[jobz, uplo, n, vec, lda, val, work, lwork, info]=lapack_ssyev(jobz, uplo, n, a, lda, w, work, lwork, info);
		end;
	end;

repo: Aero-matlab


text:
function [ y, m, d ] = inflate_common ( y, m, d )

%*****************************************************************************80
%
%% INFLATE_COMMON "inflates" dates in the Common Calendar transition month.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 March 2001
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, integer M, integer D, the YMD date.
%
%    Output, integer Y, integer M, integer D, the YMD date.
%
  if ( y == 1582 )
    if ( m == 10 )
      if ( 5 <= d )
        d = d + 10;
      end
    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [value,isterminal,direction,eventDesc] = ma_NodeEvents(T,Y, node, bodyInfo, refBody, celBodyData)
    if(strcmpi(T, 'getInfo'))
        value = node;
        isterminal = bodyInfo;
        direction = refBody;
        eventDesc = [];

        return
    end

    if(~isempty(refBody) && bodyInfo.id ~= refBody.id)
        value = [];
        isterminal = [];
        direction = [];
        eventDesc = {};

        return
    end

    rVectECI = Y(1:3);

    topLevelBodyInfo = getTopLevelCentralBody(celBodyData);
    bodyInfoSun = topLevelBodyInfo;
    [rVectECI,vVectECI] = getAbsPositBetweenSpacecraftAndBody(T, rVectECI, bodyInfoSun, bodyInfo, celBodyData);

    [lat, ~, ~] = getLatLongAltFromInertialVect(T, -rVectECI, bodyInfo, vVectECI);

    value = lat;
    isterminal = 1;

    if(strcmpi(node,'asc'))
        direction = 1;
    else
        direction = -1;
    end

    eventDesc = [upper(node), ' Node'];
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function f = vis_derived_filter(model, tree)

% AUTORIGHTS
% -------------------------------------------------------
% Copyright (C) 2011-2012 Ross Girshick
%
% This file is part of the voc-releaseX code
% (http://people.cs.uchicago.edu/~rbg/latent/)
% and is available under the terms of an MIT-like license
% provided in COPYING. Please retain this notice and
% COPYING if you use this file (or a portion of it) in
% your project.
% -------------------------------------------------------

conf = voc_config();

% indexes into info from get_detection_trees.cc
% replace with tree_mat_to_struct
N_PARENT      = 1;
N_IS_LEAF     = 2;
N_SYMBOL      = 3;
N_RULE_INDEX  = 4;
N_RHS_INDEX   = 5;
N_X           = 6;
N_Y           = 7;
N_L           = 8;
N_DS          = 9;
N_DX          = 10;
N_DY          = 11;
N_SCORE       = 12;
N_LOSS        = 13;
N_SZ          = 14;

rx = tree(N_X, 1);
ry = tree(N_Y, 1);
rl = tree(N_L, 1);

f = zeros([0 0 conf.features.dim]);
off_x = 0;
off_y = 0;

for i = 2:size(tree,2)
  s = tree(N_SYMBOL, i);
  if model.symbols(s).type == 'T'
    x = off_x + tree(N_X, i) - rx;
    y = off_y + tree(N_Y, i) - ry;
    l = tree(N_L, i) - rl;

    pad = [abs(min(0, [y x])) 0];
    f = padarray(f, pad, 0, 'pre');
    if pad(1) > 0
      off_y = off_y + pad(1);
    end
    if pad(2) > 0
      off_x = off_x + pad(2);
    end

    w = model_get_block(model, model.filters(model.symbols(s).filter));
    wsz = size(w);
    fsz = size(f);
    req_fsz = [off_y + y + wsz(1), off_x + x + wsz(2), wsz(3)];
    pad = max(0, req_fsz - fsz);
    f = padarray(f, pad, 0, 'post');
    f(off_y+1+y:off_y+1+y+wsz(1)-1, off_x+1+x:off_x+1+x+wsz(2)-1, :) = ...
      f(off_y+1+y:off_y+1+y+wsz(1)-1, off_x+1+x:off_x+1+x+wsz(2)-1, :) + w;
  end
end

visualizeHOG(max(0, f));

author: rbgirshick

repo: voc-dpm

save_path: github-repos/MATLAB/rbgirshick-voc-dpm


text:
function geometry_test0755 ( )

%*****************************************************************************80
%
%% TEST0755 tests POLYGON_1_2D, POLYGON_X_2D, POLYGON_Y_2D, POLYGON_XX_2D, POLYGON_XY_2D and POLYGON_YY_2D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 February 2009
%
%  Author:
%
%    John Burkardt
%
  n = 4;
  dim_num = 2;

  v = [ ...
    0.0, 0.0; ...
    1.0, 0.0; ...
    1.0, 1.0; ...
    0.0, 1.0 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST0755\n' );
  fprintf ( 1, '  For a polygon in 2D:\n' );
  fprintf ( 1, '  POLYGON_1_2D integrates 1\n' );
  fprintf ( 1, '  POLYGON_X_2D integrates X\n' );
  fprintf ( 1, '  POLYGON_Y_2D integrates Y\n' );
  fprintf ( 1, '  POLYGON_XX_2D integrates X*X\n' );
  fprintf ( 1, '  POLYGON_XY_2D integrates X*Y\n' );
  fprintf ( 1, '  POLYGON_YY_2D integrates Y*Y\n' );

  r8mat_transpose_print ( dim_num, n, v, '  The polygon vertices:' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  F(X,Y)    Integral\n' );
  fprintf ( 1, '\n' );

  result = polygon_1_2d ( n, v );
  fprintf ( 1, '    1  %f\n', result );

  result = polygon_x_2d ( n, v );
  fprintf ( 1, '    X  %f\n', result );

  result = polygon_y_2d ( n, v );
  fprintf ( 1, '    Y  %f\n', result );

  result = polygon_xx_2d ( n, v );
  fprintf ( 1, '  X*X  %f\n', result );

  result = polygon_xy_2d ( n, v );
  fprintf ( 1, '  X*Y  %f\n', result );

  result = polygon_yy_2d ( n, v );
  fprintf ( 1, '  Y*Y  %f\n', result );

  return
end


text:
function [ output_args ] = Makeparadigm_MI( varargin )
% Makeparadigm_MI (Experimental paradigm):
%
% Description:
%   Basic motor imagery experiment paradigm using psychtoolbox.
%   It shows a cross, an arrow, and blank screen alternately.
%
% Example:
%   Makeparadigm_MI({'time_cross',1.5;'time_sti',60;'time_blank',5;'num_trial',1;'num_class',3});
%
% Input: (Nx2 size, cell-type)
%   time_cross  - time for concentration, showing a cross [s]
%   time_sti    - time for a stimulus, showing an arrow (right, left, or down) [s]
%   time_blank  - time for rest, showing nothing but gray screen [s]
%   num_trial   - number of trials per class
%   num_class   - number of class you want, 1 to 3, (right, left, and foot)
%   num_screen  - number of screen to show
%   screen_size - size of window showing experimental stimulus
%                 'full' or matrix (e.g.[0 0 300 300])
%

opt=opt_cellToStruct(varargin{:});

%% default setting
if ~isfield(opt,'time_sti'),    time_sti=4;      else time_sti=opt.time_sti;      end
if ~isfield(opt,'time_cross'),  time_cross=3;    else time_cross=opt.time_cross;  end
if ~isfield(opt,'time_blank'),  time_blank=3;    else time_blank=opt.time_blank;  end
if ~isfield(opt,'num_trial'),   num_trial=50;    else num_trial=opt.num_trial;    end
% if ~isfield(opt,'time_jitter'), time_jitter=0.1; else time_jitter=opt.time_jitter;end
if ~isfield(opt,'num_class'),   num_class=3;     else num_class=opt.num_class;    end
% screenNumber=2;

%% screen setting
screens=Screen('Screens');
if ~isfield(opt,'num_screen'),screenNumber=max(screens); else screenNumber=opt.num_screen; end
if ~isfield(opt,'screen_size'),screen_size='full'; else screen_size=opt.screen_size; end
% if ~isfield(opt,'screen_type'),screen_type='window'; else screen_type=opt.screen_type; end

%% beep setting
beep='on';
freq=22000;
beepLengthSecs=0.5;

%% trigger setting
global IO_ADDR IO_LIB;
IO_ADDR=hex2dec('D010');
IO_LIB=which('inpoutx64.dll');

%% jittering
a=-1;b=1;
% association with + time jittering
jitter = a + (b-a).*rand(num_trial*num_class,1);
jitter=jitter*0.5;

%% image load
currentFile = mfilename( 'fullpath' );
[pathstr,~,~] = fileparts( currentFile );
img_right=imread(fullfile(pathstr, "..", '\Stimulus\right.jpg'));
img_left=imread(fullfile(pathstr, "..",'\Stimulus\\left.jpg'));
img_down=imread(fullfile(pathstr, "..",'\Stimulus\down.jpg'));
% img_cross=imread('\Stimulus\cross.jpg');

%% beep sound
if strcmp(beep,'on')
    [beep,samplingRate] = MakeBeep(500,beepLengthSecs,freq);
    Snd('Open');
    sound=1;
else
    sound=0;
end
% Screen('Preference', 'SkipSyncTests', 1);

%% psychtoolbox setting
gray=GrayIndex(screenNumber);
% screenRes = [0 0 300 300];
if strcmp(screen_size,'full')
    [w, wRect]=Screen('OpenWindow',screenNumber, gray);
else
    [w, wRect]=Screen('OpenWindow',screenNumber, gray, screen_size);
end

%% order of stimulus (random)
for i=1:num_class
    a1(i,1:num_trial)=i;
end
% a1=Shuffle(a1);
% [t s]=size(a1);
% sti_stack=reshape(a1,1,t*s);
sti_stack=Shuffle(a1(:)'); % smkim

% click to start:
Screen('TextSize',w, 50);
DrawFormattedText(w,'Mouse click to start MI experiment \n\n (Press s to pause, esc to stop)','center','center',[0 0 0]);
Screen('Flip', w);
GetClicks(w);
ppTrigger(111); % START

escapeKey = KbName('esc');
waitKey=KbName('s');
p_close=0;

%% Eyes open/closed
% run(fullfile(pathstr,"..", '\Artifact\eyesOpenClosed')) % script
run('Paradigm\Artifact\eyesOpenClosed') % script

%% fixation cross
[X,Y] = RectCenter(wRect);
FixationSize = 20;
FixCross = [X-1,Y-FixationSize,X+1,Y+FixationSize;X-FixationSize,Y-1,X+FixationSize,Y+1];

%% paradigm start
    Screen('FillRect', w, [0 0 0], FixCross');
    Screen('Flip', w);
    WaitSecs(2.5)

for num_stimulus=1:length(sti_stack)

    if num_stimulus==length(sti_stack)/2
        Screen('TextSize',w, 50);
        DrawFormattedText(w,'Rest\n\n(Pause the brain vision)','center','center',[0 0 0]);
        Screen('Flip', w);
        GetClicks(w);
        DrawFormattedText(w,'(Resume recording)','center','center',[0 0 0]);
        Screen('Flip', w);
        GetClicks(w);
        DrawFormattedText(w,'Click to continue the experiment','center','center',[0 0 0]);
        Screen('Flip', w);
        GetClicks(w);
    end

    start=GetSecs;
    while GetSecs < start+time_blank+jitter(num_stimulus)-2
        [ keyIsDown, seconds, keyCode ] = KbCheck;
        if keyIsDown
            if keyCode(escapeKey)
%                 ShowCursor;
%                 p_close=1;
                Screen('CloseAll');
                fclose('all');
                return
            elseif keyCode(waitKey)
                warning('stop')
                GetClicks(w);
                Screen('Close',tex1);

            end
        end
        pause(0.1);
    end

%     if p_close
%         break;
%     end

    Screen('FillRect', w, [0 0 0], FixCross');
    Screen('Flip', w);
    if sound
        Snd('Play',beep);
    end

    WaitSecs(1)
    ppTrigger(5);
    WaitSecs(time_cross)

    switch sti_stack(num_stimulus)
        case 1 % right class
            ppTrigger(1);
            image=img_right;
            tex1=Screen('MakeTexture', w, image );
            Screen('DrawTexture', w, tex1);
            Screen('FillRect', w, [0 0 0], FixCross');
            Screen('Flip', w);

            WaitSecs(time_sti);
            Screen('Close',tex1);
%             ppTrigger(11);

        case 2 % left class
            ppTrigger(2);
            image=img_left;
            tex1=Screen('MakeTexture', w, image );
            Screen('DrawTexture', w, tex1);
            Screen('FillRect', w, [0 0 0], FixCross');
            Screen('Flip', w);

            WaitSecs(time_sti);
            Screen('Close',tex1);
%             ppTrigger(22);

        case 3 % foot class
            ppTrigger(3);
            image=img_down;
            tex1=Screen('MakeTexture', w, image );
            Screen('DrawTexture', w, tex1);
            Screen('FillRect', w, [0 0 0], FixCross');
            Screen('Flip', w);

            WaitSecs(time_sti);
            Screen('Close',tex1);
%             ppTrigger(33);
    end
    num_stimulus

    Screen('Flip', w);
    WaitSecs(time_blank)

end

ppTrigger(222);
Screen('TextSize',w, 50);
DrawFormattedText(w, 'Thank you', 'center', 'center', [0 0 0]);
Screen('Flip', w);
WaitSecs(2);
Screen('CloseAll');
ShowCursor;
fclose('all');
Priority(0);

author: PatternRecognition

repo: OpenBMI


text:
function [shortestPaths, totalCosts] = kShortestPath(netCostMatrix, source, destination, k_paths)
% Function kShortestPath(netCostMatrix, source, destination, k_paths)
% returns the K first shortest paths (k_paths) from node source to node destination
% in the a network of N nodes represented by the NxN matrix netCostMatrix.
% In netCostMatrix, cost of 'inf' represents the 'absence' of a link
% It returns
% [shortestPaths]: the list of K shortest paths (in cell array 1 x K) and
% [totalCosts]   : costs of the K shortest paths (in array 1 x K)
%==============================================================
% Meral Shirazipour
% This function is based on Yen's k-Shortest Path algorithm (1971)
% This function calls a slightly modified function dijkstra()
% by Xiaodong Wang 2004.
% * netCostMatrix must have positive weights/costs
%==============================================================
%  DATE :           December 9 decembre 2009
%  Last Updated:    August 2 2010; January 6 2011; August 2 2011
%  ----Changes April 2 2010:----
%  1-previous version(9/12/2009)did not handle some exceptions which should
%    have returned empty matrices for the return values (added lines 20 and 29)
%  2-includes the changes proposed by Darren Rowland
%  ----Changes January 6 2011:----
%  1-fixed a bug reported by Babak Zafari that prevented from finding ALL
%    the shortest paths in large networks
%==============================================================
if source > size(netCostMatrix,1) || destination > size(netCostMatrix,1)
    warning('The source or destination node are not part of netCostMatrix');
    shortestPaths=[];
    totalCosts=[];
else
    %---------------------INITIALIZATION---------------------
    k=1;
    [path cost] = dijkstra(netCostMatrix, source, destination);
    %P is a cell array that holds all the paths found so far:
    if isempty(path)
        shortestPaths=[];
        totalCosts=[];
    else
        path_number = 1;
        P{path_number,1} = path; P{path_number,2} = cost;
        current_P = path_number;
        %X is a cell array of a subset of P (used by Yen's algorithm below):
        size_X=1;
        X{size_X} = {path_number; path; cost};

        %S path_number x 1
        S(path_number) = path(1); %deviation vertex is the first node initially

        %***********************
        % K = 1 is the shortest path returned by dijkstra():
        shortestPaths{k} = path ;
        totalCosts(k) = cost;
        %***********************

        %--------------------------------------------------------
        while (k < k_paths   &&   size_X ~= 0 )
            %remove P from X
            for i=1:length(X)
                if  X{i}{1} == current_P
                    size_X = size_X - 1;
                    X(i) = [];%delete cell
                    break;
                end
            end

            %---------------------------------------
            P_ = P{current_P,1}; %P_ is current P, just to make is easier for the notations

            %Find w in (P_,w) in set S, w was the dev vertex used to found P_
            w = S(current_P);
            for i = 1: length(P_)
                if w == P_(i)
                    w_index_in_path = i;
                end
            end

            for index_dev_vertex= w_index_in_path: length(P_) - 1   %index_dev_vertex is index in P_ of deviation vertex
                temp_netCostMatrix = netCostMatrix;
                %------
                %Remove vertices in P before index_dev_vertex and there incident edges
                for i = 1: index_dev_vertex-1
                    v = P_(i);
                    temp_netCostMatrix(v,:)=inf;
                    temp_netCostMatrix(:,v)=inf;
                end
                %------
                %remove incident edge of v if v is in shortestPaths (K) U P_  with similar sub_path to P_....
                SP_sameSubPath=[];
                index =1;
                SP_sameSubPath{index}=P_;
                for i = 1: length(shortestPaths)
                    if length(shortestPaths{i}) >= index_dev_vertex
                        if P_(1:index_dev_vertex) == shortestPaths{i}(1:index_dev_vertex)
                            index = index+1;
                            SP_sameSubPath{index}=shortestPaths{i};
                        end
                    end
                end
                v_ = P_(index_dev_vertex);
                for j = 1: length(SP_sameSubPath)
                    next = SP_sameSubPath{j}(index_dev_vertex+1);
                    temp_netCostMatrix(v_,next)=inf;
                end
                %------

                %get the cost of the sub path before deviation vertex v
                sub_P = P_(1:index_dev_vertex);
                cost_sub_P=0;
                for i = 1: length(sub_P)-1
                    cost_sub_P = cost_sub_P + netCostMatrix(sub_P(i),sub_P(i+1));
                end

                %call dijkstra between deviation vertex to destination node
                [dev_p c] = dijkstra(temp_netCostMatrix, P_(index_dev_vertex), destination);
                if ~isempty(dev_p)
                    path_number = path_number + 1;
                    P{path_number,1} = [sub_P(1:end-1) dev_p] ;  %concatenate sub path- to -vertex -to- destination
                    P{path_number,2} =  cost_sub_P + c ;

                    S(path_number) = P_(index_dev_vertex);

                    size_X = size_X + 1;
                    X{size_X} = {path_number;  P{path_number,1} ;P{path_number,2} };
                else
                    %warning('k=%d, isempty(p)==true!\n',k);
                end
            end
            %---------------------------------------
            %Step necessary otherwise if k is bigger than number of possible paths
            %the last results will get repeated !
            if size_X > 0
                shortestXCost= X{1}{3};  %cost of path
                shortestX= X{1}{1};        %ref number of path
                for i = 2 : size_X
                    if  X{i}{3} < shortestXCost
                        shortestX= X{i}{1};
                        shortestXCost= X{i}{3};
                    end
                end
                current_P = shortestX;
                %******
                k = k+1;
                shortestPaths{k} = P{current_P,1};
                totalCosts(k) = P{current_P,2};
                %******
            else
                %k = k+1;
            end
        end
    end
end
%--------------------------------------------------------

author: Sable

repo: mcbench-benchmarks


text:
function [M,indi,indj] = randblock(M,S)
% RANDBLOCK - randomize blocks of a matrix
%    R = RANDBLOCK(M,S) randomizes the matrix M by dividing M into
%    non-overlapping blocks of the size specified by S, and shuffling these
%    blocks. M can be a N-D matrix.
%    The number of elements in S should match the number of dimensions of
%    M, or S can be a scalar specifying a S-by-S-by-S-by ... block size.
%    S should contain positive integers. TThe size of M in any dimension
%    should be an integer number of times the specified size of the block in
%    that dimension (e.g., if size(M,1) equals 6, S(dim) can be 1,2,3, or 6).
%
%    [R,I,J] = RANDBLOCK(...) also returns indices I and J, so that R equals A(I)
%    and R(J) equals A.
%
%    M can be a numerical or cell array.
%
%    Examples:
%      % Shuffle blocks of 3 elements of a 15-element vector
%        M = 1:15 ;
%        randblock(M,3)
%
%      % Scramble a 2D matrix
%        M = reshape(1:24,6,[]) ;
%        randblock(M,[3 2])   % randomize the position of the four 3-by-2 blocks
%        randblock(M,2)       % randomize the position of the six 2-by-2 blocks
%
%      % Scramble a 3D volume
%        M = reshape(1:64,[4 8 2]) ;
%        randblock(M,[2 4 2]) % randomize the position of the four
%                             % 2-by-4-by-2 sub-volumes
%
%      % Scramble a cell matrix
%        M = {'1','a','bb','1c' ; '2a',[3 4 5],'2c','2dd'} ;
%        randblock(M,[2 2])   % randomize the position of the four 2-by-2 blocks
%
%      % Scramble a RGB image Z, and retrieve the original using the
%      % indices
%        Z = peaks(200) ; Z = cat(3,Z,flipud(Z), -(circshift(Z.',[100,100]))) ;
%        Z = (Z - min(Z(:))) ; Z = Z./ max(Z(:)) ;
%        [Z2,I,J] = randblock(Z,[25,25,size(Z,3)]) ; % Scramble 25-by-25 blocks
%        subplot(2,2,1) ; image(Z) ; title('Original image') ;
%        subplot(2,2,2) ; image(Z2) ; title('25x25 scrambled image') ;
%        subplot(2,2,3) ; image(Z2(J)) ; title('Z2(J) = original image') ;
%        subplot(2,2,4) ; image(Z(I)) ; title('Z(I) = scrambled image') ;
%        set (gcf,'Name','[Z2,I,J] = randblock(Z,[25,25,size(Z,3)])') ;
%
%  See also RAND, RANDPERM, MAT2CELL
%       and RANDSWAP, SHAKE on the File Exchange

% for Matlab R13 and up (works on R2010b)
% version 2.3 (mar 2011)
% (c) Jos van der Geest
% email: jos@jasen.nl

% History:
% Created: nov 2007
% 2.0 (dec 2007) implemented cell matrix capabilities
% 2.1 (dec 2007) added several examples for the File Exchange
% 2.2 (dec 2007) modified image example, remove minor m-lint messages
% 2.3 (mar 2011) fixed a small error in the examples; tested in 2010b

% check if the size of the blocks is specified for each dimension of M
NdimM = ndims(M) ;
sizM = size(M) ;  % size of the matrix
nS = numel(S) ;

if nS==1,
    S = repmat(S,1,NdimM) ; % scalar expansion
elseif nS ~= NdimM,
    error('Number of elements of S should equal the number of dimensions of M.') ;
end

% S should contain positive integers only
if ~isnumeric(S) || any(S ~= fix(S)) || any(S<1),
    error('S should be a numeric array with positive integer values.') ;
end

% vectors are 2D and S can be a scalar in this case:
% the singleton dimension of M should be taken care of.
if NdimM==2 && any(sizM==1) && nS==1 && all(S~=1),
    S(sizM==1) = 1 ;
end

if all(S==1),
    % Each element of M is a single block. We can just randomize the linear
    % indices.
    indj = randperm(numel(M)) ;
else

    nb = sizM ./ S ;  % how many blocks in each dimension, this should be an integer
    B = cell(NdimM,1) ;
    % for each dimension of M:
    for i=1:NdimM,
        % check if the number of blocks is an integer
        if nb(i) ~= fix(nb(i)),
            error(['Size mismatch: the size of the matrix M (= %d) in dimension %d\n' ...
                    'is not an integer number of times the specified blocksize (= %d).'],...
                sizM(i),i,S(i)) ;
        else
            % expand the size in that dimension. B is used for mat2cell
            B{i} = repmat(S(i),1,nb(i)) ;
        end
    end

    % M can be a numerical or cell array. By randomizing the (linear)
    % indices, we do not have to worry about the actual contents of M.
    indj = reshape(1:numel(M),sizM) ;
    % convert to a cell array. Each cell contains a block.
    C = mat2cell(indj,B{:}) ;
    % randomize the order of the cells
    C(randperm(numel(C))) = C ;
    % convert back from cell array
    indj = cell2mat(C) ;
end

% use these randomized block indices to shuffle M
M(indj) = M ;

if nargout>1,
    indi(indj) = 1:numel(indj) ;
    indi = reshape(indi,sizM) ;
    indj = reshape(indj,sizM) ;
end


text:
% This is the testing demo of FFDNet for denoising noisy grayscale images corrupted by
% AWGN.
%
% To run the code, you should install Matconvnet first. Alternatively, you can use the
% function `vl_ffdnet_matlab` to perform denoising without Matconvnet.
%
% "FFDNet: Toward a Fast and Flexible Solution for CNN based Image Denoising"
%  2018/03/23
% If you have any question, please feel free to contact with me.
% Kai Zhang (e-mail: cskaizhang@gmail.com)

%clear; clc;
format compact;
global sigmas; % input noise level or input noise level map
addpath(fullfile('utilities'));

folderModel = 'models';
folderTest  = 'H:\matlabH';
folderResult= 'results';
imageSets   = {'BSD68','Set12'}; % testing datasets
setTestCur  = imageSets{1};      % current testing dataset

showResult  = 1;
useGPU      = 1; % CPU or GPU. For single-threaded (ST) CPU computation, use "matlab -singleCompThread" to start matlab.
pauseTime   = 0;

imageNoiseSigma = 25;  % image noise level
inputNoiseSigma = 25;  % input noise level

folderResultCur       =  fullfile(folderResult, [setTestCur,'_',num2str(imageNoiseSigma),'_',num2str(inputNoiseSigma)]);
if ~isdir(folderResultCur)
    mkdir(folderResultCur)
end

load(fullfile('data','FFDNet_gray','FFDNet_gray-epoch-1.mat'));
net.layers = net.layers(1:end-1);
net = vl_simplenn_tidy(net);

% for i = 1:size(net.layers,2)
%     net.layers{i}.precious = 1;
% end

if useGPU
    net = vl_simplenn_move(net, 'gpu') ;
end

% read images
ext         =  {'*.jpg','*.png','*.bmp'};
filePaths   =  [];
for i = 1 : length(ext)
    filePaths = cat(1,filePaths, dir(fullfile(folderTest,setTestCur,ext{i})));
end

% PSNR and SSIM
PSNRs = zeros(1,length(filePaths));
SSIMs = zeros(1,length(filePaths));

for i = 1:length(filePaths)

    % read images
    label = imread(fullfile(folderTest,setTestCur,filePaths(i).name));
    [w,h,~]=size(label);
    if size(label,3)==3
        label = rgb2gray(label);
    end

    [~,nameCur,extCur] = fileparts(filePaths(i).name);
    label = im2double(label);

    % add noise
    randn('seed',0);
    noise = imageNoiseSigma/255.*randn(size(label));
    input = single(label + noise);

    if mod(w,2)==1
        input = cat(1,input, input(end,:)) ;
    end
    if mod(h,2)==1
        input = cat(2,input, input(:,end)) ;
    end

    % tic;
    if useGPU
        input = gpuArray(input);
    end

    % set noise level map
    sigmas = inputNoiseSigma/255; % see "vl_simplenn.m".

    % perform denoising
    res    = vl_simplenn(net,input,[],[],'conserveMemory',true,'mode','test'); % matconvnet default
    % res    = vl_ffdnet_concise(net, input);    % concise version of vl_simplenn for testing FFDNet
    % res    = vl_ffdnet_matlab(net, input); % use this if you did  not install matconvnet; very slow

    % output = input - res(end).x; % for 'model_gray.mat'
    output = res(end).x;

    if mod(w,2)==1
        output = output(1:end-1,:);
        input  = input(1:end-1,:);
    end
    if mod(h,2)==1
        output = output(:,1:end-1);
        input  = input(:,1:end-1);
    end

    if useGPU
        output = gather(output);
        input  = gather(input);
    end
    % toc;
    % calculate PSNR, SSIM and save results
    [PSNRCur, SSIMCur] = Cal_PSNRSSIM(im2uint8(label),im2uint8(output),0,0);
    if showResult
        imshow(cat(2,im2uint8(input),im2uint8(label),im2uint8(output)));
        title([filePaths(i).name,'    ',num2str(PSNRCur,'%2.2f'),'dB','    ',num2str(SSIMCur,'%2.4f')])
        %imwrite(im2uint8(output), fullfile(folderResultCur, [nameCur, '_' num2str(imageNoiseSigma,'%02d'),'_' num2str(inputNoiseSigma,'%02d'),'_PSNR_',num2str(PSNRCur*100,'%4.0f'), extCur] ));
        drawnow;
        pause(pauseTime)
    end
    disp([filePaths(i).name,'    ',num2str(PSNRCur,'%2.2f'),'dB','    ',num2str(SSIMCur,'%2.4f')])
    PSNRs(i) = PSNRCur;
    SSIMs(i) = SSIMCur;
end

disp([mean(PSNRs),mean(SSIMs)]);

author: cszn

repo: FFDNet

save_path: github-repos/MATLAB/cszn-FFDNet


text:
function d = DistanceToSE3(mat)
% *** CHAPTER 3: RIGID-BODY MOTIONS ***
% Takes mat: A 4x4 matrix.
% Returns the Frobenius norm to describe the distance of mat from the SE(3)
% manifold.
% Compute the determinant of matR, the top 3x3 submatrix of mat. If
% det(matR) <= 0, return a large number. If det(matR) > 0, replace the top
% 3x3 submatrix of mat with matR' * matR, and set the first three entries
% of the fourth column of mat to zero. Then return norm(mat - I).
% Example Inputs:
%
% clear; clc;
% mat = [1.0, 0.0,   0.0,  1.2;
%        0.0, 0.1, -0.95,  1.5;
%        0.0, 1.0,   0.1, -0.9;
%        0.0, 0.0,   0.1, 0.98];
% d = DistanceToSE3(mat)
%
% Output:
% d =
%     0.1349

[R, p] = TransToRp(mat);
if det(R) > 0
	d = norm([R' * R, [0; 0; 0]; mat(4, :)] - eye(4), 'fro');
else
    d = 1e+9;
end
end

author: ShuoYangRobotics

repo: QuadrupedSim


text:
function cn=adodbcn(cnstr,cto)
% cn = adodbcn(cnstr,[cto])
%
% Connects to ADO OLEDB using the Microsoft ActiveX Data Source Control
%
% Inputs:
%   cnstr,  str containing information for connecting to data source
%   cto,    CommandTimeout in seconds (default=60 seconds if unspecified)
%
% Output:
%   cn,     connection
%
% Notes: Refer to adodb_demo.m and adodbcnstr.m for example connection
% strings for Oracle, SQL, and MS Access databases
%
% This code basis on Tim Myers code (oledb*.m) and use ADO OLE DB instead
% of OWC - Office Web Component
%
% Martin Furlan
% martin.furlan@iskra-ae.com
% January 2007

%create activeX control
cn = actxserver('ADODB.Connection');
set(cn,'CursorLocation',3);

%Open connection
invoke(cn,'Open', cnstr);

%Specify connection timeout if provided
if nargin>1
    set(cn,'CommandTimeout',cto);
else
    set(cn,'CommandTimeout',60);    %default
end


text:
classdef AddMassToTankActionOptimVar < AbstractOptimizationVariable
    %AddMassToTankActionOptimVar Summary of this class goes here
    %   Detailed explanation goes here

    properties
        varObj(1,1) AddMassToTankAction = AddMassToTankAction();

        lb(1,1) double = 0;
        ub(1,1) double = 0;

        useTf(1,1) = false;
    end

    methods
        function obj = AddMassToTankActionOptimVar(varObj)
            obj.varObj = varObj;
            obj.varObj.optVar = obj;

            obj.id = rand();
        end

        function x = getXsForVariable(obj)
            x = [];

            if(obj.useTf)
                x = obj.varObj.massToAdd;
            end
        end

        function [lb, ub] = getBndsForVariable(obj)
            lb = obj.lb(logical(obj.useTf));
            ub = obj.ub(logical(obj.useTf));
        end

        function [lb, ub] = getAllBndsForVariable(obj)
            lb = obj.lb;
            ub = obj.ub;
        end

        function setBndsForVariable(obj, lb, ub)
            obj.lb = lb;
            obj.ub = ub;
        end

        function useTf = getUseTfForVariable(obj)
            useTf = logical(obj.useTf);
        end

        function setUseTfForVariable(obj, useTf)
            obj.useTf = logical(useTf);
        end

        function updateObjWithVarValue(obj, x)
            obj.varObj.massToAdd = x;
        end

        function nameStrs = getStrNamesOfVars(obj, evtNum, varLocType)
            nameStrs = {sprintf('Event %i Add Mass To Tank Action', evtNum)};
        end
    end
end


text:
function trsoutput = trs_tCG_cached(problem, trsinput, options, storedb, key)
% Truncated (Steihaug-Toint) Conjugate-Gradient method with caching.
%
% minimize <eta,grad> + .5*<eta,Hess(eta)>
% subject to <eta,eta>_[inverse precon] <= Delta^2
%
% function trsoutput = trs_tCG_cached(problem, trsinput, options, storedb, key)
%
% trs_tCG_cached stores information (when options.trscache == true)
% which can help avoid redundant computations (using tCG_rejectedstep)
% upon step rejection by trustregions compared to trs_tCG at the cost of
% using extra memory.
%
% Inputs:
%   problem: Manopt optimization problem structure
%   trsinput: structure with the following fields:
%       x: point on the manifold problem.M
%       fgradx: gradient of the cost function of the problem at x
%       Delta = trust-region radius
%   options: structure containing options for the subproblem solver
%   storedb, key: manopt's caching system for the point x
%
% Options specific to this subproblem solver:
%   kappa (0.1)
%       kappa convergence tolerance.
%       kappa > 0 is the linear convergence target rate: trs_tCG_cached
%       terminates early if the residual was reduced by a factor of kappa.
%   theta (1.0)
%       theta convergence tolerance.
%       1+theta (theta between 0 and 1) is the superlinear convergence
%       target rate. trs_tCG_cached terminates early if the residual
%       was reduced by a power of 1+theta.
%   mininner (1)
%       Minimum number of inner iterations.
%   maxinner (problem.M.dim())
%       Maximum number of inner iterations.
%   trscache (true)
%       Set to false if no caching for the trs_tCG_cached is desired. It is
%       default true to improve computation time if there are many step
%       rejections in trustregions. Setting trscache to false can reduce
%       memory usage.
%   memorytCG_warningtol (1000)
%       Tolerance memory value in MB before issuing warning when
%       trscache = true.
%       The default is 1GB but this value can be increased depending on the
%       user's machine. To disable the warning completely use:
%       warning('off', 'manopt:trs_tCG_cached:memory')
%
% Output: the structure trsoutput contains the following fields:
%   eta: approximate solution to the trust-region subproblem at x
%   Heta: Hess f(x)[eta] -- this is necessary in the outer loop, and it
%       is often naturally available to the subproblem solver at the
%       end of execution, so that it may be cheaper to return it here.
%   limitedbyTR: true if a boundary solution is returned
%   printstr: logged information to be printed by trustregions.
%   stats: structure with the following statistics:
%           numinner: number of inner loops before returning
%           hessvecevals: number of Hessian calls issued
%           memorytCG_MB: memory of store_iters and store_last in MB
%
% Stored Information:
%   store_iters: a struct array with enough information to compute the next
%       step upon step rejection when the algorithm exits due to negative
%       curvature or trust-region radius violation.
%   store_last: an additional struct to store_iters to compute the next
%       step upon step rejection when the algorithm exits but not due to
%       negative curvature or trust-region radius violation
%
%
% trs_tCG_cached can also be called in the following way (by trustregions)
% to obtain part of the header to print and an initial stats structure:
%
% function trsoutput = trs_tCG_cached([], [], options)
%
% In this case trsoutput contains the following fields:
%   printheader: subproblem header to be printed before the first pass of
%       trustregions
%   initstats: struct with initial values for stored stats in subsequent
%       calls to trs_tCG_cached. Used in the first call to savestats
%       in trustregions to initialize the info struct properly.
%
% See also: trustregions trs_tCG tCG_rejectedstep

% This file is part of Manopt: www.manopt.org.
% This code is an adaptation to Manopt of the original GenRTR code:
% RTR - Riemannian Trust-Region
% (c) 2004-2007, P.-A. Absil, C. G. Baker, K. A. Gallivan
% Florida State University
% School of Computational Science
% (http://www.math.fsu.edu/~cbaker/GenRTR/?page=download)
% See accompanying license file.
% The adaptation was executed by Nicolas Boumal.
%
% Change log:
%
%   VL June 24, 2022:
%       trs_tCG_cached by default stores information at each iteration
%       compared to trs_tCG.
%       This can be useful for the next call to trs_tCG_cached and the work
%       is passed to tCG_rejectedstep rather than the normal tCG loop.

% See trs_tCG for references to relevant equations in
% [CGT2000] Conn, Gould and Toint: Trust-region methods, 2000.

% trustregions only wants header and default values for stats.
if nargin == 3 && isempty(problem) && isempty(trsinput)
    trsoutput.printheader = '';
    if options.verbosity == 2
        trsoutput.printheader = sprintf('%9s   %9s   %9s   %s', ...
                            'numinner', 'hessvec', 'numstored', ...
                            'stopreason');
    elseif options.verbosity > 2
        trsoutput.printheader = sprintf('%9s   %9s   %9s   %9s   %s', ...
                            'numinner', 'hessvec', 'numstored', ...
                            'memtCG_MB', 'stopreason');
    end
    trsoutput.initstats = struct('numinner', 0, 'hessvecevals', 0, ...
                   'memorytCG_MB', 0);
    return;
end

if isfield(options, 'useRand') && options.useRand
    warning('manopt:trs_tCG_cached:rand', ...
            ['options.useRand = true but @trs_tCG_cached ignores it.\n' ...
             'You may set options.subproblemsolver = @trs_tCG;\n', ...
             'Alternatively, set options.useRand = false;']);
end

x = trsinput.x;
Delta = trsinput.Delta;
grad = trsinput.fgradx;

inner   = @(u, v) problem.M.inner(x, u, v);
lincomb = @(a, u, b, v) problem.M.lincomb(x, a, u, b, v);
tangent = @(u) problem.M.tangent(x, u);

% Set local defaults here
localdefaults.kappa = 0.1;
localdefaults.theta = 1.0;
localdefaults.mininner = 1;
localdefaults.maxinner = problem.M.dim();
localdefaults.trscache = true;
localdefaults.memorytCG_warningtol = 1000;

% Merge local defaults with user options, if any
if ~exist('options', 'var') || isempty(options)
    options = struct();
end
options = mergeOptions(localdefaults, options);

% If the previous step was rejected and we want to use caching,
if ~trsinput.accept && options.trscache
    % Then check if there is cached information for the current point.
    store = storedb.get(key);
    if isfield(store, 'store_iters')
        % If so, use that cache to produce the same output as would have
        % been produced by running the code below (after the 'return'), but
        % without issuing Hessian-vector calls that were already issued.
        trsoutput = tCG_rejectedstep(problem, trsinput, options, store);
        return;
    end
end

% returned boolean to trustregions. true if we are limited by the TR
% boundary (returns boundary solution). Otherwise false.
limitedbyTR = false;

theta = options.theta;
kappa = options.kappa;

eta = problem.M.zerovec(x);
Heta = problem.M.zerovec(x);
r = grad;
e_Pe = 0;

r_r = inner(r, r);
norm_r = sqrt(r_r);
norm_r0 = norm_r;

% Precondition the residual.
z = getPrecon(problem, x, r, storedb, key);

% Compute z'*r.
z_r = inner(z, r);
d_Pd = z_r;

% Initial search direction (we maintain -delta in memory, called mdelta, to
% avoid a change of sign of the tangent vector.)
mdelta = z;
e_Pd = 0;

% If the Hessian or a linear Hessian approximation is in use, it is
% theoretically guaranteed that the model value decreases strictly
% with each iteration of trs_tCG. Hence, there is no need to monitor the model
% value. But, when a nonlinear Hessian approximation is used (such as the
% built-in finite-difference approximation for example), the model may
% increase. It is then important to terminate the trs_tCG iterations and return
% the previous (the best-so-far) iterate. The variable below will hold the
% model value.
%
% This computation could be further improved based on Section 17.4.1 in
% Conn, Gould, Toint, Trust Region Methods, 2000.
% If we make this change, then also modify trustregions to gather this
% value from trs_tCG rather than recomputing it itself.
model_fun = @(eta, Heta) inner(eta, grad) + .5*inner(eta, Heta);
model_value = 0;

% Pre-assume termination because j == end.
stopreason_str = 'maximum inner iterations';

% Track certain iterations in case step is rejected.
% store_iters tracks candidate etas with increasing squared
% norm relevant when limitedbyTR = true, or when <eta, Heta> <= 0
store_iters = struct('normsq', [], 'numinner', [], 'e_Pe', [], ...
    'd_Pd', [], 'e_Pd', [], 'd_Hd', [], 'eta', [], 'Heta', [], ...
    'mdelta', [], 'Hmdelta', []);

max_normsq = 0;

% only need to compute memory for one item in store_iters in Megabytes(MB)
peritermemory_MB = 0;

% total cached memory stored in MB
memorytCG_MB = 0;

% number of iterations where trs_tCG_cached stores information. This value
% will be length(store_iters) plus 1 if store_last is used.
numstored = 0;

% string that is printed by trustregions. For printing
% per-iteration information
printstr = '';

% Begin inner/trs_tCG loop.
for j = 1 : options.maxinner

    % This call is the computationally expensive step.
    Hmdelta = getHessian(problem, x, mdelta, storedb, key);

    % Compute curvature (often called kappa).
    d_Hd = inner(mdelta, Hmdelta);

    % Note that if d_Hd == 0, we will exit at the next "if" anyway.
    alpha = z_r/d_Hd;
    % <neweta,neweta>_P =
    % <eta,eta>_P + 2*alpha*<eta,delta>_P + alpha*alpha*<delta,delta>_P
    e_Pe_new = e_Pe + 2.0*alpha*e_Pd + alpha*alpha*d_Pd;

    if options.debug > 2
        fprintf('DBG:   (r,r)  : %e\n', r_r);
        fprintf('DBG:   (d,Hd) : %e\n', d_Hd);
        fprintf('DBG:   alpha  : %e\n', alpha);
    end

    if options.trscache
        % Selectively store info in store_iter.
        % next_smallest = (1/4^n Delta)^2 with n the smallest integer such
        % that max_normsq <= next_smallest.
        % We use this condition to only store relevant iterations in case
        % of rejection in trustregions.
        if max_normsq > 0
            next_smallest = (1/16)^floor(-(1/4)*(log2(max_normsq) - ...
                                                 log2(Delta^2))) * Delta^2;
        else
            next_smallest = 0;
        end

        if d_Hd <= 0 || e_Pe_new >= next_smallest
            numstored = numstored + 1;

            store_iters(numstored) = struct('normsq', e_Pe_new, 'numinner', ...
                             j, 'e_Pe', e_Pe, 'd_Pd', d_Pd, 'e_Pd', e_Pd,...
                             'd_Hd', d_Hd, 'eta', eta, 'Heta', Heta, ...
                             'mdelta', mdelta, 'Hmdelta', Hmdelta);
            max_normsq = e_Pe_new;

            % getSize for one entry in store_iters which will be the same
            % for all others.
            if peritermemory_MB == 0
                peritermemory_MB = getsize(store_iters(numstored))/1024^2;
            end

            memorytCG_MB = memorytCG_MB + peritermemory_MB;

            if memorytCG_MB > options.memorytCG_warningtol
                warning('manopt:trs_tCG_cached:memory', ...
                [sprintf('trs_tCG_cached will cache %.2f [MB] for at least one iteration of trustregions until a step is accepted.', memorytCG_MB) ...
                'If memory is limited turn off caching by options.trscache = false.\n' ...
                'To disable this warning: warning(''off'', ''manopt:trs_tCG_cached:memory'')']);
             end

        end
    end

    % Check against negative curvature and trust-region radius violation.
    % If either condition triggers, we bail out.
    if d_Hd <= 0 || e_Pe_new >= Delta^2
        % want
        %  ee = <eta,eta>_prec,x
        %  ed = <eta,delta>_prec,x
        %  dd = <delta,delta>_prec,x
        % Note (Nov. 26, 2021, NB): numerically, it might be better to call
        %   tau = max(real(roots([d_Pd, 2*e_Pd, e_Pe-Delta^2])));
        % This should be checked.
        % Also, we should safe-guard against 0/0: could happen if grad = 0.

        % store new struct containing all the required info in store_iter
        tau = (-e_Pd + sqrt(e_Pd*e_Pd + d_Pd*(Delta^2-e_Pe))) / d_Pd;
        if options.debug > 2
            fprintf('DBG:     tau  : %e\n', tau);
        end
        eta = lincomb(1,  eta, -tau,  mdelta);

        % If only a nonlinear Hessian approximation is available, this is
        % only approximately correct, but saves an additional Hessian call.
        Heta = lincomb(1, Heta, -tau, Hmdelta);

        % Technically, we may want to verify that this new eta is indeed
        % better than the previous eta before returning it (this is always
        % the case if the Hessian approximation is linear, but I am unsure
        % whether it is the case or not for nonlinear approximations.)
        % At any rate, the impact should be limited, so in the interest of
        % code conciseness (if we can still hope for that), we omit this.

        limitedbyTR = true;

        if d_Hd <= 0
            stopreason_str = 'negative curvature';
        else
            stopreason_str = 'exceeded trust region';
        end
        break;
    end

    % No negative curvature and eta_prop inside TR: accept it.
    e_Pe = e_Pe_new;
    new_eta  = lincomb(1,  eta, -alpha,  mdelta);

    % If only a nonlinear Hessian approximation is available, this is
    % only approximately correct, but saves an additional Hessian call.
    % TODO: this computation is redundant with that of r, L241. Clean up.
    new_Heta = lincomb(1, Heta, -alpha, Hmdelta);

    % Verify that the model cost decreased in going from eta to new_eta. If
    % it did not (which can only occur if the Hessian approximation is
    % nonlinear or because of numerical errors), then we return the
    % previous eta (which necessarily is the best reached so far, according
    % to the model cost). Otherwise, we accept the new eta and go on.
    new_model_value = model_fun(new_eta, new_Heta);
    if new_model_value >= model_value
        stopreason_str = 'model increased';
        break;
    end

    eta = new_eta;
    Heta = new_Heta;
    model_value = new_model_value; %% added Feb. 17, 2015

    % Update the residual.
    r = lincomb(1, r, -alpha, Hmdelta);

    % Compute new norm of r.
    r_r = inner(r, r);
    norm_r = sqrt(r_r);

    % Check kappa/theta stopping criterion.
    % Note that it is somewhat arbitrary whether to check this stopping
    % criterion on the r's (the gradients) or on the z's (the
    % preconditioned gradients). [CGT2000], page 206, mentions both as
    % acceptable criteria.
    if j >= options.mininner && norm_r <= norm_r0*min(norm_r0^theta, kappa)
        % Residual is small enough to quit
        if kappa < norm_r0^theta
            stopreason_str = 'reached target residual-kappa (linear)';
        else
            stopreason_str = 'reached target residual-theta (superlinear)';
        end
        break;
    end

    % Precondition the residual.
    z = getPrecon(problem, x, r, storedb, key);

    % Save the old z'*r.
    zold_rold = z_r;
    % Compute new z'*r.
    z_r = inner(z, r);

    % Compute new search direction.
    beta = z_r/zold_rold;
    mdelta = lincomb(1, z, beta, mdelta);

    % Since mdelta is passed to getHessian, which is the part of the code
    % we have least control over from here, we want to make sure mdelta is
    % a tangent vector up to numerical errors that should remain small.
    % For this reason, we re-project mdelta to the tangent space.
    % In limited tests, it was observed that it is a good idea to project
    % at every iteration rather than only every k iterations, the reason
    % being that loss of tangency can lead to more inner iterations being
    % run, which leads to an overall higher computational cost.
    mdelta = tangent(mdelta);

    % Update new P-norms and P-dots [CGT2000, eq. 7.5.6 & 7.5.7].
    e_Pd = beta*(e_Pd + alpha*d_Pd);
    d_Pd = z_r + beta*beta*d_Pd;

end  % of trs_tCG loop

if options.trscache
    store = storedb.get(key);
    store.store_iters = store_iters;
    if ~limitedbyTR
        % Store extra information since we did not exit because we were
        % limited by TR (model value increased or kappa/theta stopping
        % criterion satisfied)
        store_last = struct('numinner', j, 'stopreason_str', ...
            stopreason_str, 'eta', eta, 'Heta', Heta);
        memorytCG_MB = memorytCG_MB + getsize(store_last)/1024^2;

        if memorytCG_MB > options.memorytCG_warningtol
            warning('manopt:trs_tCG_cached:memory', ...
            [sprintf('trs_tCG_cached will cache %.2f [MB] for at least one iteration of trustregions until a step is accepted.', memorytCG_MB) ...
            'If memory is limited turn off caching by options.trscache = false.\n' ...
            'If more memory can be used without problem increase options.memorytCG_warningtol accordingly.\n' ...
            'To disable this warning: warning(''off'', ''manopt:trs_tCG_cached:memory'')']);
        end
        store.store_last = store_last;

        numstored = numstored + 1;
    end

    storedb.set(store, key);
end
stats = struct('numinner', j, 'hessvecevals', j, ...
                    'memorytCG_MB', memorytCG_MB);

if options.verbosity == 2
    printstr = sprintf('%9d   %9d   %9d   %s', j, j, numstored, ...
                stopreason_str);
elseif options.verbosity > 2
    printstr = sprintf('%9d   %9d   %9d   %9.2f   %s', j, j, numstored, ...
                memorytCG_MB, stopreason_str);
end

trsoutput.eta = eta;
trsoutput.Heta = Heta;
trsoutput.limitedbyTR = limitedbyTR;
trsoutput.printstr = printstr;
trsoutput.stats = stats;
end

author: NicolasBoumal

repo: manopt


text:
function [key,sdet]=posdef(A)
[nr,nc]=size(A); sdet=[];
for i=1:nr
   sdet=[sdet,det(A(1:i,1:i))];
end
key=1;
if any(sdet<=0), key=0; end


text:
% demo for creating an SPM M/EEG dataset from arbitrary data using
% conversion of simple Fieldtrip raw data struct.
% SPM8 internal format is quite complex but is transparent to the user
% as meeg object's methods take care of maintaining its consistency.
% The most straightforward way to convert arbitrary data that is available
% as an ASCII file or *.mat file with some variables to SPM8 is to create
% a quite simple Fieldtrip raw data struct and then use SPM's
% spm_eeg_ft2spm to convert this struct to SPM8 file. Missing information
% can then be supplemented using meeg methods and SPM functions.
% Fieldtrip raw struct must contain the following fields:
% .fsample - sampling rate (Hz)
% .trial - cell array of matrices with identical dimensions channels x time
% .time - cell array of time vectors (in sec), the same length as the
%         second dimension of the data. For SPM8 they must be identical.
% .label - cell array of strings, list of channel labels. Same length as
%         the first dimension of the data.
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston, Vladimir Litvak
% $Id: spm_lfp_txt2mat.m 2255 2008-09-30 15:36:59Z vladimir $

% load data (in this case a .mat file)
%--------------------------------------------------------------------------
load amy_hc_0501_R1

% define the output file name
%--------------------------------------------------------------------------
fname = 'LFP_example';

% define the sampling rate
%--------------------------------------------------------------------------
fsample = 1000;

% define epochs and create data array - this is specific for this data
% replace with your own code
%--------------------------------------------------------------------------
bins  = find(diff(data(:,3)) > 4);
bins  = bins([3,9]);
nbins = 1000*10;

% Define the channels of interest - in this case only 3,4 and 5
%--------------------------------------------------------------------------
Ic = [3 4 5];

% Create the Fieldtrip raw struct
%--------------------------------------------------------------------------
ftdata = [];

for i = 1:length(bins)
   It            = [1:nbins] + bins(i);
   ftdata.trial{i} = data(It, Ic)';
   ftdata.time{i} = [0:(nbins-1)]./fsample;
end

ftdata.fsample = fsample;
ftdata.label = colheaders(Ic);
ftdata.label = ftdata.label(:);

% Convert the ftdata struct to SPM M\EEG dataset
%--------------------------------------------------------------------------
D = spm_eeg_ft2spm(ftdata, fname);

% Examples of providing additional information in a script
% [] comes instead of an index vector and means that the command
% applies to all channels/all trials.
%--------------------------------------------------------------------------
D = type(D, 'single');             % Sets the dataset type
D = chantype(D, [], 'LFP');        % Sets the channel type
D = conditions(D, [], 'Sound 1');  % Sets the condition label

% save
%--------------------------------------------------------------------------
save(D);

author: spm


text:
function [W, a, train_err] = train_foe(images, basis, no_experts, a)
%TRAIN_FOE Trains a Field of Experts model on a collection of images
%
%   [W, a, train_err] = train_foe(images, basis, no_experts)
%   [W, a, train_err] = train_foe(images, basis, W, a)
%
% Trains a Field of Experts model on a set of (grayscale) images. The
% potential functions of the MRF are given by the energy function of a
% product of no_experts Student-t distributions (with parameters W and a).
% The potential functions are defined over a cliques that are image patches
% of patch_size x patch_size pixels.
%
%
% (C) Laurens van der Maaten, 2009
% Delft University of Technology

    % Initialize parameters
    eta_W = .01;
    eta_a = .1;
    max_iter = 1000;
    no_images = size(images, 1);
    size_im = sqrt(size(images, 2));
    batch_size = min(200, no_images);
    momentum = .9;
    train_err = zeros(max_iter, 1);
    addpath(genpath('netlab'));

    % Initialize model
    if exist('a', 'var')
        W = no_experts;
    else
        W = randn(size(basis, 1), no_experts);
        a = repmat(0.01, [no_experts 1]);
    end
    prev_dW = zeros(size(W));
    prev_da = zeros(size(a));

    % Perform training of the model
    for iter=1:max_iter

        % Create batches
        tic
        ind = randperm(no_images);
        disp(['Iteration ' num2str(iter) '...']);

        % Loop over batches
        for b=1:batch_size:no_images

            % Get batch and initialize positive and negative data
            batch = images(ind(b:b + batch_size - 1),:);
            pos_W = zeros(size(W)); pos_a = zeros(size(a));
            neg_W = zeros(size(W)); neg_a = zeros(size(a));

            % Loop over all training images to compute gradient terms
            for i=1:batch_size

                % Get image
                im = reshape(batch(i,:), [size_im size_im]);

%                 % Check gradients
%                 options = zeros(18, 1);
%                 options(9) = 1;
%                 hmc('foe_energy_w', W(1:end), options, 'foe_energy_grad_w', im, basis, a);
%                 hmc('foe_energy_a', a, options, 'foe_energy_grad_a', im, basis, W);
%                 hmc('foe_energy', im(1:end), options, 'foe_energy_grad_x', basis, W, a, [size_im size_im 1]);

                % Compute positive part of gradient
                pos_W = pos_W + reshape(foe_energy_grad_w(W(1:end), im, basis, a), size(W));
                pos_a = pos_a + foe_energy_grad_a(a, im, basis, W);

                % Draw samples from model using hybrid Monte Carlo
                options = zeros(18, 1);
                options(1)  = 0;                    % do not print diagnostics
                options(5)  = 0;                    % do not use momentum persistence
                options(7)  = 30;                   % number of leaps (= steps in leap-frog)
                options(9)  = 0;                    % do not check gradient
                options(14) = 1;                    % number of samples to return
                options(15) = 0;                    % number of samples to omit from start
                sample = hmc('foe_energy', batch(i,:), options, 'foe_energy_grad_x', basis, W, a, [size_im size_im 1]);

                % Compute negative CD part
                neg_W = neg_W + reshape(foe_energy_grad_w(W(1:end), sample, basis, a), size(W));
                neg_a = neg_a + foe_energy_grad_a(a, sample, basis, W);
            end

            % Perform the gradient updates
            prev_dW = eta_W * (momentum * prev_dW + (1 - momentum) * (1 / batch_size) * (pos_W - neg_W));
            prev_da = eta_a * (momentum * prev_da + (1 - momentum) * (1 / batch_size) * (pos_a - neg_a));
            W = W + prev_dW;
            a = exp(log(a) + a .* prev_da);
        end

        % Evaluate sum of energies of training data
        sum_E = 0;
        for i=1:no_images
            im = reshape(images(i,:), [size_im size_im]);
            sum_E = sum_E + foe_energy(im, basis, W, a);
        end
        disp(['     mean energy of training data is ' num2str(sum_E ./ no_images)]);
        train_err(iter) = sum_E ./ no_images;
        toc
    end

author: lbasek

repo: image-denoising-benchmark


text:
function output = calloptiooqp(interfacedata)

% Standard input interface
options = interfacedata.options;
F_struc = interfacedata.F_struc;
c       = full(interfacedata.c);
K       = interfacedata.K;
lb      = full(interfacedata.lb);
ub      = full(interfacedata.ub);
Q       = interfacedata.Q;

showprogress('Calling OOQP',options.showprogress);

blow = -inf(K.l,1);
if isempty(F_struc)
    A = sparse([]);
    bupp = [];
    blow = [];
    Aeq = sparse([]);
    beq = [];
else
    Aeq = -F_struc(1:K.f,2:end);
    beq = full(F_struc(1:K.f,1));
    A =-F_struc(K.f+1:end,2:end);
    bupp = full(F_struc(K.f+1:end,1));
end

if isempty(lb)
    lb = -inf(length(c),1);
    ub = inf(length(c),1);
end

opts = options.ooqp;

H = 2*sparse(triu(Q));
A = A';
Aeq = Aeq';
if K.f == 0
    Aeq = [];
    beq = [];
end
if K.l == 0
    A = [];
    blow = [];
    bupp = [];
end
if options.savedebug
    save ooqpdebug c H lb ub Aeq beq A blow bupp opts
end

solvertime = tic;
[x,fval,stat,iter] = ooqp(H, c, A, blow, bupp, Aeq,beq,lb,ub,opts);
solvertime = toc(solvertime);

% No duals
D_struc = [];

%   Return Status:
%       0 - optimal
%       1 - not finished
%       2 - maximum iterations exceeded
%       3 - infeasible
%       4 - ooqp error
switch stat
    case 0
        problem = 0;
    case 2
        problem = 3;
    case 3
        problem = 1;
    case 1
        problem = 11;
    case {1,4}
        problem = 9;
    otherwise
        problem = -1;
end
infostr = yalmiperror(problem,'OOQP');

% Save all data sent to solver?
if options.savesolverinput
    solverinput.A = A;
    solverinput.b = b;
    solverinput.f = c;
    solverinput.lb = lb;
    solverinput.ub = ub;
else
    solverinput = [];
end

% Save all data from the solver?
if options.savesolveroutput
    solveroutput.x = x;
    solveroutput.fmin = fmin;
    solveroutput.flag = flag;
    solveroutput.output=output;
else
    solveroutput = [];
end

% Standard interface
output = createOutputStructure(x,D_struc,[],problem,infostr,solverinput,solveroutput,solvertime);

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
% Script make_example_data
% Creates example data for uniQC-demos based on the ME-data in Heunis et al.
%
%  make_example_data
%
% NOTE: This script has to be run only once, in order to prepare some
%       example files for all uniqc-demos from a publicly available imaging
%       dataset.
%
% INSTRUCTIONS:
% 1. Please download the following data
%   (we only need subject sub-001 and sub-021) separately at
%    https://dataverse.nl/dataverse/rt-me-fmri.
%
%   from
%   Heunis, Stephan, 2020, "rt-me-fMRI: A task and resting state dataset for
%   real-time, multi-echo fMRI methods development and validation",
%   https://doi.org/10.34894/R1TNL8, DataverseNL, V1.
%
% 2. Unzip the downloaded subject folders
%
% 3. Update the paths in section "User Input" below, i.e.
%       a) where the downloaded, unzipped data resides (parent folder of sub-001)
%       b) where the generated uniqc example data files shall be stored
%          (for future use)
%
% 4. Run this script
%

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

clear; close all; clc;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% User Inputs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% file path to the Heunis-example data (we only need sub-001)
% PLEASE MODIFY TO YOUR PATH
sourceDataPath = 'path/to/downloaded/me_sh_data';

% where do want to store the example data for uniQC - this is also the path
% you will be asked for when running the examples
% PLEASE MODIFY TO YOUR PATH
exampleDataPath = '/path/to/save/uniQC-examples';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Load and plot source data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
m = MrImage(fullfile(sourceDataPath, 'sub-001', 'func', ...
    'sub-001_task-fingerTapping_echo-2_bold.nii'));
m.plot();

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Create small EPI data sets
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% nifti/rest/meanfmri.nii - just the mean of a times series (3D image)
meanfmri = m.mean();
meanfmri.parameters.save.path = fullfile(exampleDataPath, 'nifti', 'rest');
meanfmri.parameters.save.fileName = 'meanfmri.nii';
meanfmri.save();

% nifti/rest/fmri_short.nii - a short section of the fMRI time series to
% speed up processing in the demos
fmri_short = m.select('t', 1:15);
fmri_short.parameters.save.path = fullfile(exampleDataPath, 'nifti', 'rest');
fmri_short.parameters.save.fileName = 'fmri_short.nii';
fmri_short.save();

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Save anatomy
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
s = MrImage(fullfile(sourceDataPath, 'sub-001', 'anat', 'sub-001_T1w.nii'));
s.parameters.save.path = fullfile(exampleDataPath, 'nifti', 'rest');
s.parameters.save.fileName = 'struct.nii';
s.save();

% uniQC creates an extra mat file to keep track of any dimensions not saved
% in the nifti file, we delete it here to simulate the use case that we
% usually start only with a nifti
delete(fullfile(exampleDataPath, 'nifti', 'rest', '*.mat'));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Create 5D file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% nifti/5D/y_5d_deformation_field.nii
% this could for example be on of the deformation fields estimated during
% segmentation (the large samplingDistance makes this a quick and dirty
% version)
segmentInput = meanfmri.remove_dims.copyobj();
segmentInput.parameters.save.keepCreatedFiles = 'processed';
segmentInput.parameters.save.fileName = '5D.nii';
segmentInput.parameters.save.path = fullfile(exampleDataPath, 'nifti', '5D');

[~, ~, deformationFields] = segmentInput.remove_dims.segment(...
    'samplingDistance', 20, 'tissueTypes', {'GM'});
movefile(deformationFields{1}.get_filename, ...
    fullfile(exampleDataPath, 'nifti', '5D', 'y_5d_deformation_field.nii'));
delete(fullfile(exampleDataPath, 'nifti', '5D', 'bias*'));
delete(fullfile(exampleDataPath, 'nifti', '5D', 'tissue*'));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Create split file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% now load the complete multi-echo data set
meFileNames = dir(fullfile(sourceDataPath, 'sub-001', 'func', 'sub-001_task-fingerTapping_echo-*_bold.nii'));
for nFiles = 1:numel(meFileNames)
    % load files into temporary array
    meTmp{nFiles} = MrImage(fullfile(sourceDataPath, 'sub-001', 'func',...
        meFileNames(nFiles).name));
    % get echo time from json data
    metaInfo = jsondecode(fileread(fullfile(sourceDataPath, 'sub-001', 'func', ...
        regexprep(meFileNames(nFiles).name, '.nii', '.json'))));
    meTmp{nFiles}.dimInfo.add_dims('TE', ...
        'samplingPoints', metaInfo.EchoTime*1000, 'units', 'ms');
end

% combine into 5D image
me = meTmp{1}.combine(meTmp);

% save using the uniQC convention
% nift/split/subset/fmri_***.nii
meSplit = me.select('t', [1, 7, 8], 'TE', [1, 2]);
meSplit.parameters.save.path = fullfile(exampleDataPath, 'nifti', 'split', 'subset');
meSplit.parameters.save.fileName = 'fmri.nii';
meSplit.save('splitDims', {'t', 'TE'});
delete(fullfile(exampleDataPath,'nifti', 'split', 'subset', '*.mat'));

% nifti/data_multi_echo/multi_echo.nii
me.parameters.save.path = fullfile(exampleDataPath, 'nifti', 'data_multi_echo');
me.parameters.save.fileName = 'multi_echo.nii';
me.save();
delete(fullfile(exampleDataPath, 'nifti', 'data_multi_echo', '*.mat'));

% save without uniQC conventions
% /nifti/split/split_residual_images/Res.nii'
for n = 1:5
    currentM = me.select('t', n, 'TE', 1).remove_dims();
    currentM.parameters.save.path = fullfile(exampleDataPath, 'nifti', 'split', 'split_residual_images');
    currentM.parameters.save.fileName = ['Res_', num2str(n, '%04.0f'), '.nii'];
    currentM.save();
end
delete(fullfile(exampleDataPath, 'nifti', 'split', 'split_residual_images', '*.mat'));

% save with uniQC dimInfo
% /nifti/split/with_dim_info/MrImage.nii'
meSplitWith = me.select('t', [1, 7, 8]);
meSplitWith.parameters.save.path = fullfile(exampleDataPath, 'nifti', 'split', 'with_dim_info');
meSplitWith.parameters.save.fileName = 'MrImage.nii';
meSplitWith.save();

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Create data for first level
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% /nifti/data_first_level/
% structural image
s.parameters.save.path = fullfile(exampleDataPath, 'nifti', 'data_first_level');
s.parameters.save.fileName = 'anatomy.nii';
s.save();

% third echo as functional image
fl = me.select('TE', 3).remove_dims();
fl.parameters.save.path = fullfile(exampleDataPath, 'nifti', 'data_first_level');
fl.parameters.save.fileName = 'single_echo.nii';
fl.save();
delete(fullfile(exampleDataPath, 'nifti', 'data_first_level', '*.mat'));

author: translationalneuromodeling

repo: tapas


text:
function varargout = efficency(X,P)
% function [con_eff,hrf_eff,X,hrfX,P] = efficency(X,PARAMS)
% a useful summary function for a design, incorporating tests of contrasts and of HRF shape estimation
%
% X
%   either a design vector (condition function) with integers indicating onsets of conditions
%   or an already-constructed design matrix
%
% PARAMS structure contains inputs
% Optional fields:
%   contrasts   required for testing contrast efficiency; see calcEfficiency.m
%   Vi          intrinsic autocorrelation matrix; empty assumes identity matrix (independence)
%   S           smoothing matrix; empty for no smoothing
%   HRF         canonical hemodynamic response function; empty assumes SPM's canonical HRF
%               should be in resolution of .1 s per element
%   nonlint     saturation threshold (ceiling) for predictors; simple nonlinearity model
%   HRFtime     number of s to estimate HRF shape for, in FIR model.
%   delta       You must enter this if you want HRF shape estimate eff and you enter a model matrix X
%               Delta is n x m, with columns for conditions, made of ones and zeros, in TR-length time bins
%
% Required fields for condition function (vector) input:
%   ISI         sampling resolution of design vector, as time between elements in s
%   TR          sampling resolution of final design matrix (volumes acquired), in s
%
% Functions Called: (from OptimizeDesign GA toolbox)
% designvector2model.m
% tor_make_deconv_mtx2.m
% calcEfficiency.m

if ~isfield(P,'contrasts'), P.contrasts = eye(size(X,2));, end
if ~isfield(P,'HRFtime'), P.HRFtime = 30;, end
if isfield(P,'delta'), delta = P.delta;, end
if isfield(P,'TR'), TR = P.TR;, end
if isfield(P,'ISI'), ISI = P.ISI;, end

% -------------------------------------------------------------------------------------------------
% * Set up design matrix
% -------------------------------------------------------------------------------------------------

% For condition function
if any(size(X) == 1)

    if ~isfield(P,'ISI'), error('ISI is required field in PARAMS input'), end
    if ~isfield(P,'TR'), error('TR is required field in PARAMS input'), end
    if ~isfield(P,'S'), P.S = [];, end
    if ~isfield(P,'HRF'), P.HRF = spm_hrf(.1);, HRF = HRF ./ max(HRF);, end
    if ~isfield(P,'nonlint'), P.nonlint = [];, end

    numsamps = ceil(length(X)*P.ISI/P.TR);

    if nargout > 1
        delta = [];
        for i = 1:max(double(X))
            delta(:,i) = (X == i);
        end
    end

    X = designvector2model(X,P.ISI,P.HRF,P.TR,numsamps,P.nonlint,P.S);

    if ~isfield(P,'Vi'), P.Vi = eye(size(X,1));, end

% For a design matrix
else
    if ~isfield(P,'S'), P.S = [];, end
    if ~isfield(P,'Vi'), P.Vi = eye(size(X,1));, end

end

if ~isempty(P.S), svi = P.S * P.Vi;, else, svi = P.Vi;, end

if ~isempty(P.contrasts)
        % -------------------------------------------------------------------------------------------------
		% * efficiency
		% -------------------------------------------------------------------------------------------------

        	xtxitx = pinv(X);                                       		% inv(X'S'SX)*(SX)'; pseudoinv of (S*X)
            contrastweights = ones(1,size(P.contrasts,1));
			[dummy,varargout{1}] = calcEfficiency(contrastweights,P.contrasts,xtxitx,svi);

else
        varargout{1} = [];
end

if isfield(P,'ISI') & isfield(P,'TR') & nargout > 1 & exist('delta') == 1

        % -------------------------------------------------------------------------------------------------
		% * HRF shape estimation efficiency
		% -------------------------------------------------------------------------------------------------

			[X2] = tor_make_deconv_mtx3(delta,round(P.HRFtime / P.TR),P.TR / P.ISI);
            if ~isempty(P.S), X2 = P.S * X2;,end

        	xtxitx = pinv(X2);                                       		% inv(X'S'SX)*(SX)'; pseudoinv of (S*X)
			[dummy,varargout{2}] = calcEfficiency([],[],xtxitx,svi);

elseif nargout > 1
            disp(['Missing ISI or TR field, or already-constructed model entered: HRF efficiency not calculated.'])
            varargout{2} = [];
end

if nargout > 2
    varargout{3} = X; varargout{4} = X2;
    P.hrfsamples = round(P.HRFtime / P.TR);
    P.hrftimeres = P.TR / P.ISI;
    varargout{5} = P;
end

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
classdef SettingsLevelSetRectangleInclusion < SettingsLevelSetCreator

    properties (Access = public)
        widthH
        widthV
    end

    methods (Access = public)

        function obj = SettingsLevelSetRectangleInclusion(varargin)
            obj.loadParams('paramsLevelSetCreator_RectangleInclusion.json')
            if nargin == 1
                obj.loadParams(varargin{1})
            end
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function line2 = normalizeLine3d(line)
%NORMALIZELINE3D Normalizes the direction vector of a 3D line.
%
%   LINE2 = normalizeVector3d(LINE);
%   Returns the normalization of the direction vector of the line, such
%   that ||LINE2(4:6)|| = 1.
%
%   See also
%     normalizePlane, normalizeVector3d
%

% ------
% Author: oqilipo
% E-mail: N/A
% Created: 2020-03-13
% Copyright 2020-2022

isLine3d = @(x) validateattributes(x,{'numeric'},...
    {'nonempty','nonnan','real','finite','size',[nan,6]});

% Check if the line is valid
p=inputParser;
addRequired(p,'line',isLine3d)
parse(p,line)

line2 = line;
line2(:,4:6) = normalizeVector3d(line2(:,4:6));

end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
%DXform Distance transform navigation class
%
% A concrete subclass of the abstract Navigation class that implements the distance
% transform navigation algorithm which computes minimum distance paths.
%
% Methods::
%  DXform       Constructor
%  plan         Compute the cost map given a goal and map
%  query        Find a path
%  plot         Display the distance function and obstacle map
%  plot3d       Display the distance function as a surface
%  display      Print the parameters in human readable form
%  char         Convert to string
%
% Properties (read only)::
%  distancemap   Distance from each point to the goal.
%  metric        The distance metric, can be 'euclidean' (default) or 'cityblock'
%
% Example::
%
%        load map1           % load map
%        goal = [50,30];     % goal point
%        start = [20, 10];   % start point
%        dx = DXform(map);   % create navigation object
%        dx.plan(goal)       % create plan for specified goal
%        dx.query(start)     % animate path from this start location
%
% Notes::
% - Obstacles are represented by NaN in the distancemap.
% - The value of each element in the distancemap is the shortest distance from the
%   corresponding point in the map to the current goal.
%
% References::
% -  Robotics, Vision & Control, Sec 5.2.1,
%    Peter Corke, Springer, 2011.
%
% See also Navigation, Dstar, PRM, distancexform.

% Copyright (C) 1993-2017, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

classdef DXform < Navigation

    properties
        metric;     % distance metric
        distancemap;   % distance transform results
    end

    methods

        function dx = DXform(world, varargin)
            %DXform.DXform Distance transform constructor
            %
            % DX = DXform(MAP, OPTIONS) is a distance transform navigation object,
            % and MAP is an occupancy grid, a representation of a planar
            % world as a matrix whose elements are 0 (free space) or 1
            % (occupied).
            %
            % Options::
            % 'goal',G      Specify the goal point (2x1)
            % 'metric',M    Specify the distance metric as 'euclidean' (default)
            %               or 'cityblock'.
            % 'inflate',K   Inflate all obstacles by K cells.
            %
            % Other options are supported by the Navigation superclass.
            %
            % See also Navigation.Navigation.

            % TODO NEEDS PROPER ARG HANDLER

            % invoke the superclass constructor
            dx = dx@Navigation(world, varargin{:});

            opt.metric = {'euclidean', 'cityblock'};
            [opt,args] = tb_optparse(opt, varargin);
            dx.metric = opt.metric;
        end

        function s = char(dx)
            %DXform.char Convert to string
            %
            % DX.char() is a string representing the state of the object in
            % human-readable form.
            %
            % See also DXform.display, Navigation.char

            % most of the work is done by the superclass
            s = char@Navigation(dx);

            % dxform specific stuff
            s = char(s, sprintf('  distance metric: %s', dx.metric));
            if ~isempty(dx.distancemap)
                s = char(s, sprintf('  distancemap: computed:'));
            else
                s = char(s, sprintf('  distancemap: empty:'));
            end
        end

        % invoked by superclass on a change of goal, mark the distancemap
        % as invalid
        function goal_change(dx, goal)

            dx.distancemap = [];
            if dx.verbose
                disp('Goal changed -> distancemap cleared');
            end
        end

        function plan(dx, varargin)
            %DXform.plan Plan path to goal
            %
            % DX.plan(GOAL, OPTIONS) plans a path to the goal given to the constructor,
            % updates the internal distancemap where the value of each element is the
            % minimum distance from the corresponding point to the goal.
            %
            % DX.plan(GOAL, OPTIONS) as above but goal is specified explicitly
            %
            % Options::
            % 'animate'    Plot the distance transform as it evolves
            %
            % Notes::
            % - This may take many seconds.
            %
            % See also Navigation.path.

            opt.animate = false;

            [opt,args] = tb_optparse(opt, varargin);

            if opt.animate
                show = 0.05;
            else
                show = 0;
            end

            if ~isempty(args) && isvec(args{1},2)
                dx.setgoal(args{1});
            end

            assert(~isempty(dx.goal), 'RTB:DXform:plan', 'no goal specified here or in constructor');

            dx.distancemap = distancexform(dx.occgridnav, dx.goal, dx.metric, show);
        end

        function plot(dx, varargin)
            %DXform.plot Visualize navigation environment
            %
            % DX.plot(OPTIONS) displays the occupancy grid and the goal distance
            % in a new figure.  The goal distance is shown by intensity which
            % increases with distance from the goal.  Obstacles are overlaid
            % and shown in red.
            %
            % DX.plot(P, OPTIONS) as above but also overlays a path given by the set
            % of points P (Mx2).
            %
            % Notes::
            % - See Navigation.plot for options.
            %
            % See also Navigation.plot.

            plot@Navigation(dx, varargin{:}, 'distance', dx.distancemap);

        end

        function n = next(dx, robot)
            if isempty(dx.distancemap)
                error('No distancemap computed, you need to plan');
            end

            % list of all possible directions to move from current cell
            directions = [
                -1 -1
                0 -1
                1 -1
                -1 0
                0 0
                1 0
                -1 1
                0 1
                1 1];

            x = robot(1); y = robot(2);

            % find the neighbouring cell that has the smallest distance
            mindist = Inf;
            mindir = [];
            for d=directions'
                % use exceptions to catch attempt to move outside the map
                try
                    if dx.distancemap(y+d(1), x+d(2)) < mindist
                        mindir = d;
                        mindist = dx.distancemap(y+d(1), x+d(2));
                    end
                catch
                end
            end

            x = x + mindir(2);
            y = y + mindir(1);

            if all([x;y] == dx.goal)
                n = [];     % indicate we are at the goal
            else
                n = [x; y];  % else return the next closest point to the goal
            end
        end % next

        function plot3d(dx, p, varargin)
            %DXform.plot3d 3D costmap view
            %
            % DX.plot3d() displays the distance function as a 3D surface with
            % distance from goal as the vertical axis.  Obstacles are "cut out"
            % from the surface.
            %
            % DX.plot3d(P) as above but also overlays a path given by the set
            % of points P (Mx2).
            %
            % DX.plot3d(P, LS) as above but plot the line with the MATLAB linestyle LS.
            %
            % See also Navigation.plot.
            surf(dx.distancemap);
            shading interp

            if nargin > 1
                % plot path if provided
                k = sub2ind(size(dx.distancemap), p(:,2), p(:,1));
                height = dx.distancemap(k);
                hold on
                if isempty(varargin)
                    varargin{1} = 'k.';
                end
                plot3(p(:,1), p(:,2), height, varargin{:})
                hold off
            end
        end
    end % methods
end % classdef

author: petercorke

repo: robotics-toolbox-matlab


text:
function rmsData=RmsFin(rmsWindow,rmsStride,foldedData)
%%
rmsData=foldedData;
for i =1:length(foldedData)
    %data is trial*channel*time
    data=permute(foldedData(i).x,[3 2 1]);
    [nTrials,nChannel,nTime]=size(data);
    for ii=1:rmsStride:nTime-rmsWindow
        window=zeros(nTrials,nChannel,rmsWindow);
        window=data(:,:,ii:ii+rmsWindow);
        %% rms
        windowSquare=window.*window;
        windowRms=sqrt(sum(windowSquare,3)/rmsWindow);
        if ii==1
            windowFin=windowRms;
            time=1;
        else
            time=cat(2,time,ii);
            windowFin=cat(3,windowFin,windowRms);
        end
    end
    windowFin=permute(windowFin,[3 2 1]);
    rmsData(i).x=windowFin;
    rmsData(i).t=time;
end


text:
function c = i4_division ( a, b )

%*****************************************************************************80
%
%% I4_DIVISION returns the result of integer division.
%
%  Discussion:
%
%    This routine computes C = A / B, where the result is rounded to the
%    integer value nearest 0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer A, B, the number to be divided,
%    and the divisor.
%
%    Output, integer C, the rounded result of the division.
%
  if ( a * b < 0.0 )
    s = -1;
  else
    s = +1;
  end

  a = abs ( a );
  b = abs ( b );
  c = s * floor ( a / b );

  return
end


text:
function writedouble(filename,x)
% write to file

fid=fopen(filename,'w');
fwrite(fid,x,'double');
fclose(fid);

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function [c1,c2]= binaryfit(phi,U,epsilon)
% compute c1 c2 for optimal binary fitting
% input:
%    U: input image
%    phi: level set function
%    epsilon: parameter for computing smooth Heaviside and dirac function
% output:
%    c1: a constant to fit the image U in the region phi>0
%    c2: a constant to fit the image U in the region phi<0

H = Heaviside(phi,epsilon); % compute the Heaveside function values

a = H .* U;
numer_1 = sum(a(:));
denom_1 = sum(H(:));
c1 = numer_1 / denom_1;

b = (1-H) .* U;
numer_2 = sum(b(:));
c = 1-H;
denom_2 = sum(c(:));
c2 = numer_2 / denom_2;

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function isInside = bst_intriangle(A, B, C, P)
% BST_INTRIANGLE: Test if the orthogonal projection of points P on a triangle (A,B,C) are inside the triangle.

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2013

% Compute vectors
v0 = repmat(C - A, size(P,1), 1);
v1 = repmat(B - A, size(P,1), 1);
v2 = bst_bsxfun(@minus, P, A);

% Compute dot products
dot00 = sum(v0.^2, 2);
dot01 = sum(v0.*v1, 2);
dot02 = sum(v0.*v2, 2);
dot11 = sum(v1.^2, 2);
dot12 = sum(v1.*v2, 2);

% Compute barycentric coordinates
invDenom = 1 ./ (dot00 .* dot11 - dot01 .* dot01);
u = (dot11 .* dot02 - dot01 .* dot12) .* invDenom;
v = (dot00 .* dot12 - dot01 .* dot02) .* invDenom;

% Check if point is in triangle
isInside = (u >= 0) & (v >= 0) & (u + v < 1);

author: brainstorm-tools


text:
function [paramsOut] = globalTrendGUI(hObject, callbackdata, inputData)%#ok<INUSL>
    theTitle ='Global trend parameters';
    defaultStruct = inputData.userData.globalTrend;

    while(true)
        mainFigure = findobj('Type', 'Figure', '-and', 'Name', inputData.name);
        userdata = get(mainFigure, 'UserData');
        if isempty(userdata) || ~isfield(userdata, 'globalTrend')
            paramsOut = struct();
        else
            paramsOut = userdata.globalTrend;
        end
        [defaultStruct, errors] = checkStructureDefaults(paramsOut, defaultStruct);

        if ~isempty(errors)
            warning('globalTrendGUI:bad parameters', getMessageString(errors)); %#ok<CTPCT>
        end

        % Creates structure for
        fNamesDefault = fieldnames(defaultStruct);
        for k = 1:length(fNamesDefault)
            textColorStruct.(fNamesDefault{k}) = 'k';
        end
        if defaultStruct.doLocal.value
            checkValue = 1;
        else
            checkValue = 0;
        end

        %% starts the while loop, sets up the uilist and creates the GUI
        closeOpenWindows(theTitle);
        geometry = {[1,4], 1, [3,1,3,1],[3,1,4]};
        geomvert = [];
        uilist={{'style', 'text', 'string', 'Global trend channels', ...
            'TooltipString', defaultStruct.globalTrendChannels.description}...
            {'style', 'edit', 'string', ...
            num2str(defaultStruct.globalTrendChannels.value), ...
            'tag', 'globalTrendChannels', 'ForegroundColor', ...
            textColorStruct.globalTrendChannels}...
            {'style', 'text', 'string', ''}...
            {'style', 'text', 'string', 'Do local', 'TooltipString', ...
            defaultStruct.doLocal.description}...
            {'style', 'checkbox', 'Value', checkValue, ...
            'tag', 'doLocal', 'ForegroundColor', textColorStruct.doLocal}...
            {'style', 'text', 'string', 'Local cutoff', ...
            'TooltipString', defaultStruct.localCutoff.description}...
            {'style', 'edit', 'string', ...
            num2str(defaultStruct.localCutoff.value), ...
            'tag', 'localCutoff', 'ForegroundColor', textColorStruct.localCutoff}...
            {'style', 'text', 'string', 'Local step size', ...
            'TooltipString', defaultStruct.localStepSize.description}...
            {'style', 'edit', 'string', ...
            num2str(defaultStruct.localStepSize.value), ...
            'tag', 'localStepSize', 'ForegroundColor', ...
            textColorStruct.localStepSize}...
            {'style', 'text', 'string', ''}};
        [~, ~, ~, paramsOut] = inputgui('geometry', geometry, ...
            'geomvert', geomvert, 'uilist', uilist, 'title', theTitle, ...
            'helpcom', 'pophelp(''pop_prepPipeline'')');
        if(isempty(paramsOut))
            break;
        end

        [paramsOut, typeErrors, fNamesErrors] = ...
            changeType(paramsOut, defaultStruct);
        %paramsOut.detrendType = typeMenuString(paramsOut.detrendType);
        mainFigure = findobj('Type', 'Figure', '-and', 'Name', inputData.name);
        userdata = get(mainFigure, 'UserData');
        userdata.globalTrend = paramsOut;
        set(mainFigure, 'UserData', userdata);
        if isempty(typeErrors)
            break;
        end
        textColorStruct = highlightErrors(fNamesErrors, ...
            fNamesDefault, textColorStruct);

        displayErrors(typeErrors); % Displays the errors and restarts GUI
    end
end

author: VisLab

repo: EEG-Clean-Tools


text:
%-------------------------------------------------------------------------------------------------------------------%
%
% IB2d is an Immersed Boundary Code (IB) for solving fully coupled non-linear
% 	fluid-structure interaction models. This version of the code is based off of
%	Peskin's Immersed Boundary Method Paper in Acta Numerica, 2002.
%
% Author: Nicholas A. Battista
% Email:  nickabattista@gmail.com
% Date Created: May 27th, 2015
% Institution: UNC-CH
%
% This code is capable of creating Lagrangian Structures using:
% 	1. Springs
% 	2. Beams (*torsional springs, non-invariant beams*)
% 	3. Target Points
%	4. Muscle-Model (combined Force-Length-Velocity model, "Hill+(Length-Tension)")
%
% One is able to update those Lagrangian Structure parameters, e.g., spring constants, resting lengths, etc
%
% There are a number of built in Examples, mostly used for teaching purposes.
%
% If you would like us to add a specific muscle model, please let Nick (nickabattista@gmail.com) know.
%
%--------------------------------------------------------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: updates the spring attributes!
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function springs_info = update_Springs(dt,current_time,xLag,yLag,springs_info)

%springs_info: col 1: starting spring pt (by lag. discretization)
%              col 2: ending spring pt. (by lag. discretization)
%              col 3: spring stiffness
%              col 4: spring resting lengths

%RL = springs_info(:,4); % resting-length vector

% Contraction Frequency
freq = 1.0;

% CHANGE RESTING LENGTH BTWN SIDES OF JELLYFISH BELL
springs_info(159:end,4) = abs( cos(freq*current_time*pi) );

%NOTE: not 2*pi*ft b/c of abs()

author: nickabattista


text:
function [BCMatrix, BCRHS] = boundaryConditionCylindrical3D(BC)
% It creates the matrix of coefficient based on the BC structureprovided
% by the user. It also generates the right hand side vector of the linear
% system of equations
%
% SYNOPSIS:
%       [BCMatrix, BCRHS] = boundaryConditionCylindrical3D(BC)
%
% PARAMETERS:
%
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

% Note: I use a for loop here fr more readability of the code!

% extract data from the mesh structure
Nxyz = BC.domain.dims;
Nx = Nxyz(1); Ntetta = Nxyz(2); Nz = Nxyz(3);
G=reshape(1:(Nx+2)*(Ntetta+2)*(Nz+2), Nx+2, Ntetta+2, Nz+2);
dx_1 = BC.domain.cellsize.x(1);
dx_end = BC.domain.cellsize.x(end);
dtetta_1 = BC.domain.cellsize.y(1);
dtetta_end = BC.domain.cellsize.y(end);
dz_1 = BC.domain.cellsize.z(1);
dz_end = BC.domain.cellsize.z(end);
rp = repmat(BC.domain.cellcenters.x, 1, Nz);

% number of boundary nodes (axact number is 2[(m+1)(n+1)*(n+1)*(p+1)+(m+1)*p+1]:
nb = 8*((Nx+1)*(Ntetta+1)+(Nx+1)*(Nz+1)+(Ntetta+1)*(Nz+1));

% define the vectors to be used for the creation of the sparse matrix
ii = zeros(nb,1);
jj = zeros(nb,1);
s = zeros(nb,1);

% define the RHS column vector
BCRHS = zeros((Nx+2)*(Ntetta+2)*(Nz+2), 1);

% assign value to the corner nodes (useless cells)
q = 1:8;
ii(q) = BC.domain.corners; jj(q) = BC.domain.corners;
s(q) = 1; BCRHS(BC.domain.corners) = 0;

% assign values to the edges (useless cells)
q = q(end)+(1:length(BC.domain.edges));
ii(q) = BC.domain.edges; jj(q) = BC.domain.edges;
s(q) = 1; BCRHS(BC.domain.edges) = 0;

% Assign values to the boundary condition matrix and the RHS vector based
% on the BC structure
if (BC.top.periodic ==0) && (BC.bottom.periodic == 0)
    % top boundary
    j=Ntetta+2;
    i=2:Nx+1;
    k=2:Nz+1;
    q = q(end)+(1:Nx*Nz);
    ii(q) = G(i,j,k);  jj(q) = G(i,j,k);  s(q) = BC.top.b/2 + BC.top.a./(dtetta_end*rp);
    q = q(end)+(1:Nx*Nz);
    ii(q) = G(i,j,k);  jj(q) = G(i,j-1,k); s(q) = BC.top.b/2 - BC.top.a./(dtetta_end*rp);
    BCRHS(G(i,j,k)) = BC.top.c;

    % Bottom boundary
    j=1;
    i=2:Nx+1;
    k=2:Nz+1;
    q = q(end)+(1:Nx*Nz);
    ii(q) = G(i,j,k);  jj(q) = G(i,j+1,k);  s(q) = -(BC.bottom.b/2 + BC.bottom.a./(dtetta_1*rp)); % consider the reverse direction of normal
    q = q(end)+(1:Nx*Nz);
    ii(q) = G(i,j,k);  jj(q) = G(i,j,k); s(q) = -(BC.bottom.b/2 - BC.bottom.a./(dtetta_1*rp)); % consider the reverse direction of normal
    BCRHS(G(i,j,k)) = -(BC.bottom.c);
elseif (BC.top.periodic ==1) || (BC.bottom.periodic == 1) % periodic
    % top boundary
    j=Ntetta+2;
    i=2:Nx+1;
    k=2:Nz+1;
    q = q(end)+(1:Nx*Nz);
    ii(q) = G(i,j,k);  jj(q) = G(i,j,k);  s(q) = 1;
    q = q(end)+(1:Nx*Nz);
    ii(q) = G(i,j,k);  jj(q) = G(i,j-1,k);  s(q) = -1;
    q = q(end)+(1:Nx*Nz);
    ii(q) = G(i,j,k);  jj(q) = G(i,1,k); s(q) = dtetta_end/dtetta_1;
    q = q(end)+(1:Nx*Nz);
    ii(q) = G(i,j,k);  jj(q) = G(i,2,k); s(q) = -dtetta_end/dtetta_1;
    BCRHS(G(i,j,k)) = 0;

    % Bottom boundary
    j=1;
    i=2:Nx+1;
    k=2:Nz+1;
    q = q(end)+(1:Nx*Nz);
    ii(q) = G(i,j,k);  jj(q) = G(i,j,k);  s(q) = 1;
    q = q(end)+(1:Nx*Nz);
    ii(q) = G(i,j,k);  jj(q) = G(i,j+1,k);  s(q) = 1;
    q = q(end)+(1:Nx*Nz);
    ii(q) = G(i,j,k);  jj(q) = G(i,Ntetta+1,k); s(q) = -1;
    q = q(end)+(1:Nx*Nz);
    ii(q) = G(i,j,k);  jj(q) = G(i,Ntetta+2,k); s(q) = -1;
    BCRHS(G(i,j,k)) = 0;
end

if (BC.right.periodic == 0) && (BC.left.periodic == 0)
    % Right boundary
    i=Nx+2;
    j=2:Ntetta+1;
    k=2:Nz+1;
    q = q(end)+(1:Ntetta*Nz);
    ii(q) = G(i,j,k);  jj(q) = G(i,j,k);  s(q) = BC.right.b/2 + BC.right.a/dx_end;
    q = q(end)+(1:Ntetta*Nz);
    ii(q) = G(i,j,k);  jj(q) = G(i-1,j,k); s(q) = BC.right.b/2 - BC.right.a/dx_end;
    BCRHS(G(i,j,k)) = BC.right.c;

    % Left boundary
    i = 1;
    j=2:Ntetta+1;
    k=2:Nz+1;
    q = q(end)+(1:Ntetta*Nz);
    ii(q) = G(i,j,k);  jj(q) = G(i+1,j,k);  s(q) = -(BC.left.b/2 + BC.left.a/dx_1); % consider the reverse direction of normal
    q = q(end)+(1:Ntetta*Nz);
    ii(q) = G(i,j,k);  jj(q) = G(i,j,k); s(q) = -(BC.left.b/2 - BC.left.a/dx_1); % consider the reverse direction of normal
    BCRHS(G(i,j,k)) = -(BC.left.c);
elseif (BC.right.periodic == 1) || (BC.left.periodic == 1) % periodic
    % Right boundary
    i=Nx+2;
    j=2:Ntetta+1;
    k=2:Nz+1;
    q = q(end)+(1:Ntetta*Nz);
    ii(q) = G(i,j,k);  jj(q) = G(i,j,k);  s(q) = 1;
    q = q(end)+(1:Ntetta*Nz);
    ii(q) = G(i,j,k);  jj(q) = G(i-1,j,k);  s(q) = -1;
    q = q(end)+(1:Ntetta*Nz);
    ii(q) = G(i,j,k);  jj(q) = G(1,j,k); s(q) = dx_end/dx_1;
    q = q(end)+(1:Ntetta*Nz);
    ii(q) = G(i,j,k);  jj(q) = G(2,j,k); s(q) = -dx_end/dx_1;
    BCRHS(G(i,j,k)) = 0;

    % Left boundary
    i = 1;
    j=2:Ntetta+1;
    k=2:Nz+1;
    q = q(end)+(1:Ntetta*Nz);
    ii(q) = G(i,j,k);  jj(q) = G(i,j,k);  s(q) = 1;
    q = q(end)+(1:Ntetta*Nz);
    ii(q) = G(i,j,k);  jj(q) = G(i+1,j,k);  s(q) = 1;
    q = q(end)+(1:Ntetta*Nz);
    ii(q) = G(i,j,k);  jj(q) = G(Nx+1,j,k); s(q) = -1;
    q = q(end)+(1:Ntetta*Nz);
    ii(q) = G(i,j,k);  jj(q) = G(Nx+2,j,k); s(q) = -1;
    BCRHS(G(i,j,k)) = 0;
end

if (BC.front.periodic == 0) && (BC.back.periodic == 0)
    % Front boundary
    k=Nz+2;
    i = 2:Nx+1;
    j=2:Ntetta+1;
    q = q(end)+(1:Nx*Ntetta);
    ii(q) = G(i,j,k);  jj(q) = G(i,j,k);  s(q) = BC.front.b/2 + BC.front.a/dz_end;
    q = q(end)+(1:Nx*Ntetta);
    ii(q) = G(i,j,k);  jj(q) = G(i,j,k-1); s(q) = BC.front.b/2 - BC.front.a/dz_end;
    BCRHS(G(i,j,k)) = BC.front.c;

    % Back boundary
    k=1;
    i = 2:Nx+1;
    j=2:Ntetta+1;
    q = q(end)+(1:Nx*Ntetta);
    ii(q) = G(i,j,k);  jj(q) = G(i,j,k+1);  s(q) = -(BC.back.b/2 + BC.back.a/dz_1); % consider the reverse direction of normal
    q = q(end)+(1:Nx*Ntetta);
    ii(q) = G(i,j,k);  jj(q) = G(i,j,k); s(q) = -(BC.back.b/2 - BC.back.a/dz_1); % consider the reverse direction of normal
    BCRHS(G(i,j,k)) = -(BC.back.c);
elseif (BC.front.periodic == 1) || (BC.back.periodic == 1) % periodic
    % Front boundary
    k=Nz+2;
    i = 2:Nx+1;
    j=2:Ntetta+1;
    q = q(end)+(1:Nx*Ntetta);
    ii(q) = G(i,j,k);  jj(q) = G(i,j,k);  s(q) = 1;
    q = q(end)+(1:Nx*Ntetta);
    ii(q) = G(i,j,k);  jj(q) = G(i,j,k-1);  s(q) = -1;
    q = q(end)+(1:Nx*Ntetta);
    ii(q) = G(i,j,k);  jj(q) = G(i,j,1); s(q) = dz_end/dz_1;
    q = q(end)+(1:Nx*Ntetta);
    ii(q) = G(i,j,k);  jj(q) = G(i,j,2); s(q) = -dz_end/dz_1;
    BCRHS(G(i,j,k)) = 0;

    % Back boundary
    k=1;
    i = 2:Nx+1;
    j=2:Ntetta+1;
    q = q(end)+(1:Nx*Ntetta);
    ii(q) = G(i,j,k);  jj(q) = G(i,j,k);  s(q) = 1;
    q = q(end)+(1:Nx*Ntetta);
    ii(q) = G(i,j,k);  jj(q) = G(i,j,k+1);  s(q) = 1;
    q = q(end)+(1:Nx*Ntetta);
    ii(q) = G(i,j,k);  jj(q) = G(i,j,Nz+1); s(q) = -1;
    q = q(end)+(1:Nx*Ntetta);
    ii(q) = G(i,j,k);  jj(q) = G(i,j,Nz+2); s(q) = -1;
    BCRHS(G(i,j,k)) = 0;
end

% Build the sparse matrix of the boundary conditions
BCMatrix = sparse(ii(1:q(end)), jj(1:q(end)), s(1:q(end)), ...
    (Nx+2)*(Ntetta+2)*(Nz+2), (Nx+2)*(Ntetta+2)*(Nz+2));

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool