text:
%%%%%%%VERSION 2
%%ANOTHER DESCRIBTION OF GABOR FILTER

%The Gabor filter is basically a Gaussian (with variances sx and sy along x and y-axes respectively)
%modulated by a complex sinusoid (with centre frequencies U and V along x and y-axes respectively)
%described by the following equation
%%
%                            -1     x' ^     y'  ^
%%% G(x,y,theta,f) =  exp ([----{(----) 2+(----) 2}])*cos(2*pi*f*x');
%                             2    sx'       sy'
%%% x' = x*cos(theta)+y*sin(theta);
%%% y' = y*cos(theta)-x*sin(theta);

%% Describtion :

%% I : Input image
%% Sx & Sy : Variances along x and y-axes respectively
%% f : The frequency of the sinusoidal function
%% theta : The orientation of Gabor filter

%% G : The output filter as described above
%% gabout : The output filtered image

%%  Author : Ahmad poursaberi  e-mail : a.poursaberi@ece.ut.ac.ir
%%          Faulty of Engineering, Electrical&Computer Department,Tehran
%%          University,Iran,June 2004

function [G,gabout] = gaborfilter(I,Sx,Sy,f,theta);

if isa(I,'double')~=1
    I = double(I);
end

for x = -fix(Sx):fix(Sx)
    for y = -fix(Sy):fix(Sy)
        xPrime = x * cos(theta) + y * sin(theta);
        yPrime = y * cos(theta) - x * sin(theta);
        G(fix(Sx)+x+1,fix(Sy)+y+1) = exp(-.5*((xPrime/Sx)^2+(yPrime/Sy)^2))*cos(2*pi*f*xPrime);
    end
end

Imgabout = conv2(I,double(imag(G)),'same');
Regabout = conv2(I,double(real(G)),'same');

gabout = sqrt(Imgabout.*Imgabout + Regabout.*Regabout);


text:
function points = intersectLineCircle(line, circle)
%INTERSECTLINECIRCLE Intersection point(s) of a line and a circle.
%
%   INTERS = intersectLineCircle(LINE, CIRCLE);
%   Returns a 2-by-2-by-N array, containing on each row the coordinates of
%   an intersection point for each line-circle pair, i.e. INTERS(:,:,k)
%   contains the intersections between LINE(k,:) and CIRCLE(k,:).
%
%   If a line-circle pair does not intersect, the corresponding results are
%   set to NaN.
%
%   Example
%     % base point
%     center = [10 0];
%     % create vertical line
%     l1 = [center 0 1];
%     % circle
%     c1 = [center 5];
%     pts = intersectLineCircle(l1, c1)
%     pts =
%     10   -5
%     10    5
%     % draw the result
%     figure; clf; hold on;
%     axis([0 20 -10 10]);
%     drawLine(l1);
%     drawCircle(c1);
%     drawPoint(pts, 'rx');
%     axis equal;
%
%   See also
%   lines2d, circles2d, intersectLines, intersectCircles
%
%   References
%   http://local.wasp.uwa.edu.au/~pbourke/geometry/sphereline/
%   http://mathworld.wolfram.com/Circle-LineIntersection.html
%

% ------
% Authors: David Legland, JuanPi Carbajal
% E-mail: david.legland@inrae.fr, ajuanpi+dev@gmail.com
% Created: 2011-01-14, using Matlab 7.9.0.529 (R2009b)
% Copyright 2011-2022 INRA - Cepia Software Platform

% check size of inputs
nLines = size(line, 1);
nCircles = size(circle, 1);
if nLines ~= nCircles
  error ('matGeom:geom3d:invalidArguments', ...
      'Requires same number of lines and circles');
end

% center parameters
center = circle(:, 1:2);
radius = circle(:, 3);

% line parameters
dp = line(:, 1:2) - center;
vl = line(:, 3:4);

% coefficients of second order equation
a = sum(line(:, 3:4).^2, 2);
b = 2 * sum(dp .* vl, 2);
c = sum(dp.^2, 2) - radius.^2;

% discriminant
delta = b .^ 2 - 4 * a .* c;

points = nan(2, 2, nCircles);

valid = delta >= 0;

if any(valid)
    % compute roots (as a N-by-N-by-2 array)
    u = bsxfun(@plus, -b(valid), bsxfun(@times, [-1 1], sqrt(delta(valid))));
    u = bsxfun(@rdivide, u, a(valid)) / 2;

    if sum(valid) == 1
        points = [...
            line(1:2) + u(:,1) .* line(3:4); ...
            line(1:2) + u(:,2) .* line(3:4)];
    else
        tmp = [...
            line(valid, 1:2) + u(:,1) .* line(valid, 3:4) ...
            line(valid, 1:2) + u(:,2) .* line(valid, 3:4)].';
	    points(:, :, valid) = permute(reshape(tmp, [2, 2, nCircles]), [2 1 3]);
    end
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function eztools()
    %  put up eztools menu on current plot
    %
    %
    %  Richard G. Cobb  3/96

    %%%%%%%%%%%%%%%%%% Begin Figure Menu %%%%%%%%%%%%

    FigMenu = uimenu('Label','Figure');

    uimenu(FigMenu,'Label','Page Setup','Callback','pageset');
    uimenu(FigMenu,'Label','Labels','Separator','on','Callback','ezlabels');
    uimenu(FigMenu,'Label','Legend/Lines','Separator','on','Callback','ezlegend');

    uimenu(FigMenu,'Label','XY Locator','Separator','on','Callback','crosshair');
    uimenu(FigMenu,'Label','Add Arrows','Callback','ezarrow','Separator','on');

    % Colormap

    ColorEdit = uimenu(FigMenu,'Label','Edit ColorMap','Separator','on','Callback','makecmap(gcf,''fig'')');

    ColorMenu = uimenu(FigMenu,'Label','ColorMap','Separator','on');

    uimenu(ColorMenu,'Label','Default','Callback',...
        'colormap(''default'')');
    uimenu(ColorMenu,'Label','HSV','Callback',...
        'colormap(''hsv'')');
    uimenu(ColorMenu,'Label','summer','Callback',...
        'colormap(''summer'')');
    uimenu(ColorMenu,'Label','Gray','Callback',...
        'colormap(''gray'')');
    uimenu(ColorMenu,'Label','Hot','Callback',...
        'colormap(''hot'')');
    uimenu(ColorMenu,'Label','Cool','Callback',...
        'colormap(''cool'')');
    uimenu(ColorMenu,'Label','Copper','Callback',...
        'colormap(''copper'')');
    uimenu(ColorMenu,'Label','Pink','Callback',...
        'colormap(''pink'')');
    uimenu(ColorMenu,'Label','Jet','Callback',...
        'colormap(''jet'')');

    uimenu(FigMenu,'Label','Refresh','Separator','on','Callback','refresh');
    uimenu(FigMenu,'Label','WYSIWYG','Separator','on','Callback','wysiwyg');

    uimenu(FigMenu,'Label','Clear Figure',...
        'Callback','delete(findobj(''type'',''axes''))');

    %%%%%%%%%%%%%%%%%% Begin AXIS Menu %%%%%%%%%%%%

    AxMenu = uimenu('Label','Axis');
    uimenu(AxMenu,'Label','EZ-Axes','Callback','ezaxes');

    uimenu(AxMenu,'Label','Clear','Callback','cla');
    uimenu(AxMenu,'Label','Grid','Callback','grid');
    if(strcmp(axis('state'),'manual'))
        uimenu(AxMenu,'Label','Auto',...
            'Callback','axiscall(5,0)');
    else
        uimenu(AxMenu,'Label','Freeze',...
            'Callback','axiscall(5,0)');
    end
    AspectMenu = uimenu(AxMenu,'Label','Aspect');
    uimenu(AspectMenu,'Label','Normal','Callback','axis(''normal'')');
    uimenu(AspectMenu,'Label','Square','Callback','axis(''square'')');
    uimenu(AspectMenu,'Label','Image','Callback','axis(''image'')');
    uimenu(AspectMenu,'Label','Equal','Callback','axis(''equal'')');

    Xaxis = uimenu(AxMenu,'Label','X Opts');
    Yaxis = uimenu(AxMenu,'Label','Y Opts');
    Zaxis = uimenu(AxMenu,'Label','Z Opts');

    if(strcmp(get(gca,'XScale'),'linear'))
        uimenu(Xaxis,'Label','Log','Callback','axiscall(1,''X'')');
    else
        uimenu(Xaxis,'Label','Linear','Callback','axiscall(1,''X'')');
    end
    uimenu(Xaxis,'Label','Auto Min','Callback','axiscall(2,''X'')');
    uimenu(Xaxis,'Label','Auto Max','Callback','axiscall(3,''X'')');
    if(strcmp(get(gca,'YScale'),'linear'))
        uimenu(Yaxis,'Label','Log','Callback','axiscall(1,''Y'')');
    else
        uimenu(Yaxis,'Label','Linear','Callback','axiscall(1,''Y'')');
    end
    uimenu(Yaxis,'Label','Auto Min','Callback','axiscall(2,''Y'')');
    uimenu(Yaxis,'Label','Auto Max','Callback','axiscall(3,''Y'')');
    if(strcmp(get(gca,'ZScale'),'linear'))
        uimenu(Zaxis,'Label','Log','Callback','axiscall(1,''Z'')');
    else
        uimenu(Zaxis,'Label','Linear','Callback','axiscall(1,''Z'')');
    end
    uimenu(Zaxis,'Label','Auto Min','Callback','axiscall(2,''Z'')');
    uimenu(Zaxis,'Label','Auto Max','Callback','axiscall(3,''Z'')');

    if(ishold)
        uimenu(AxMenu,'Label','Hold off','Callback','axiscall(4,0)');
    else
        uimenu(AxMenu,'Label','Hold on','Callback','axiscall(4,0)');
    end

    uimenu(AxMenu,'Label','Axes Vis','Callback','axiscall(12,[])');

    %  Zoom Menu

    %uimenu(AxMenu,'Label','Zoom Box','Callback','subzoom');

    uimenu(AxMenu,'Label','Zoom ON','Callback','zoomrb on','Checked','off');
    uimenu(AxMenu,'Label','Zoom OFF','Callback','zoomrb off','Checked','off');

    uimenu(AxMenu,'Label','X-Y Sliders','Callback','axiscall(9,[])',...
        'Checked','off');

    uimenu(AxMenu,'Label','3-D Viewer','Callback','viewer','Separator','on');
    %New Zoom InteractiveMouse
    %uimenu(AxMenu,'Label','Mouse Zoom','Separator','on','Callback','interactivemouse');

    Zoom3d = uimenu(AxMenu,'Label','3-D Zoom');
    maglevel = 2;
    set(Zoom3d,'UserData',[1;1;1;maglevel]);
    uimenu(Zoom3d,'Label','Zoom in','Callback','zoom3d(1)');
    uimenu(Zoom3d,'Label','Zoom out','Callback','zoom3d(2)');
    uimenu(Zoom3d,'Label','X Axis','Checked','on','Callback',...
        'axiscall(6,''X'')','Separator','on');
    uimenu(Zoom3d,'Label','Y Axis','Checked','on','Callback',...
        'axiscall(6,''Y'')');
    uimenu(Zoom3d,'Label','Z Axis','Checked','on','Callback',...
        'axiscall(6,''Z'')');
    uimenu(Zoom3d,'Label',['Mag Level -' num2str(maglevel) '-'],'Callback',...
        'axiscall(7,0)','Separator','on');
    set(gcf,'WindowButtonDownFcn','tsel');

author: CelsoReyes


text:
function [IFF]=openiff(fid,LEN)
% OPENIFF is an auxillary function to SOPEN for
% opening of IFF files
%
% Use SOPEN instead of OPENIFF
%
% See also: fopen, SOPEN
%
% References:

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

%	$Id: openiff.m 2205 2009-10-27 12:18:15Z schloegl $
%	Copyright (C) 2004,2005,2007 by Alois Schloegl <a.schloegl@ieee.org>
%    	This is part of the BIOSIG-toolbox http://biosig.sf.net/

FLAG.CLOSE = 0;
if ischar(fid)
        fid = fopen(fid,'r');
        FLAG.CLOSE = 1;
end;

if nargin<2,
        LEN = 8;
end;

IFF = [];
K   = 0;
K1  = 0;

[tmp,c] = fread(fid,[1,4],'uint8');
while ((LEN>0) | isnan(LEN)) & (c>0),
        tag     = char(tmp);
        tagsize = fread(fid,1,'uint32');        % which size
        tagsize0= tagsize + rem(tagsize,2);
        filepos = ftell(fid);
        LEN = LEN - 8;

%       fprintf(1,'tag: %6s\tpos: %8i\tsize: %8i\n',tag,filepos,tagsize);

        if 0,
                VAL = openiff(fid,tagsize);
        elseif strcmp(tag,'FORM')
                [tmp,c] = fread(fid,[1,4],'uint8');tag,
                VAL = setfield([],char(tmp),openiff(fid,tagsize-4));
        elseif 0,strcmp(tag,'RIFF')
                VAL = openiff(fid,tagsize);
        elseif strcmp(tag,'RIFF')
                [tmp,c] = fread(fid,[1,4],'uint8');
                %val = fread(fid,tagsize-4,'uint8');
                val = openiff(fid,tagsize-4);
                VAL = setfield([],char(tmp),val);
        elseif strcmp(tag,'MThd')
                VAL.MThd = fread(fid,tagsize,'uchar');
                VAL.MIDI = openiff(fid,NaN);
                %LEN = NaN;
                %VAL.MIDI = tmp;

        elseif strcmp(tag,'LIST');	% CNT_RIFF (EEP 3.1)
                [tmp,c] = fread(fid,[1,4],'uint8');
                VAL = setfield([],char(tmp),openiff(fid,tagsize-4));

        elseif strcmp(tag,'chan')	% CNT_RIFF (EEP 3.1)
                VAL = fread(fid,[1,tagsize/2],'uint16');
        elseif strcmp(tag,'info')	% CNT_RIFF (EEP 3.1)
                VAL = fread(fid,[1,tagsize],'*char');
        elseif strcmp(tag,'eeph')	% CNT_RIFF (EEP 3.1)
                VAL = fread(fid,[1,tagsize],'*char');
        elseif strncmp(tag,'ep  ',4),	% CNT_RIFF (EEP 3.1)
                VAL = fread(fid,[1,tagsize/4],'uint32');
        elseif strcmp(tag,'hdrl')
                [tmp,c] = fread(fid,[1,4],'uint8');
                VAL = setfield([],char(tmp),openiff(fid,tagsize-4));

        elseif strcmp(tag,'CAT ')
                VAL = openiff(fid,tagsize);

        elseif strcmp(tag,'data')
        %        LEN = fread(fid,1,'uint32');
                VAL = fread(fid,[1,LEN],'*uint8');
                %[tmp,c] = fread(fid,[1,4],'*char');
                %LEN = fread(fid,1,'int32');
                %VAL.data = fread(fid,LEN,'uint8');
                %VAL = openiff(fid,tagsize);
        elseif strncmp(tag,'(c)',3)
                tag = 'Copyright';
                VAL = char(fread(fid,tagsize,'uchar')');
                %VAL.CopyRight = char(VAL);
        else
                if 0,tagsize<1024*2,
                        VAL = fread(fid,tagsize,'uchar');
                else
			VAL = [];
                        VAL.handle = ftell(fid);
                        VAL.size = tagsize;
                        status = fseek(fid,tagsize,'cof');
                        if status, return; end;
                end
        end;

        if strcmp(tag(3:4),'dc')
                K = K+1;
                ix = (tag(1:2)-48)*[16;1]+1;
                IFF.dc{K,ix} = VAL;
        elseif strcmp(tag(3:4),'wb')
                K1 = K1+1;
                ix = (tag(1:2)-48)*[16;1]+1;
                IFF.wb{K1,ix} = VAL;
        else
                %try,
	                IFF = setfield(IFF,deblank(tag),VAL);
                %catch,
                %end;
        end;
        status = fseek(fid,filepos+tagsize0,'bof');
        LEN = LEN - tagsize0;
        [tmp,c] = fread(fid,[1,4],'uint8');
end;

if ~isfield(IFF,'MThd'), % do not check MIDI files
        if ((LEN ~= 0) & c),
                fprintf(2,'Warning OPENIFF: LEN=%i %i %i %i %s  %i\n',LEN,filepos,tagsize0,ftell(fid),char(tmp),c);
        end;
end;

if FLAG.CLOSE,
        fclose(fid);
end;

author: PatternRecognition

repo: OpenBMI


text:
function [ind,t0,s0,t0close,s0close] = crossing(S,t,level,imeth)
% CROSSING find the crossings of a given level of a signal
%   ind = CROSSING(S) returns an index vector ind, the signal
%   S crosses zero at ind or at between ind and ind+1
%   [ind,t0] = CROSSING(S,t) additionally returns a time
%   vector t0 of the zero crossings of the signal S. The crossing
%   times are linearly interpolated between the given times t
%   [ind,t0] = CROSSING(S,t,level) returns the crossings of the
%   given level instead of the zero crossings
%   ind = CROSSING(S,[],level) as above but without time interpolation
%   [ind,t0] = CROSSING(S,t,level,par) allows additional parameters
%   par = {'none'|'linear'}.
%	With interpolation turned off (par = 'none') this function always
%	returns the value left of the zero (the data point thats nearest
%   to the zero AND smaller than the zero crossing).
%
%	[ind,t0,s0] = ... also returns the data vector corresponding to
%	the t0 values.
%
%	[ind,t0,s0,t0close,s0close] additionally returns the data points
%	closest to a zero crossing in the arrays t0close and s0close.
%
%	This version has been revised incorporating the good and valuable
%	bugfixes given by users on Matlabcentral. Special thanks to
%	Howard Fishman, Christian Rothleitner, Jonathan Kellogg, and
%	Zach Lewis for their input.

% Steffen Brueckner, 2002-09-25
% Steffen Brueckner, 2007-08-27		revised version

% Copyright (c) Steffen Brueckner, 2002-2007
% brueckner@sbrs.net

% check the number of input arguments
error(nargchk(1,4,nargin));

% check the time vector input for consistency
if nargin < 2 || isempty(t)
	% if no time vector is given, use the index vector as time
    t = 1:length(S);
elseif length(t) ~= length(S)
	% if S and t are not of the same length, throw an error
    error('t and S must be of identical length!');
end

% check the level input
if nargin < 3
	% set standard value 0, if level is not given
    level = 0;
end

% check interpolation method input
if nargin < 4
    imeth = 'linear';
end

% make row vectors
t = t(:)';
S = S(:)';

% always search for zeros. So if we want the crossing of
% any other threshold value "level", we subtract it from
% the values and search for zeros.
S   = S - level;

% first look for exact zeros
ind0 = find( S == 0 );

% then look for zero crossings between data points
S1 = S(1:end-1) .* S(2:end);
ind1 = find( S1 < 0 );

% bring exact zeros and "in-between" zeros together
ind = sort([ind0 ind1]);

% and pick the associated time values
t0 = t(ind);
s0 = S(ind);

if strcmp(imeth,'linear')
    % linear interpolation of crossing
    for ii=1:length(t0)
        if abs(S(ind(ii))) > eps(S(ind(ii)))
            % interpolate only when data point is not already zero
            NUM = (t(ind(ii)+1) - t(ind(ii)));
            DEN = (S(ind(ii)+1) - S(ind(ii)));
            DELTA =  NUM / DEN;
            t0(ii) = t0(ii) - S(ind(ii)) * DELTA;
            % I'm a bad person, so I simply set the value to zero
            % instead of calculating the perfect number ;)
            s0(ii) = 0;
        end
    end
end

% Addition:
% Some people like to get the data points closest to the zero crossing,
% so we return these as well
[CC,II] = min(abs([S(ind-1) ; S(ind) ; S(ind+1)]),[],1);
ind2 = ind + (II-2); %update indices

t0close = t(ind2);
s0close = S(ind2);

author: Sable

repo: mcbench-benchmarks


text:
function [ n1, n2 ] = year_to_dominical_common ( y )

%*****************************************************************************80
%
%% YEAR_TO_DOMINICAL_COMMON: dominical numbers, Common calendar.
%
%  Discussion:
%
%    The Julian calendar calculations are used through the year 1582,
%    and the Gregorian thereafter.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, the year.
%
%    Output, integer N1, N2, the dominical numbers for the year.
%    If Y is a leap year, then N1 applies before March 1, and N2 after.
%    If Y is not a leap year, then N1 applies throughout the year,
%    and N2 is returned as N1.
%
  if ( y <= 1582 )
    [ n1, n2 ] = year_to_dominical_julian ( y );
  else
    [ n1, n2 ] = year_to_dominical_gregorian ( y );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [K, P] = biasVardistPsi1Compute(biaskern, vardist, Z)

% BIASVARDISTPSI1COMPUTE one line description
% FORMAT
% DESC description
% RETURN K : description
% RETURN P : description
% ARG biasKern : the kernel structure associated with the white kernel.
% ARG vardist : description
% ARG Z : description
%
%
% SEEALSO : others
%
%
% COPYRIGHT : Michalis K. Titsias, 2009
%

% VARGPLVM

K = repmat(biaskern.variance,size(vardist.means,1),size(Z,1));

P = [];

author: SheffieldML

repo: vargplvm


text:
function [s1,s2] = paramSize(net)
    s1 = 0;s2=0;
    p = net.getParamIndex('dimred_filter') ;
    for i = 1:p-1
        s1 = s1 + prod(size(net.params(i).value));
    end
    s1 = s1 * 4/1024/1024;
    p2 = net.getParamIndex('predcls_filter');
    for i = p:p2
        s2 = s2 + prod(size(net.params(i).value));
    end
    s2 = s2 * 4/1024/1024;
end

author: HuiZeng

repo: Grid-Anchor-based-Image-Cropping


text:
function results = vl_test_grad(varargin)
% VL_TEST_GRAD
vl_test_init ;

function s = setup()
s.I = rand(150,253) ;
s.I_small = rand(2,2) ;

function test_equiv(s)
vl_assert_equal(gradient(s.I), vl_grad(s.I)) ;

function test_equiv_small(s)
vl_assert_equal(gradient(s.I_small), vl_grad(s.I_small)) ;

function test_equiv_forward(s)
Ix = diff(s.I,2,1) ;
Iy = diff(s.I,2,1) ;

vl_assert_equal(gradient(s.I_small), vl_grad(s.I_small)) ;

author: Cloud-CV

repo: object-proposals


text:
% Inference on a conditional Gaussian model

% Make the following polytree, where all arcs point down

% 1   2
%  \ /
%   3
%  / \
% 4   5

N = 5;
dag = zeros(N,N);
dag(1,3) = 1;
dag(2,3) = 1;
dag(3, [4 5]) = 1;

ns = [2 1 2 1 2];

dnodes = 1;
%onodes = [1 5];
bnet = mk_bnet(dag, ns, 'discrete', dnodes, 'observed', dnodes);

bnet.CPD{1} = tabular_CPD(bnet, 1);
for i=2:N
  bnet.CPD{i} = gaussian_CPD(bnet, i);
end

engine = {};
engine{end+1} = jtree_inf_engine(bnet);
engine{end+1} = pearl_inf_engine(bnet, 'protocol', 'parallel');

[time, engine] = cmp_inference_static(bnet, engine, 'maximize', 0, 'check_ll', 0, ...
				      'singletons_only', 0, 'observed', [1 3]);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function Data = ParamStructData(S)

%FUNCTION:
%   This function takes a data struct, where every field contains a (N x M)
%   matrix of numbers, and store the names of each field, the size of each
%   field, and the total number of elements. This is then later used to
%   pack everything into a struct.
%
%INPUTS:
%   S = the data struct of interest. Each field must contain a matrix of
%       numbers (no strings or other weird data types)
%
%OUTUTS:
%   Data = struct with the added fields
%       Data.Names = a (L x 1) cell array of field names
%       Data.Size = a (L x 2) matrix of the size of these fields
%       Data.Length = a scalar number of elements in S

    Names = fieldnames(S);
    Sizes = zeros(length(Names),2);
    Idx = zeros(length(Names),2);

    for i=1:length(Names)
       Sizes(i,:) = size(S.(Names{i}));
    end

    %Figure out the start and end indicies of each field.
    L = Sizes(:,1).*Sizes(:,2);
    Idx(1,1) = 1;
    for i=2:length(L)
       Idx(i-1,2) = Idx(i-1,1)-1+L(i-1);
       Idx(i,1) = Idx(i-1,2) + 1;
    end
    Idx(end,2) = sum(L);

    %Store the results:
    Data.Names = Names;
    Data.Sizes = Sizes;
    Data.Idx = Idx;

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function ringMaskM = getSurfaceRing(structNum,radius,planC)
%function mask3M = getSurfaceRing(structNum,radius,planC)
%
%This function returns surface ring of radius "radius" for structNum
%
%APA, 05/19/2020
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

if ~exist('planC','var')
    global planC
end

% Get Surface mask
% mask3M = getStructSurface(structNum,planC);

% Get surface mask of contracted unionStrNum
xyDownsampleIndex = 1; % no downsampling
contractedStrM = getSurfaceContract(structNum,radius/2,xyDownsampleIndex,planC);
expandedStrM = getSurfaceExpand(structNum,radius/2,xyDownsampleIndex,planC);

ringMaskM = expandedStrM - contractedStrM;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function lstOut = pyrLst2OP(lstIn1, lstIn2, fun)
%
%
%        lstOut = pyrLst2OP(lstIn1, lstIn2, fun)
%
%
%        Input:
%           -lstIn1: a list of image pyramids (first operand of fun)
%           -lstIn2: a list of image pyramids (second operand of fun)
%           -fun: a function to apply to lstIn1 and lstIn2
%
%        Output:
%           -lstOut: the result of the function
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%

n = length(lstIn1);

lstOut = [];

for i=1:n
    p = fun(lstIn1(i), lstIn2(i));
    lstOut = [lstOut, p];
end

end

author: banterle


text:
%% Analyzing Neural Time Series Data
% Matlab code for Chapter 33
% Mike X Cohen
%
% This code accompanies the book, titled "Analyzing Neural Time Series Data"
% (MIT Press). Using the code without following the book may lead to confusion,
% incorrect data analyses, and misinterpretations of results.
% Mike X Cohen assumes no responsibility for inappropriate or incorrect use of this code.

%% Figure 33.1

figure

subplot(121)
plot(normpdf(-4:.001:4))
axis tight

subplot(122)
a = randn(1000,1);
hist(a,50)
% this histogram was also used in figure 33.3

disp([ 'p_n = ' num2str(sum(a>2)/1000) ])
disp([ 'p_z = ' num2str(1-normcdf(2))  ])

%% Figure 33.5/6

% These figures are generated in the code for chapter 34

%% Figure 33.8

% introduction to bwlabeln and bwconncomp

% create 2D smoothing kernel
[xi,yi] = meshgrid(-10:10,-10:10);
zi = xi.^2+yi.^2;
zi = 1-(zi/max(zi(:)));

% create a random smoothed map
map = conv2(randn(100),zi,'same');

% threshold map at an arbitrary value
mapt = map;
mapt(abs(map)<range(map(:))/4) = 0;

% get labeled map via bwlabeln (in image processing toolbox)
[mapl,nblobs] = bwlabeln(mapt);
% mapl is a map of size <map> that contains numbers at each cluster.
% To extract information from clusters:
clustcount = zeros(1,nblobs);
clustsum   = zeros(1,nblobs);
for i=1:nblobs
    clustcount(i) = sum(mapl(:)==i);
    clustsum(i)   = sum(map(mapl(:)==i));
end

% bwconncomp works slightly differently, but will give similar information
blobinfo = bwconncomp(mapt);
% blobinfo is a structure that contains coordinates for each cluster rather than a map
% To extract information from clusters:
clustcount = zeros(1,nblobs);
clustsum   = zeros(1,nblobs);
for i=1:nblobs
    clustcount(i) = numel(blobinfo.PixelIdxList{i});
    clustsum(i)   = sum(map(blobinfo.PixelIdxList{i}));
end
% cluster count can be done faster using cellfun:
clustercount = cellfun(@numel,blobinfo.PixelIdxList);

figure
subplot(131)
imagesc(map), axis square
title('original')

subplot(132)
imagesc(mapt), axis square
title('thresholded')

subplot(133)
imagesc(mapl), axis square
title([ 'labeled (' num2str(nblobs) ' clusters in total)' ])

%% Figure 33.10

% This script calls the function fdr.m, which is included in the online
% code and was downloaded (and slightly modified) in summer 2012 from http://www-personal.umich.edu/~nichols/FDR/FDR.m

nsigs = round(linspace(1,500,80));
nnons = round(linspace(1,500,100));

fdrpvals = zeros(20,length(nsigs),length(nnons));

for iteri=1:20
    for i=1:length(nsigs)
        for j=1:length(nnons)

            pvals = [ rand(1,nsigs(i))*.05 rand(1,nnons(j))*.5+.05 ];
            temp = fdr(pvals,.05);
            if isempty(temp)
                fdrpvals(iteri,i,j) = NaN;
            else
                fdrpvals(iteri,i,j) = temp;
            end
        end
    end
end

fdrpvals = squeeze(nanmean(fdrpvals));

figure
imagesc(fdrpvals)
set(gca,'clim',[0 .05])
xlabel('number of non-significant p-values')
ylabel('number of significant p-values')

figure
subplot(211)
plot(nanmean(fdrpvals,1))
xlabel('number of non-significant p-values')
ylabel('critical p-value')
set(gca,'ylim',[0 .05])

subplot(212)
plot(nanmean(fdrpvals,2))
xlabel('number of significant p-values')
ylabel('critical p-value')
set(gca,'ylim',[0 .05])

%% end.

author: mikexcohen

repo: AnalyzingNeuralTimeSeries


text:
function x = FMdemod(y,beta)

[r, c] = size(y);
if r*c == 0
    y = [];
    return;
end;
if (r == 1)
    y = y(:);
    len = c;
else
    len = r;
end;

Fc = 40000;
Fs = 8000*20;
pi2 = 2*pi;

    [num, den] = butter(5, Fc * 2 / Fs);

    sen = 2*pi*beta*100;

    %pre-process the filter.
    if abs(den(1)) < eps
        error('First denominator filter coefficient must be non-zero.');
    else
        num = num/den(1);
        if (length(den) > 1)
            den = - den(2:length(den)) / den(1);
        else
            den = 0;
        end;
        num = num(:)';
        den = den(:)';
    end;
    len_den = length(den);
    len_num = length(num);

    x = y;
    y = 2 * y;

    ini_phase = pi/2;
    for ii = 1 : size(y, 2)
        z1 = zeros(length(den), 1);
        s1 = zeros(len_num, 1);
        intgl = 0;

        memo = 0;
        for i = 1:size(y, 1)
            %start with the zero-initial condition integer.
            vco_out = cos(pi2 * intgl+ini_phase);
            if len_num > 1
                s1 = [y(i, ii) * vco_out; s1(1:len_num-1)];
            else
                s1 = y(i, ii);
            end
            tmp = num * s1 + den * z1;
            if len_den > 1
                z1 = [tmp; z1(1:len_den-1)];
            else
                z1 = tmp;
            end;
            intgl = rem(((tmp*sen + Fc)/ Fs + intgl), 1);
            x(i, ii) = tmp;
        end;
    end;
    x = x;
x = decimate(x,20);


text:
function [B,C]=nearestKronProd(A,m1,m2,n1,n2)
%%NEARESTKRONPRODSVD Given an (m1*m2)X(n1*n2) matrix, determine the set of
%           matrices B and C that minimize norm(A-kron(B,C),'fro'). B is
%           m1Xn1 and C is m2Xn2. B and C are not unique. One can always
%           replace them with x*B and (1/x)*C for some arbitrary scalar x.
%
%INPUTS: A An (m1*m2)X(n1*n2) real matrix which is an m1Xn1 collection of
%          blocks of size m2Xn2. It has the form
%          A=[A_{1,1} ... A_{1,n1};
%             ...     ...      ...;
%             A_{m1,1} ... A_{m1,n1}];
% m1,m2,n1,n2 The dimensions specifying A and its submatrices as described
%          above. These specify the dimensions of B and C.
%
%OUTPUTS: B The m1Xn1 left matrix in the Kronecker product.
%         C The m2Xn2 right matrix in the Kronecker product.
%
%This function implements the algorithm described in Section 12.3.6 of [1].
%
%EXAMPLE:
%Here, we compute A explicitely from a Kronecker product of B and C
%matrices. The nearest Kronecker product obtained by this function should
%thus have zero error, which, within finite precision limits, is
%demonstrated.
% m1=2;
% m2=3;
% n1=5;
% n2=4;
% B=randn(m1,n1);
% C=randn(m2,n2);
% A=kron(B,C);
%
% [BK,CK]=nearestKronProd(A,m1,m2,n1,n2);
% norm(A-kron(BK,CK),'fro')%The error will be close to 0.
%
%REFERENCES:
%[1] G. H. Golub and C. F. Van Loan, Matrix Computations, 4th ed.
%    Baltimore: Johns Hopkins University Press, 2013.
%
%September 2019 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

[sigma,U,V]=KroneckerProdSVD(A,m1,n1,m2,n2);

sqrtSigma1=sqrt(sigma(1));

B=sqrtSigma1*U(:,:,1);
C=sqrtSigma1*V(:,:,1);

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function plotResults(price,x,N)
% Copyright (c) 2012, MathWorks, Inc.
%% Create a new figure
figure;

%% Top plot of price data
subplot(2,1,1);
plot(price)
% Add labels
xlabel('Time (hrs)');
ylabel('Price ($/MWh)');
title('Electricity Prices');

%% Bottom plot is flow rates
subplot(2,1,2);
plot(x(1:N),'-*b');
hold on;
plot(x(N+1:end),'-*r');
% Add labels
xlabel('Flow (CFS)');
ylabel('Time (hrs)');
title('Spill and Turbine Flow Rates');
% Add legend in upper left (NorthWest) corner
legend('Turbine Flow','Spill Flow','Location','NorthWest');


text:
function check = uniform_discrete_check ( a, b )

%*****************************************************************************80
%
%% UNIFORM_DISCRETE_CHECK checks the parameters of the Uniform discrete CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer A, B, the parameters of the PDF.
%    A <= B.
%
%    Output, logical CHECK, is true if the parameters are legal.
%
  if ( b < a )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'UNIFORM_DISCRETE_CHECK - Fatal error!\n' );
    fprintf ( 1, '  B < A.\n' );
    check = 0;
  end

  check = 1;

  return
end


text:
function errores = errors

source = evalin('base','source');
decout = evalin('base','decout');
N = evalin('base','N');

errores = sum(abs(source - decout));


text:
function [m] = m2m(m)
% Convert length from British metres to American meters.
m = m;


text:
function [ n_data, x, fx ] = arcsin_values ( n_data )

%*****************************************************************************80
%
%% ARCSIN_VALUES returns some values of the arc sine function.
%
%  Discussion:
%
%    In Mathematica, the function can be evaluated by:
%
%      ArcSin[x]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 June 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz, Irene Stegun,
%    Handbook of Mathematical Functions,
%    National Bureau of Standards, 1964,
%    ISBN: 0-486-61272-4,
%    LC: QA47.A34.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Cambridge University Press, 1999,
%    ISBN: 0-521-64314-7,
%    LC: QA76.95.W65.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real X, the argument of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 12;

  fx_vec = [ ...
    -0.10016742116155979635, ...
     0.00000000000000000000, ...
     0.10016742116155979635, ...
     0.20135792079033079146, ...
     0.30469265401539750797, ...
     0.41151684606748801938, ...
     0.52359877559829887308, ...
     0.64350110879328438680, ...
     0.77539749661075306374, ...
     0.92729521800161223243, ...
     1.1197695149986341867, ...
     1.5707963267948966192 ];

  x_vec = [ ...
    -0.1, ...
     0.0, ...
     0.1, ...
     0.2, ...
     0.3, ...
     0.4, ...
     0.5, ...
     0.6, ...
     0.7, ...
     0.8, ...
     0.9, ...
     1.0 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    x = 0.0;
    fx = 0.0;
  else
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end


text:
% ttest2_cell() - compute unpaired t-test. Allow fast computation of
%                 multiple t-test using matrix manipulation.
%
% Usage:
%    >> [F df] = ttest2_cell( { a b } );
%    >> [F df] = ttest2_cell(a, b);
%    >> [F df] = ttest2_cell(a, b, 'inhomogenous');
%
% Inputs:
%   a,b       = data consisting of UNPAIRED arrays to be compared. The last
%               dimension of the data array is used to compute the t-test.
%   'inhomogenous' = use computation for the degree of freedom using
%                    inhomogenous variance. By default the computation of
%                    the degree of freedom is done with homogenous
%                    variances.
%
% Outputs:
%   T   - T-value
%   df  - degree of freedom (array)
%
% Example:
%   a = { rand(1,10) rand(1,10)+0.5 }
%   [T df] = ttest2_cell(a)
%   signif = 2*tcdf(-abs(T), df(1))
%
%   % for comparison, the same using the Matlab t-test function
%   [h p ci stats] = ttest2(a{1}', a{2}');
%   [ stats.tstat' p]
%
%   % fast computation (fMRI scanner volume 100x100x100 and 10 control
%   % subjects and 12 test subjects). The computation itself takes 0.5
%   % seconds instead of  half an hour using the standard approach (1000000
%   % loops and Matlab  t-test function)
%   a = rand(100,100,100,10); b = rand(100,100,100,10);
%   [F df] = ttest_cell({ a b });
%
% Author: Arnaud Delorme, SCCN/INC/UCSD, La Jolla, 2005
%         (thank you to G. Rousselet for providing the formula for
%         inhomogenous variances).
%
% Reference:
%   Schaum's outlines in statistics (3rd edition). 1999. Mc Graw-Hill.
%   Howel, Statistical Methods for Psychology. 2009. Wadsworth Publishing.

% Copyright (C) Arnaud Delorme
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function [tval, df] = ttest2_cell(a,b,c) % assumes equal variances

    if nargin < 1
        help ttest2_cell;
        return;
    end;

    homogenous = 'homogenous';
    if nargin > 1 && isstr(b)
        homogenous = b;
    end;
    if nargin > 2 && isstr(c)
        homogenous = c;
    end;
    if iscell(a),
        b = a{2};
        a = a{1};
    end;
    if ~strcmpi(homogenous, 'inhomogenous') && ~strcmpi(homogenous, 'homogenous')
        error('Value for homogenous parameter can only be ''homogenous'' or ''inhomogenous''');
    end;

    nd    = myndims(a);
    na    = size(a, nd);
    nb    = size(b, nd);
    meana = mymean(a, nd);
    meanb = mymean(b, nd);

    if strcmpi(homogenous, 'inhomogenous')
        % inhomogenous variance from Howel, 2009, "Statistical Methods for Psychology"
        % thank you to G. Rousselet for providing these formulas
        m  = meana - meanb;
        s1 = var(a,0,nd) ./ na;
        s2 = var(b,0,nd) ./ nb;
        se = sqrt(s1 + s2);
        sd = sqrt([s1.*na, s2.*nb]);
        tval = m ./ se;

        df = ((s1 + s2).^2) ./ ((s1.^2 ./ (na-1) + s2.^2 ./ (nb-1)));
    else
        sda   = mystd(a, [], nd);
        sdb   = mystd(b, [], nd);
        sp    = sqrt(((na-1)*sda.^2+(nb-1)*sdb.^2)/(na+nb-2));
        tval  = (meana-meanb)./sp/sqrt(1/na+1/nb);
        df    = na+nb-2;
    end;

    % check values againg Matlab statistics toolbox
    % [h p ci stats] = ttest2(a', b');
    % [ tval stats.tstat' ]

function val = myndims(a)
    if ndims(a) > 2
        val = ndims(a);
    else
        if size(a,1) == 1,
            val = 2;
        elseif size(a,2) == 1,
            val = 1;
        else
            val = 2;
        end;
    end;

function res = mymean( data, varargin) % deal with complex numbers
    res = mean( data, varargin{:});
    if ~isreal(data)
        res = abs( res );
    end;

function res = mystd( data, varargin) % deal with complex numbers
    if ~isreal(data)
        res = std( abs(data), varargin{:});
    else
        res = sqrt(sum( bsxfun(@minus, data, mean( data, varargin{2})).^2, varargin{2})/(size(data,varargin{2})-1)); % 8 percent speedup
        %res = std( data, varargin{:});
    end;


text:
function [ a, b] = p05_interval ( )

%*****************************************************************************80
%
%% P05_INTERVAL returns a starting interval for optimization for problem 5.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 February 2002
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real A, B, two points defining an interval in which
%    the local minimizer should be sought.
%
  a = 0.0002;
  b = 2.0;

  return
end


text:
function estimatedLabels = GMLVQ_classify(Data, model)

%GMLVQ_classify.m - classifies the given data with the given model
%
% estimatedLabels = GMLVQ_classify(Data, model)
%
% Input and output arguments:
%  trainSet          (struct)
%                    (matrix) matrix with training samples in its rows
%  model    : GMLVQ model with prototypes w their labels c_w and the matrix omega
%
%  estimatedLabels       : the estimated labels
%
%  example for usage:
%  trainSet = [1,2,3;4,5,6;7,8,9];
%  trainLab = [1;1;2];
%  GMLVQ_model=gmlvq_core(trainSet,trainLab); % minimal parameters required
%  estimatedTrainLabels = GMLVQ_classify(trainSet, GMLVQ_model);
%  trainError = mean( trainLab ~= estimatedTrainLabels );
%
% Kerstin Bunte
% kerstin.bunte@googlemail.com
% Mon Nov 05 09:05:52 CEST 2012
%
% Conditions of GNU General Public License, version 2 and BSD License apply.
% See file 'license-gpl2.txt' and 'BSD_license.txt' enclosed in this package.
% Programs are not for use in critical applications!
%

% Contributed to SOM Toolbox vs2, December 3rd, 2012 by Alexander Schulz
% Copyright (c) Kerstin Bunte
% http://www.cis.hut.fi/projects/somtoolbox/

if isstruct(Data)
    Data = Data.data;
end

if ~isfield(model, 'w')
    model.c_w = som_label2num(model.labels);
    dist = computeDistance(Data, model.codebook, model);
else
    dist = computeDistance(Data, model.w, model);
end

[~, index] = min(dist,[],2);

estimatedLabels = model.c_w(index);

function distance = computeDistance(X, W, model)
nb_samples = size(X,1);
distance = zeros(nb_samples,length(model.c_w));

% tic;
% for i = 1:size(W,1)
%     delta = X - ones(P,1) * W(i,:);
%     delta = bsxfun(@minus, X, W(i,:));
%     delta(isnan(delta)) = 0;
    % Hadamard product: to skip unnecessary calculation between two different examples
%     distance(1:nb_samples,i) = sum( ((delta*model.omega'*model.omega).*delta) ,2 );
%     wi = W(i,:);
%     distance(1:nb_samples,i) = sum( (( (X-wi(ones(nb_samples,1),:)) *model.omega'*model.omega).*(X-wi(ones(nb_samples,1),:))) ,2 );
% end
% disp(toc);
if isfield(model,'psis')
    if length(model.psis)~=size(W,1)
        classes = unique(model.c_w);
        for i = 1:size(W,1)
            matrixIdx = classes==model.c_w(i);
%             delta = X-ones(nb_samples,1)*W(i,:);
%             distance(1:nb_samples,i) = sum( ((delta* model.psis{matrixIdx}'*model.psis{matrixIdx}).*delta) ,2 );
            distance(1:nb_samples,i) = sum((bsxfun(@minus, X, W(i,:))*model.psis{matrixIdx}').^2, 2);
        end
    else
        for i = 1:size(W,1)
%             delta = X-ones(nb_samples,1)*W(i,:);
%             distance(1:nb_samples,i) = sum( ((delta* model.psis{i}'*model.psis{i}).*delta) ,2 );
            distance(1:nb_samples,i) = sum((bsxfun(@minus, X, W(i,:))*model.psis{i}').^2, 2);
        end
    end
else
    if isfield(model,'lambda')
        for i = 1:size(W,1)
%             delta = X-ones(nb_samples,1)*W(i,:);
%             distance(1:nb_samples,i) = sum( (( delta *model.lambda).*delta) ,2 );
            delta = bsxfun(@minus, X, W(i,:));
            distance(1:nb_samples,i) = sum(bsxfun(@times,delta.^2,model.lambda), 2);
%             distance(1:nb_samples,i) = sum( (( delta *model.lambda).*delta) ,2 );
        end
    end
    if isfield(model,'omega')
        % tic;
        for i = 1:size(W,1)
%             delta = X-ones(nb_samples,1)*W(i,:);
%             distance(1:nb_samples,i) = sum( (( delta *model.omega'*model.omega).*delta) ,2 );
            distance(1:nb_samples,i) = sum((bsxfun(@minus, X, W(i,:))*model.omega').^2, 2);
        end
        % disp(toc);
    end
end

author: ilarinieminen

repo: SOM-Toolbox


text:
% NANMIN provides a replacement for MATLAB's nanmin that is almost
% compatible.

function [varargout] = nanmin(varargin)

% the standard MATLAB min function ignores NaNs from at least version 7.0.4 onwards
[varargout{1:nargout}] = min(varargin{:});

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
filename = 'BulkCaseSym';%'CantileverSquareNewFine';%'LshapeTriFine';%'CantileverSquareNew';'ArchTriFine';%'CantileverSquare';%'ArchTriFine';%'CantileverSquare';%'Lshape';%'LshapeTriFine';%'CantileverSquareSmall';%'';%'Lshape';'CantileverSquare';'Lshape';%'CantileverSquare';%'LshapeFine';%'Bridge_quad_coarse';%'Arch_quad_coarse';%'BridgeCool_Quadrilateral_Bilinear_Structured_Coarse';%'Bridge';%'CantileverSquareSmall';
ptype = 'MACRO';
initial_case = 'given';
m1 = 0.0101;
m2 = 0.0101;
%cost = {'compliance'};
cost = {'stressNorm'};
%cost = {'stressNorm','compliance'};
%weights = [0.55,0.45];
weights = 1;
constraint = {'volumeConstraint'};
filterType = 'P1';
constraint_case = 'EQUALITY';

Vfrac_initial = 0.3;
optimality_initial = 1e-5;
constr_initial = 1e-5;

Vfrac_final = 0.3;
optimality_final = 1e-5;
constr_final = 1e-5;

stressNormExponent_initial = 2;
stressNormExponent_final = 32;

optimizer = 'DualNestedInPrimal';
optimizerUnconstrained = 'PROJECTED GRADIENT';

designVariable = 'MicroParams';
ub = 0.989;
lb = 0.011;
homegenizedVariablesComputer = 'ByVademecum';
%
%vademecumFileName = 'SuperEllipseQMax';
%vademecumFileName = 'SuperEllipseQ2';
vademecumFileName = 'SuperEllipseQOptAnalytic';
%
% designVariable = 'Density';
% homegenizedVariablesComputer = 'ByInterpolation';
% method = 'SIMPALL';
% materialType = 'ISOTROPIC';

kfrac = 2;
nsteps = 100;

plotting = false;
printing = true;
monitoring = true;
monitoring_interval = 1;
maxiter = 1000;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [ parent_pop ] = evaluate_pop( parent_pop, obj_func )
%   This function evaluates a whole population with the supplied
%   objective function handle obj_func.

global nreal ;
global nobj ;
global ncon ;

obj_col = nreal + 1: nreal + nobj ;
parent_pop(:,obj_col) = 0 ;
parent_pop = obj_func(parent_pop);

if(ncon > 0)
    cv_col = nreal + nobj + ncon + 1 ;
    parent_pop(:,cv_col) = 0;
    cons = parent_pop(:,nreal+nobj+1:nreal+nobj+ncon);
    % parent_pop(:,cv_col) = sum(((cons < 0.0) .* cons),2);
    parent_pop(:,cv_col) = sum(bsxfun(@times, ...
                                    bsxfun(@lt, cons, 0.0), cons),2);
end
end

author: chudur-budur


text:
function out = feval(S, varargin)
%FEVAL   Evaluate the operator of a SPINOP (1, 2 or 3) at a CHEBFUN (1, 2, 3) or
%at a CHEBMATRIX.
%   OUT = FEVAL(S, U) for a CHEBFUN (1, 2, or 3) or CHEBMATRIX U applies S to U.
%
%   OUT = FEVAL(S, U, V, ...) for CHEBFUN'S (1, 2, 3) U, V, ... applies S to
%   the CHEBFUN'S.
%
% See also SPINOPERATOR/SUBSREF.

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

L = S.lin;
N = S.nonlin;

% CASE 1: S(U).
if ( nargin == 2 )

    u = varargin{1};

    if ( isa(u, 'chebfun') || isa(u, 'chebfun2') || isa(u, 'chebfun3') )

        out = L(u) + N(u);

    elseif ( isa(u, 'chebmatrix' ) )

        u = u.blocks.';
        out = L(u{:}) + N(u{:});

    end

% CASE 2: S(U,V,...).
else

    out = L(varargin{:}) + N(varargin{:});

end

% Convert to a CHEBMATRIX if it's a CHEBFUN2V or a CHEBFUN3V:
if ( isa(out, 'chebfun2v') || isa(out, 'chebfun3v') )

    u = chebmatrix(out(1));
    for k = 2:size(out, 1)
        u(k,1) = out(k);
    end
    out = u;

end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [X, fX, i] = kafbox_minimize(X, f, length, varargin)

% Minimize a differentiable multivariate function.
%
% Usage: [X, fX, i] = minimize(X, f, length, P1, P2, P3, ... )
%
% where the starting point is given by "X" (D by 1), and the function named in
% the string "f", must return a function value and a vector of partial
% derivatives of f wrt X, the "length" gives the length of the run: if it is
% positive, it gives the maximum number of line searches, if negative its
% absolute gives the maximum allowed number of function evaluations. You can
% (optionally) give "length" a second component, which will indicate the
% reduction in function value to be expected in the first line-search (defaults
% to 1.0). The parameters P1, P2, P3, ... are passed on to the function f.
%
% The function returns when either its length is up, or if no further progress
% can be made (ie, we are at a (local) minimum, or so close that due to
% numerical problems, we cannot get any closer). NOTE: If the function
% terminates within a few iterations, it could be an indication that the
% function values and derivatives are not consistent (ie, there may be a bug in
% the implementation of your "f" function). The function returns the found
% solution "X", a vector of function values "fX" indicating the progress made
% and "i" the number of iterations (line searches or function evaluations,
% depending on the sign of "length") used.
%
% The Polack-Ribiere flavour of conjugate gradients is used to compute search
% directions, and a line search using quadratic and cubic polynomial
% approximations and the Wolfe-Powell stopping criteria is used together with
% the slope ratio method for guessing initial step sizes. Additionally a bunch
% of checks are made to make sure that exploration is taking place and that
% extrapolation will not be unboundedly large.
%
% See also: checkgrad
%
% Copyright (C) 2001 - 2006 by Carl Edward Rasmussen (2006-09-08).

INT = 0.1;    % don't reevaluate within 0.1 of the limit of the current bracket
EXT = 3.0;                  % extrapolate maximum 3 times the current step-size
MAX = 20;                         % max 20 function evaluations per line search
RATIO = 10;                                       % maximum allowed slope ratio
SIG = 0.1; RHO = SIG/2; % SIG and RHO are the constants controlling the Wolfe-
% Powell conditions. SIG is the maximum allowed absolute ratio between
% previous and new slopes (derivatives in the search direction), thus setting
% SIG to low (positive) values forces higher precision in the line-searches.
% RHO is the minimum allowed fraction of the expected (from the slope at the
% initial point in the linesearch). Constants must satisfy 0 < RHO < SIG < 1.
% Tuning of SIG (depending on the nature of the function to be optimized) may
% speed up the minimization; it is probably not worth playing much with RHO.

% The code falls naturally into 3 parts, after the initial line search is
% started in the direction of steepest descent. 1) we first enter a while loop
% which uses point 1 (p1) and (p2) to compute an extrapolation (p3), until we
% have extrapolated far enough (Wolfe-Powell conditions). 2) if necessary, we
% enter the second loop which takes p2, p3 and p4 chooses the subinterval
% containing a (local) minimum, and interpolates it, unil an acceptable point
% is found (Wolfe-Powell conditions). Note, that points are always maintained
% in order p0 <= p1 <= p2 < p3 < p4. 3) compute a new search direction using
% conjugate gradients (Polack-Ribiere flavour), or revert to steepest if there
% was a problem in the previous line-search. Return the best value so far, if
% two consecutive line-searches fail, or whenever we run out of function
% evaluations or line-searches. During extrapolation, the "f" function may fail
% either with an error or returning Nan or Inf, and minimize should handle this
% gracefully.

if max(size(length)) == 2, red=length(2); length=length(1); else red=1; end
if length>0, S='Linesearch'; else S='Function evaluation'; end

i = 0;                                            % zero the run length counter
ls_failed = 0;                             % no previous line search has failed
[f0 df0] = feval(f, X, varargin{:});          % get function value and gradient
fX = f0;
i = i + (length<0);                                            % count epochs?!
s = -df0; d0 = -s'*s;           % initial search direction (steepest) and slope
x3 = red/(1-d0);                                  % initial step is red/(|s|+1)

while i < abs(length)                                      % while not finished
  i = i + (length>0);                                      % count iterations?!

  X0 = X; F0 = f0; dF0 = df0;                   % make a copy of current values
  if length>0, M = MAX; else M = min(MAX, -length-i); end

  while 1                             % keep extrapolating as long as necessary
    x2 = 0; f2 = f0; d2 = d0; f3 = f0; df3 = df0;
    success = 0;
    while ~success && M > 0
      try
        M = M - 1; i = i + (length<0);                         % count epochs?!
        [f3 df3] = feval(f, X+x3*s, varargin{:});
        if isnan(f3) || isinf(f3) || any(isnan(df3)+isinf(df3)), error(''), end
        success = 1;
      catch                                % catch any error which occured in f
        x3 = (x2+x3)/2;                                  % bisect and try again
      end
    end
    if f3 < F0, X0 = X+x3*s; F0 = f3; dF0 = df3; end         % keep best values
    d3 = df3'*s;                                                    % new slope
    if d3 > SIG*d0 || f3 > f0+x3*RHO*d0 || M == 0  % are we done extrapolating?
      break
    end
    x1 = x2; f1 = f2; d1 = d2;                        % move point 2 to point 1
    x2 = x3; f2 = f3; d2 = d3;                        % move point 3 to point 2
    A = 6*(f1-f2)+3*(d2+d1)*(x2-x1);                 % make cubic extrapolation
    B = 3*(f2-f1)-(2*d1+d2)*(x2-x1);
    x3 = x1-d1*(x2-x1)^2/(B+sqrt(B*B-A*d1*(x2-x1))); % num. error possible, ok!
    if ~isreal(x3) || isnan(x3) || isinf(x3) || x3 < 0 % num prob | wrong sign?
      x3 = x2*EXT;                                 % extrapolate maximum amount
    elseif x3 > x2*EXT                  % new point beyond extrapolation limit?
      x3 = x2*EXT;                                 % extrapolate maximum amount
    elseif x3 < x2+INT*(x2-x1)         % new point too close to previous point?
      x3 = x2+INT*(x2-x1);
    end
  end                                                       % end extrapolation

  while (abs(d3) > -SIG*d0 || f3 > f0+x3*RHO*d0) && M > 0  % keep interpolating
    if d3 > 0 || f3 > f0+x3*RHO*d0                         % choose subinterval
      x4 = x3; f4 = f3; d4 = d3;                      % move point 3 to point 4
    else
      x2 = x3; f2 = f3; d2 = d3;                      % move point 3 to point 2
    end
    if f4 > f0
      x3 = x2-(0.5*d2*(x4-x2)^2)/(f4-f2-d2*(x4-x2));  % quadratic interpolation
    else
      A = 6*(f2-f4)/(x4-x2)+3*(d4+d2);                    % cubic interpolation
      B = 3*(f4-f2)-(2*d2+d4)*(x4-x2);
      x3 = x2+(sqrt(B*B-A*d2*(x4-x2)^2)-B)/A;        % num. error possible, ok!
    end
    if isnan(x3) || isinf(x3)
      x3 = (x2+x4)/2;               % if we had a numerical problem then bisect
    end
    x3 = max(min(x3, x4-INT*(x4-x2)),x2+INT*(x4-x2));  % don't accept too close
    [f3 df3] = feval(f, X+x3*s, varargin{:});
    if f3 < F0, X0 = X+x3*s; F0 = f3; dF0 = df3; end         % keep best values
    M = M - 1; i = i + (length<0);                             % count epochs?!
    d3 = df3'*s;                                                    % new slope
  end                                                       % end interpolation

  if abs(d3) < -SIG*d0 && f3 < f0+x3*RHO*d0          % if line search succeeded
    X = X+x3*s; f0 = f3; fX = [fX' f0]';                     % update variables
    fprintf('%s %6i;  Value %4.6e\r', S, i, f0);
    s = (df3'*df3-df0'*df3)/(df0'*df0)*s - df3;   % Polack-Ribiere CG direction
    df0 = df3;                                               % swap derivatives
    d3 = d0; d0 = df0'*s;
    if d0 > 0                                      % new slope must be negative
      s = -df0; d0 = -s'*s;                  % otherwise use steepest direction
    end
    x3 = x3 * min(RATIO, d3/(d0-realmin));          % slope ratio but max RATIO
    ls_failed = 0;                              % this line search did not fail
  else
    X = X0; f0 = F0; df0 = dF0;                     % restore best point so far
    if ls_failed || i > abs(length)         % line search failed twice in a row
      break;                             % or we ran out of time, so we give up
    end
    s = -df0; d0 = -s'*s;                                        % try steepest
    x3 = 1/(1-d0);
    ls_failed = 1;                                    % this line search failed
  end
end
fprintf('\n');

repo: kafbox


text:
function varargout = MagDesignPanel(varargin)
% MAGNETDESIGNPANEL MATLAB code for MagDesignPanel.fig
%      MAGNETDESIGNPANEL, by itself, creates a new MAGNETDESIGNPANEL or raises the existing
%      singleton*.
%
%      H = MAGNETDESIGNPANEL returns the handle to a new MAGNETDESIGNPANEL or the handle to
%      the existing singleton*.
%
%      MAGNETDESIGNPANEL('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in MAGNETDESIGNPANEL.M with the given input arguments.
%
%      MAGNETDESIGNPANEL('Property','Value',...) creates a new MAGNETDESIGNPANEL or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before MagDesignPanel_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to MagDesignPanel_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help MagDesignPanel

% Last Modified by GUIDE v2.5 17-Apr-2014 16:33:01

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

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

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

handles.Simuh=varargin{1};
% Load Pulse element list
try
    handles.MagElemListStruct=DoParseXML([handles.Simuh.MRiLabPath filesep 'Macro' filesep 'MagElem' filesep 'MagElem.xml']);
    guidata(hObject,handles);
    Root=DoConvStruct2Tree(handles.MagElemListStruct);
    handles.MagElemListTree=uitree('v0', 'Root',Root,'SelectionChangeFcn', {@ChkMagElemListTreeNode,handles});
    set(handles.MagElemListTree.getUIContainer,'Units','normalized');
    set(handles.MagElemListTree.getUIContainer,'Position',[0.0,0.4,0.2,0.58]);
    handles.Attrh1=[];
catch ME
    errordlg('Mag element list can not be loaded!');
    close(hObject);
    return;
end

%Load tabs for Mag Simulation
try
    handles.MagSimuAttrStruct=DoParseXML([handles.Simuh.MRiLabPath filesep 'Macro' filesep 'MagElem' filesep 'MagSimuAttr.xml']);
catch ME
    errordlg('MagSimuAttr.xml file is missing or can not be loaded!');
    close(hObject);
    return;
end

handles.Setting_tabgroup=uitabgroup(handles.Setting_uipanel);
guidata(handles.MagDesignPanel_figure,handles);
for i=1:length(handles.MagSimuAttrStruct.Children)
    eval(['handles.' handles.MagSimuAttrStruct.Children(i).Name '_tab=uitab( handles.Setting_tabgroup,' '''title'',' '''' handles.MagSimuAttrStruct.Children(i).Name ''',''Units'',''normalized'');']);
    eval(['DoEditValue(handles,handles.' handles.MagSimuAttrStruct.Children(i).Name '_tab,handles.MagSimuAttrStruct.Children(' num2str(i) ').Attributes,2,[0.2,0.15,0.0,0.15,0.15,0.05]);']);
    handles=guidata(handles.MagDesignPanel_figure);
end

% Create grid and display structure
VMmg.xdimres=str2num(get(handles.Attrh2.('XDimRes'),'String'));
VMmg.ydimres=str2num(get(handles.Attrh2.('YDimRes'),'String'));
VMmg.zdimres=str2num(get(handles.Attrh2.('ZDimRes'),'String'));

Mxdims=size(VObj.Rho);
[VMmg.xgrid,VMmg.ygrid,VMmg.zgrid]=meshgrid((-(Mxdims(2)-1)/2)*VObj.XDimRes:VMmg.xdimres:((Mxdims(2)-1)/2)*VObj.XDimRes,...
    (-(Mxdims(1)-1)/2)*VObj.YDimRes:VMmg.ydimres:((Mxdims(1)-1)/2)*VObj.YDimRes,...
    (-(Mxdims(3)-1)/2)*VObj.ZDimRes:VMmg.zdimres:((Mxdims(3)-1)/2)*VObj.ZDimRes);

VMmg.axes=handles.MagField_axes;
handles.IV=struct(...
    'xslice',(handles.Simuh.SV.Slice - (Mxdims(2)-1)/2)* VObj.XDimRes,...
    'yslice',(handles.Simuh.CV.Slice - (Mxdims(1)-1)/2)* VObj.YDimRes,...
    'zslice',(handles.Simuh.AV.Slice - (Mxdims(3)-1)/2)* VObj.ZDimRes ...
    );
view(handles.MagField_axes,3);

%Open current loaded Mag XML file
if isfield(handles.Simuh,'MagXMLFile')
    import javax.swing.*
    import javax.swing.tree.*;
    handles.MagStruct=DoParseXML(handles.Simuh.MagXMLFile);
    handles.MagXMLFile=handles.Simuh.MagXMLFile;
    DoWriteXML2m(handles.MagStruct,[handles.MagXMLFile(1:end-3) 'm']);
    %     guidata(hObject, handles);
    Root=DoConvStruct2Tree(handles.MagStruct);
    handles.MagTree=uitree('V0','Root',Root,'SelectionChangeFcn', {@ChkNode,handles});
    set(handles.MagTree.getUIContainer,'Units','normalized');
    set(handles.MagTree.getUIContainer,'Position',[0.205,0.4,0.2,0.58]);
    handles.MagTreeModel=DefaultTreeModel(Root);
    handles.MagTree.setModel(handles.MagTreeModel);
    try
        DoPlotMagFie(handles);
        set(handles.Forwards_pushbutton,'Enable','on');
        set(handles.Backwards_pushbutton,'Enable','on');
        set(handles.Display_pushbutton,'Enable','on');
        handles=guidata(hObject);
    catch me
        error_msg{1,1}='ERROR!!! Displaying dB0 field aborted.';
        error_msg{2,1}=me.message;
        errordlg(error_msg);
    end

end

% Choose default command line output for MagDesignPanel
handles.input=varargin;
handles.output=hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes MagDesignPanel wait for user response (see UIRESUME)
% uiwait(handles.MagDesignPanel_figure);

function ChkMagElemListTreeNode(tree, value, handles)

handles=guidata(handles.MagDesignPanel_figure);
SelNode=tree.SelectedNodes;
SelNode=SelNode(1);
Level=SelNode.getValue;
if ~ischar(Level)
    LevelChg=diff(Level,1);
    Level([LevelChg; -1]==1)=[];
end

Node='handles.MagElemListStruct';
if Level=='0'
    eval(['handles.MagElemNode=' Node ';']);
else
    for i=1:length(Level)
        Node=[Node '.Children(' num2str(Level(i)) ')'];
    end
    eval(['handles.MagElemNode=' Node ';']);
end
guidata(handles.MagDesignPanel_figure, handles);

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

% Get default command line output from handles structure
varargout{1} = hObject;

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

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

import javax.swing.*
import javax.swing.tree.*;

[filename,pathname,filterindex]=uigetfile({'Mag*.xml','XML-files (*.xml)'},'MultiSelect','off');
if filename~=0
    if ~isempty(handles.Attrh1)
        Attrh1name=fieldnames(handles.Attrh1);
        for i=1:length(Attrh1name)
            delete(handles.Attrh1.(Attrh1name{i}));
        end
        handles.Attrh1=[];
    end
    handles.MagStruct=DoParseXML(fullfile(pathname,filename));
    handles.MagXMLFile=fullfile(pathname,filename);
    guidata(hObject, handles);
    Root=DoConvStruct2Tree(handles.MagStruct);
    handles.MagTree=uitree('v0', 'Root',Root,'SelectionChangeFcn', {@ChkNode,handles});
    set(handles.MagTree.getUIContainer,'Units','normalized');
    set(handles.MagTree.getUIContainer,'Position',[0.205,0.4,0.2,0.58]);
    handles.MagTreeModel=DefaultTreeModel(Root);
    handles.MagTree.setModel(handles.MagTreeModel);
    % Add searchig path
    path(path,pathname);
else
    errordlg('No Mag is loaded!');
    return;
end
guidata(hObject, handles);

function ChkNode(tree, value, handles)

handles=guidata(handles.MagDesignPanel_figure);
SelNode=tree.SelectedNodes;
SelNode=SelNode(1);
Level=SelNode.getValue;
if ~ischar(Level)
    LevelChg=diff(Level,1);
    Level([LevelChg; -1]==1)=[];
end

Node='handles.MagStruct';
if Level=='0'
    eval(['handles.MagNode=' Node ';']);
else
    for i=1:length(Level)
        Node=[Node '.Children(' num2str(Level(i)) ')'];
    end
    eval(['handles.MagNode=' Node ';']);
end
handles.MagNodeLvl=Node;
handles.SelNode=SelNode;

if ~isempty(handles.Attrh1)
    Attrh1name=fieldnames(handles.Attrh1);
    for i=1:length(Attrh1name)
        delete(handles.Attrh1.(Attrh1name{i}));
    end
    handles.Attrh1=[];
end

set(handles.MagAttri_uipanel,'Title',['''' char(SelNode.getName) ''' Mag Attribute'],'Unit','normalized');
if SelNode.isLeaf | SelNode.isRoot
    if strcmp(char(SelNode.getName),'Specials')
        DoEditValue(handles,handles.MagAttri_uipanel,handles.MagNode.Attributes,1,[0.25,0.25,0.0,0.1,0.1,0.05]);
    else
        DoEditValue(handles,handles.MagAttri_uipanel,handles.MagNode.Attributes,1,[0.25,0.25,0.0,0.1,0.1,0.05]);
    end
else
    guidata(handles.MagDesignPanel_figure, handles);
end

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

import javax.swing.*
import javax.swing.tree.*;

try
    switch handles.MagNode.Name
        case 'MRiLabMag'
            errordlg('MRiLabMag Root can not be deleted!');
            return;
    end

    eval([handles.MagNodeLvl '=[];']);
    nN=handles.SelNode.getNextSibling;
    while ~isempty(nN)
        Level=nN.getValue;
        Level(end)=[];
        nN.setValue(Level);
        nN=nN.getNextSibling;
    end
    guidata(hObject, handles);
    pause(0.1);

    % remove node
    handles.MagTree.setSelectedNode(handles.SelNode.getParent);
    handles.MagTreeModel.removeNodeFromParent(handles.SelNode);

    DoWriteXML(handles.MagStruct,handles.MagXMLFile);
catch me
    error_msg{1,1}='ERROR!!! Node operation aborted.';
    error_msg{2,1}=me.message;
    errordlg(error_msg);

end

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

import javax.swing.*
import javax.swing.tree.*;

try
    ChkMagElemListTreeNode(handles.MagElemListTree, [], handles); % In case MagElem.xml updated
    handles=guidata(handles.MagDesignPanel_figure);
    switch handles.MagNode.Name
        case {'MRiLabMag'}
            % add node
            SelNode=handles.MagElemListTree.SelectedNodes;
            SelNode=SelNode(1);
            AddedNode=SelNode.clone;
            switch char(AddedNode.getName)
                case {'MagElem'}
                    errordlg([char(AddedNode.getName) ' can not be added!']);
                    return;
            end
            handles.MagTreeModel.insertNodeInto(AddedNode,handles.SelNode,handles.SelNode.getChildCount());
            handles.MagTree.setSelectedNode(AddedNode); % expand to show added child
            handles.MagTree.setSelectedNode(handles.SelNode); % insure additional nodes are added to parent

            nP=AddedNode.getPreviousSibling;
            if ~isempty(nP)
                Level=nP.getValue;
                AddedNode.setValue([Level; Level(end)+1]);
                eval([handles.MagNodeLvl '.Children(' num2str(Level(end)+1) ')=handles.MagElemNode;']);
                DoWriteXML(handles.MagStruct,handles.MagXMLFile);
            else
                nP=AddedNode.getParent;
                Level=nP.getValue;
                AddedNode.setValue([Level; 1]);
                eval([handles.MagNodeLvl '.Children=handles.MagElemNode;']);
                DoWriteXML(handles.MagStruct,handles.MagXMLFile);
            end
        otherwise
            errordlg(['Mag Element can not be added under ' handles.MagNode.Name '!']);
            return;
    end
    guidata(hObject, handles);
catch me
    error_msg{1,1}='ERROR!!! Node operation aborted.';
    error_msg{2,1}=me.message;
    errordlg(error_msg);

end

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

% update Mag XML file
if ~isempty(handles.Attrh1)
    Attrh1name=fieldnames(handles.Attrh1);
    for i=1:2:length(Attrh1name)
        switch get(handles.Attrh1.(Attrh1name{i+1}),'Style')
            case 'edit'
                handles.MagNode.Attributes((i+1)/2).Value=get(handles.Attrh1.(Attrh1name{i+1}),'String');
            case 'checkbox'
                if get(handles.Attrh1.(Attrh1name{i+1}),'Value')
                    handles.MagNode.Attributes((i+1)/2).Value='^1';
                else
                    handles.MagNode.Attributes((i+1)/2).Value='^0';
                end
            case 'popupmenu'
                handles.MagNode.Attributes((i+1)/2).Value(2)=num2str(get(handles.Attrh1.(Attrh1name{i+1}),'Value'));
        end
    end
    eval([handles.MagNodeLvl '=handles.MagNode;']);
    DoWriteXML(handles.MagStruct,handles.MagXMLFile);
    % update associated m function
    DoWriteXML2m(handles.MagStruct,[handles.MagXMLFile(1:end-3) 'm']);
end

% update MagSimuAttr.xml
tabs=get(handles.Setting_tabgroup,'Children');
for j=1:length(handles.MagSimuAttrStruct.Children)
    Attrh2=get(tabs(j),'Children');
    for i=1:2:length(Attrh2)
        if ~iscell(get(Attrh2(end-i),'String'))
            handles.MagSimuAttrStruct.Children(j).Attributes((i+1)/2).Value=get(Attrh2(end-i),'String');
        else
            handles.MagSimuAttrStruct.Children(j).Attributes((i+1)/2).Value(2)=num2str(get(Attrh2(end-i),'Value'));
        end
    end
end
DoWriteXML(handles.MagSimuAttrStruct,[handles.Simuh.MRiLabPath filesep 'Macro' filesep 'MagElem' filesep 'MagSimuAttr.xml']);

guidata(hObject, handles);

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

import javax.swing.*
import javax.swing.tree.*;

try
    switch handles.MagNode.Name
        case 'MRiLabMag'
            errordlg('MRiLabMag Root can not be copied!');
            return;
    end
    handles.CopSelNode=handles.SelNode.clone;
    handles.CopMagNode=handles.MagNode;
    guidata(hObject, handles);
catch me
    error_msg{1,1}='ERROR!!! Node operation aborted.';
    error_msg{2,1}=me.message;
    errordlg(error_msg);

end

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

import javax.swing.*
import javax.swing.tree.*;

try
    switch handles.MagNode.Name
        case {'MRiLabMag'}
            % paste node
            handles.MagTreeModel.insertNodeInto(handles.CopSelNode,handles.SelNode,handles.SelNode.getChildCount());
            handles.MagTree.setSelectedNode(handles.CopSelNode); % expand to show added child
            handles.MagTree.setSelectedNode(handles.SelNode); % insure additional nodes are added to parent

            nP=handles.CopSelNode.getPreviousSibling;
            if ~isempty(nP)
                Level=nP.getValue;
                handles.CopSelNode.setValue([Level; Level(end)+1]);
                eval([handles.MagNodeLvl '.Children(' num2str(Level(end)+1) ')=handles.CopMagNode;']);
                DoWriteXML(handles.MagStruct,handles.MagXMLFile);
            else
                nP=handles.CopSelNode.getParent;
                Level=nP.getValue;
                handles.CopSelNode.setValue([Level; 1]);
                eval([handles.MagNodeLvl '.Children=handles.CopMagNode;']);
                DoWriteXML(handles.MagStruct,handles.MagXMLFile);
            end
        otherwise
            errordlg(['Mag Element can not be added under ' handles.MagNode.Name '!']);
            return;
    end
    handles.CopSelNode=[];
    handles.CopMagNode=[];
    guidata(hObject, handles);
catch me
    error_msg{1,1}='ERROR!!! Node operation aborted.';
    error_msg{2,1}=me.message;
    errordlg(error_msg);

end

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

set(hObject,'Enable','off');
set(hObject,'String','Calc...');
pause(0.01);
try
    DoPlotMagFie(handles);
    set(handles.Forwards_pushbutton,'Enable','on');
    set(handles.Backwards_pushbutton,'Enable','on');
    set(handles.Display_pushbutton,'Enable','on');
catch me
    error_msg{1,1}='Error occurred at displaying dB0 field.';
    error_msg{2,1}=me.message;
    errordlg(error_msg);
end
set(hObject,'Enable','on');
set(hObject,'String','Execute');

% --------------------------------------------------------------------
function SavePanel_uipushtool_ClickedCallback(hObject, eventdata, handles)
% hObject    handle to SavePanel_uipushtool (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
DoSaveSnapshot(handles.MagDesignPanel_figure);

% --- Executes when user attempts to close MagDesignPanel_figure.
function MagDesignPanel_figure_CloseRequestFcn(hObject, eventdata, handles)
% hObject    handle to MagDesignPanel_figure (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: delete(hObject) closes the figure
delete(hObject);

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

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

global VObj;
TAttr=get(handles.Attrh2.('Plane'),'String');
switch TAttr{get(handles.Attrh2.('Plane'),'Value')}
    case 'XY'
        handles.IV.zslice=handles.IV.zslice+VObj.ZDimRes;
    case 'XZ'
        handles.IV.yslice=handles.IV.yslice+VObj.YDimRes;
    case 'YZ'
        handles.IV.xslice=handles.IV.xslice+VObj.XDimRes;
end
DoUpdateSlice(handles.MagField_axes,handles.dB0,handles.IV,'Mag');
guidata(hObject, handles);

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

global VObj;
TAttr=get(handles.Attrh2.('Plane'),'String');
switch TAttr{get(handles.Attrh2.('Plane'),'Value')}
    case 'XY'
        handles.IV.zslice=handles.IV.zslice-VObj.ZDimRes;
    case 'XZ'
        handles.IV.yslice=handles.IV.yslice-VObj.YDimRes;
    case 'YZ'
        handles.IV.xslice=handles.IV.xslice-VObj.XDimRes;
end
DoUpdateSlice(handles.MagField_axes,handles.dB0,handles.IV,'Mag');
guidata(hObject, handles);

% --------------------------------------------------------------------
function Cursor_uitoggletool_ClickedCallback(hObject, eventdata, handles)
% hObject    handle to Cursor_uitoggletool (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
if ~isfield(handles,'dcm')
    handles.dcm=datacursormode(handles.MagDesignPanel_figure);
    datacursormode on
    set(handles.dcm,'updatefcn',{@Cursor_update,handles});
else
    if strcmp(get(handles.dcm,'Enable'),'on')
        set(handles.dcm,'Enable','off');
    else
        set(handles.dcm,'Enable','on');
    end
end
guidata(hObject, handles);

function output_txt = Cursor_update(obj,event_obj,handles)
% Display the position of the data cursor
% obj          Currently not used (empty)
% event_obj    Handle to event object
% output_txt   Data cursor text string (string or cell array of strings).

global VObj;
global VCtl;
global VMmg;

handles=guidata(handles.MagDesignPanel_figure);
pos = get(event_obj,'Position');
pos(1)=round((pos(1)+VCtl.ISO(1)*VObj.XDimRes)/VMmg.xdimres);
pos(2)=round((pos(2)+VCtl.ISO(2)*VObj.YDimRes)/VMmg.ydimres);
pos(3)=round((pos(3)+VCtl.ISO(3)*VObj.ZDimRes)/VMmg.zdimres);

output_txt = {['X: ',num2str(pos(1),4)],...
    ['Y: ',num2str(pos(2),4)],...
    ['Z: ',num2str(pos(3),4)],...
    ['Value: ',num2str(handles.dB0(pos(2),pos(1),pos(3)),4)]};

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

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

prompt = {'Enter NEW mag name:','Enter NEW mag note:'};
dlg_title = 'Input for creating new mag';
num_lines = 1;
def = {'Mag_CustomHead','Customized dB0 Field Head Magnet'};
answer = inputdlg(prompt,dlg_title,num_lines,def);
if isempty(answer)
    warndlg('No new mag file is created!');
    return
end
if isempty(answer{1})
    warndlg('No new mag file is created!');
    return
end
if ~strcmp(answer{1}(1:4),'Mag_')
    warndlg('Mag name must have a prefix ''Mag_'', please rename your mag!');
    return
end
magpath=uigetdir(pwd,'Specify a saving path for new mag file.');
if magpath==0
    warndlg('You have to specify a saving path!');
    return;
end

% Create new mag file
mkdir([magpath filesep answer{1}]); % make folder
copyfile([handles.Simuh.MRiLabPath filesep 'Config' filesep 'Mag' filesep 'Head' filesep 'Mag_LinearHead' filesep 'Mag_LinearHead.xml'],...
    [magpath filesep answer{1} filesep answer{1} '.xml']); % copy mag temple
% Add searchig path
path(path,[magpath filesep answer{1}]);
% Load new mag file
import javax.swing.*;
import javax.swing.tree.*;
if ~isempty(handles.Attrh1)
    Attrh1name=fieldnames(handles.Attrh1);
    for i=1:length(Attrh1name)
        delete(handles.Attrh1.(Attrh1name{i}));
    end
    handles.Attrh1=[];
end
handles.MagStruct=DoParseXML([magpath filesep answer{1} filesep answer{1} '.xml']);
for i=1:length(handles.MagStruct.Attributes)
    if strcmp(handles.MagStruct.Attributes(i).Name,'Name')
        handles.MagStruct.Attributes(i).Value = answer{1};
    elseif strcmp(handles.MagStruct.Attributes(i).Name,'Notes')
        handles.MagStruct.Attributes(i).Value = answer{2};
    end
end
DoWriteXML(handles.MagStruct,[magpath filesep answer{1} filesep answer{1} '.xml']); %update new mag file
handles.MagXMLFile=[magpath filesep answer{1} filesep answer{1} '.xml'];
DoWriteXML2m(handles.MagStruct,[handles.MagXMLFile(1:end-3) 'm']);
guidata(hObject, handles);
Root=DoConvStruct2Tree(handles.MagStruct);
handles.MagTree=uitree('V0','Root',Root,'SelectionChangeFcn', {@ChkNode,handles});
set(handles.MagTree.getUIContainer,'Units','normalized');
set(handles.MagTree.getUIContainer,'Position',[0.205,0.4,0.2,0.58]);
handles.MagTreeModel=DefaultTreeModel(Root);
handles.MagTree.setModel(handles.MagTreeModel);
cla(handles.MagField_axes);
guidata(hObject, handles);

% --- Executes on mouse press over figure background, over a disabled or
% --- inactive control, or over an axes background.
function MagDesignPanel_figure_WindowButtonDownFcn(hObject, eventdata, handles)
% hObject    handle to MagDesignPanel_figure (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

tempaxes = gca;
try
    if strcmp(get(gcf,'Selectiontype'),'alt')
        figure 			% Create a new figure
        ax_new = axes;
        copyaxes(tempaxes, ax_new);
    end
catch me
end

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

assignin('base', 'dB0', handles.dB0);
MU_Matrix_Display('dB0','Magnitude');

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function out = void(f, out)

% out = void(f, out)
%
% This is a helper function to enable more complex anonymous functions.
%
% While some functions don't return anything, it can be useful to return
% something after calling them. This function calls the input function and
% then returns either 1 (the default) or the second input, which can be any
% value of the user's choosing. Consider the following anonymous function:
%
% plot_it = @() {figure(1), ...
%                clf(), ...
%                plot(randn(1, 50), '.'), ...
%                axis('off')};
%
% plot_it() % This throws an error.
%
% Clearly, this should open a figure, clear it, plot some points, and turn
% off the axis. However, the axis('off') command doesn't return anything,
% and something is necessary for that position of the cell array, so MATLAB
% throws an error. However, if we just pass the function to the "void"
% command, as below, it works just fine, returning a 1 that we simply don't
% use but that fills the cell array as MATLAB requires.
%
% plot_it = @() {figure(1), ...
%                clf(), ...
%                plot(randn(1, 50), '.'), ...
%                void(@() axis('off'))};
%
% plot_it() % This works.
%
% Tucker McClure
% Copyright 2013 The MathWorks, Inc.

    if nargin == 1, out = 1; end
    f();

end

author: Sable


text:
function testConvolutionalFilters
% testConvolutionalFilters.m This script compares convolutional filter
% responses as currently implemented to those submitted to IBSI-2 to
% ensure continued compliance.
%
% AI 12/09/22

init_ML_DICOM

%Path to "gold standard" response maps
CERRPath = getCERRPath;
CERRPathSlashes = strfind(getCERRPath,filesep);
topLevelCERRDir = CERRPath(1:CERRPathSlashes(end-1));
stdDir = fullfile(topLevelCERRDir,'Unit_Testing','data_for_cerr_tests',...
    'IBSI2_synthetic_phantoms','Results');

%Create temp output dir
testDir = fullfile(topLevelCERRDir,'Unit_Testing','tests_for_cerr');
tmpDir = fullfile(testDir,'IBSI2');
mkdir(tmpDir)

%Compute filter response maps (IBSI 2 -phase1)
runIBSI2benchmarkFilters(tmpDir,'all');

%Assess deviation from standard
configC = { '1a1','1a2','1a3','1a4','1b1','2a1','2b1','2c1','3a1','3a2',...
    '3a3','3b1','3b2','3b3','3c1','3c2','3c3','4a1','4a2','4b1','4b2',...
    '5a1','5a2','6a1','6a2'};
assertTOL = 1e-5;

disp(['========= Maximum difference for filt config. =========='])
for n = 1:length(configC)
    currResponse3M = niftiread(fullfile(tmpDir,[configC{n},'.nii']));
    stdResponse3M = niftiread(fullfile(stdDir,[configC{n},'.nii']));
    diff3M = currResponse3M - stdResponse3M;
    maxDiff = max(diff3M(:));
    assertAlmostEqual(currResponse3M,stdResponse3M,assertTOL);
    disp([configC{n},': ', sprintf('%0.1g',maxDiff)])

end

rmdir(tmpDir,'s')

disp('========= testConvolutionalFilters succeeded. ==========')

end


text:
function ls = logsum(xx,dim)
% ls = logsum(x,dim)
%
% returns the log of sum of logs
% computes ls = log(sum(exp(x),dim))
% but in a way that tries to avoid underflow/overflow
%
% basic idea: shift before exp and reshift back
% log(sum(exp(x))) = alpha + log(sum(exp(x-alpha)));
%
% This program was originally written by Sam Roweis

if(length(xx(:))==1) ls=xx; return; end

xdims=size(xx);
if(nargin<2)
  dim=find(xdims>1);
end

alpha = max(xx,[],dim)-log(realmax)/2;
repdims=ones(size(xdims)); repdims(dim)=xdims(dim);
ls = alpha+log(sum(exp(xx-repmat(alpha,repdims)),dim));

repo: DeepFi


text:
function varargout = sqrt(varargin)
%SQRT   Square root.
%   SQRT(F) returns the square root of a positive CHEBFUN2 F.

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

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

end


text:
%GETDATA Get data fields in mapping
%
%    DATA = GETDATA(W)
%
% Get content of data field of W
%
%    DATA = GETDATA(W,N)
%
% Get content of N-th cell in W.DATA.
%
%    DATA = GETDATA(W,FIELD)
%
% Get content of desired field if W.DATA is structure

author: marianux

repo: ecg-kit

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


text:
clear all;
addpath(genpath('../CoreModules'));
n_epoch=25;
dataset_name='cifar';
network_name='cnn';
use_gpu=(gpuDeviceCount>0) %use gpu or not
if use_gpu
    %Requires Neural Network Toolbox to use it.
    opts.use_nntoolbox=license('test','neural_network_toolbox')
end
%function handle to prepare your data
PrepareDataFunc=@PrepareData_CIFAR_CNN;
%function handle to initialize the network
NetInit=@net_init_cifar_cnn;

%automatically select learning rates
use_selective_sgd=1;
%select a new learning rate every n epochs
ssgd_search_freq=10;
learning_method=@sgd2;%training method: @sgd,@adagrad,@rmsprop,@adam,@sgd2
opts.parameters.weightDecay=5e-4;
sgd_lr=1.0;%This is unnecessary as lr is automatically determined.

Main_Template(); %call training template

author: yechengxi

repo: LightNet

save_path: github-repos/MATLAB/yechengxi-LightNet


text:
function boolVal=isProbablePrime(n,numRuns)
%%ISPROBABLEPRIME Determine whether a positive integer is prime, with an
%                 algorithm that has a small probablity of falsely
%                 declaring a non-prime to be a prime, but which is faster
%                 than dbrute force division.
%
%INPUTS: n A real positive integer. This value must be <=sqrt(flintmax).
%          The algorithm gets slower for smaller numbers.
%  numRuns An optional parameter specifying how many times the algorithm is
%          rerun when it initially detects a probable prime. The default if
%          omitted or an empty matrix is passed is 25. The more runs, the
%          more conficence one can have in the solution. In Chapter 4.5.4
%          of [1], it is noted that the algorithm might have a pase
%          positive up to 1/4 of the time, so the likelihood of having 25
%          of them is about (1/4)^25=8.8818e-16 is very slow.
%
%OUTPUTS: boolVal This is true if n is a probable prime and false if n is
%                 definitely not a prime.
%
%This function implements Algorithm P in Section 4.5.4 of [1].
%
%EXAMPLE:
%Here, we test a known Mersenne prime and a nearby number that is not
%prime.
% boolVal0=isProbablePrime(2^19-1)
% boolVal1=isProbablePrime(2^19-3)
%One should get boolVal0=true and boolVal1=false.
%
%REFERENCES:
%[1] D. E. Knuth, The Art of Computer Programming. Vol. 2: Seminumerical
%    Algorithms, 3rd Edition, Boston: Addison-Wesley, 2011.
%
%October 2018 David Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(n<0||~isreal(n)||n~=fix(n))
    error('n must be a real non-negative integer.')
end

if(n>sqrt(flintmax))
   error('The algorithm will not work with n>sqrt(flintmax).')
end

if(nargin<2||isempty(numRuns))
    numRuns=25;
end

%If divisible by 2.
if(mod(n,2)==0)
   boolVal=false;
   return
end

%We now want to find an odd q such that n-1=2^k*q
q=n-1;
k=0;
pow2Val=1;
while(mod(q,2)==0)
    pow2Val=pow2Val*2;
    q=q/2;
    k=k+1;
end

for curRun=1:numRuns
    %Step P1
    x=randi(n-2)+1;%1<x<n
    %Step P2
    %This calculation method helps avoid overflow.
    j=0;
    %This does y=mod(x^q,n), but in a manner that avoids numerical
    %overflows unless x is very big.
    y=powModuloC(x,q,n);

    while(1)
        %Step P3
        if(y==n-1||y==1&&j==0)
            boolVal=true;
            break;
        end

        %Step P4
        if(y==1&&j>0)
            boolVal=false;
            return
        end

        j=j+1;
        if(j<k)
            y=mod(y^2,n);
        else
            boolVal=false;
            return;
        end
    end
end
end

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


text:
function alpha = p06_alpha ( )

%*****************************************************************************80
%
%% P06_ALPHA returns ALPHA for problem 6.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 July 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real ALPHA, the value of ALPHA.
%
  alpha = 2.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function scores = iForestAnomaly(X, t, psi)
%iForestAnomaly Use the path length to calculate the anomaly score
%   X: a m by n matrix, each row is a data point
%   t: the number of trees to construct
%   psi: the sub-sampling size

forest = iForest(X, t, psi);

[m ~] = size(X);
scores = zeros(m, 1);
for i=1:m
    x = X(i, :);
    for tree=forest
        path_length = 0;
        treenode = tree;
        while ~strcmp(class(treenode), 'iTreeLeaf')
            if treenode.compare(x)
                treenode = treenode.Right;
            else
                treenode = treenode.Left;
            end
            path_length = path_length + 1;
        end
        scores(i) = scores(i) + path_length + adjustment(treenode.Size);
    end
    scores(i) = 2^(-scores(i) / length(forest) / adjustment(m));
end
end

function value=adjustment(n)
%The average path length of an unsuccessful BST search
    if n<=1
        value = 0;
    else
        value = 2 * (log(n-1)+0.5772156649) - 2*(n-1) / n;
    end
end

author: dsmi-lab-ntust

repo: AnomalyDetectionToolbox


text:
function s = safeStr(s)
% Change punctuation characters to they print properly

s = strrep(s, '\', '/');
s = strrep(s, '_', '-');


text:
function [M0,M1,L,M2] = spm_mfa_bi_multi(S,C)
% Bilinear form for multiple Gibb's ensembles
% FORMAT [M0,M1,L,M2] = spm_mfa_bi_multi(S,C)
%--------------------------------------------------------------------------
% S(s)   - MFA system specification structure for s ensembles
% Required fields:
%     S(i).M:  [1x1 struct]    - dynamic model structure
%     S(i).J0: [n x n double]  - Jacobian
%     S(i).J1: {1 x M.m cell}  - dJ0/du
%     S(i).L : [l x n double]  - d<y>/dp
%     S(i).u:  [n x m double]  - probability modes
%     S(i).v:  [m x n double]  - v*u = 1
%     S(i).X:  [n x d double]  - evaluation points of state [d]-space
%     S(i).x:  {1 x d cell}    - range of state [d]-space
%     S(i).p0: [n x 1 sparse]  - expansion point
%
% C{s x s cell} - coupling cell = dP/d<y> (change in parameters of S(i).M with
%                                   mean outputs of S(j).M - [p x l double])
%
% M0  [ns + 1 x ns + 1double]  - 1st order Bilinear matrix dq/dt;
% M1  {M.m x s}                - 2nd order Bilinear matrix dM0/du
% M2  {s x s}                  - 2nd order Bilinear matrix dM0/dC
% L   [1s x ns + 1 double]     - output matrix          <y> = L*q;
%
% Transformed probability states:  q = [1; v*(p(X) - p0)];
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_mfa_bi_multi.m 4936 2012-09-18 19:47:55Z karl $

% indices
%--------------------------------------------------------------------------
s     = length(S);                      % number of ensembles
m     = length(S(1).J1);                % number of inputs
I     = {(1:size(S(1).v,1)) + 1};       % ensemble indices
for i = 2:s
    I{i} = (1:size(S(i).v,1)) + I{i - 1}(end);
end
q     = I{s}(end);
M0    = sparse(q,q);
M1    = {};
L     = sparse(0,q);

% compute M0 = df/dx. M1 = d(df/dx)/du amd L = dy/dx
%--------------------------------------------------------------------------
for i = 1:s

    % Ist order bilinear operator M0
    %======================================================================
    M0(I{i},1)       = S(i).v*S(i).J0*S(i).p0;
    M0(I{i},I{i})    = S(i).v*S(i).J0*S(i).u;

    % 2nd order bilinear operators M1{i} - dM0/du
    %======================================================================
    for j = 1:m
        M            = sparse(q,q);
        M(I{i},1)    = S(i).v*S(i).J1{j}*S(i).p0;
        M(I{i},I{i}) = S(i).v*S(i).J1{j}*S(i).u;
        M1{m,i}      = M;
    end

    % output matrix <y>
    %======================================================================
    j         = (1:size(S(i).L,1)) + size(L,1);
    L(j,1)    = S(i).L*S(i).p0;
    L(j,I{i}) = S(i).L*S(i).u;

end

% return unless coupling matrices are required
%--------------------------------------------------------------------------
if nargout < 4
    return
end

% 2nd order bilinear operators M2{i} - dM0/dC
%==========================================================================
M2    = cell(s,s);
dP    = 1e-1;
for i = 1:s

    % target ensemble df/dP = dJ/dP*p0
    %----------------------------------------------------------------------
    M     = S(i).M;
    p     = length(M.pE);
    for k = 1:p
        M.pE      = S(i).M.pE;
        M.pE(k)   = M.pE(k) + dP;
        dJdP      = (spm_mfa(M,S(i).x) - S(i).J0)/dP;
        dfdP(:,k) = S(i).v*dJdP*S(i).p0;
    end

    for j = 1:s

        % source ensemble dP/dx = dP/d<y>*d<y>/dx = C*L
        %------------------------------------------------------------------
        dfdx      = sparse(q,q);
        if ~isempty(C{i,j})
            dPdx            = C{i,j}*S(j).L*S(j).u;
            dfdx(I{i},I{j}) = dfdP*dPdx;
        end

        % M2 = df/dx = df/dP*dP/dx
        %------------------------------------------------------------------
        M2{i,j} = dfdx;

    end
end

author: spm