text:
function combo_test21 ( )

%*****************************************************************************80
%
%% COMBO_TEST21 tests NPART_TABLE and PART_TABLE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 January 2011
%
%  Author:
%
%    John Burkardt
%
  maxn = 10;
  maxpart = 5;

  fprintf ( 1, ' \n' );
  fprintf ( 1, 'COMBO_TEST21\n' );
  fprintf ( 1, '  NPART_TABLE tabulates partitions\n' );
  fprintf ( 1, '  of N with NPART parts;\n' );
  fprintf ( 1, '  PART_TABLE tabulates partitions of N.\n' );

  offset = 1;

  p = npart_table ( maxn, maxpart );

  p2 = part_table ( maxn );

  fprintf ( 1, ' \n' );
  fprintf ( 1, '    I P(I)  P(I,0) P(I,1) P(I,2) P(I,3) P(I,4) P(I,5)\n' );
  fprintf ( 1, ' \n' );

  for i = 0 : maxn
    fprintf ( 1, '%5d%5d', i, p2(i+offset) );
    for j = 0 : maxpart
      fprintf ( 1, '%5d', p(i+offset,j+offset) );
    end
    fprintf ( 1, '\n' );
  end

  return
end


text:
function J=uGradient2D(xTar,lRx,M,includeV)
%%UGRADIENT2D Determine the gradient of a direction cosine measurement in
%          2D with respect to 2D position. Relativity and atmospheric
%          effects are not taken into account.
%
%INPUTS: xTar The 2XN set of target position vectors in the global
%          coordinate system with [x;y] components for which a gradient is
%          desired.
%      lRx The 2X1 position vector of the receiver. If omitted, the
%          receiver is placed at the origin.
%        M A 2X2 rotation matrix from the global Coordinate system to the
%          orientation of the coordinate system at the receiver. If
%          omitted, it is assumed to be the identity matrix.
% includeV An optional boolean value indicating whether a second direction
%          cosine component should be included. The u direction cosine is
%          one parts of a 2D unit vector. The default if this parameter is
%          omitted or an empty matrix is passed is false.
%
%OUTPUTS: J A 1X2XN set of N Jacobian matrices (or a 2X2XN set if includeV
%           is true) where the row represents u (and v) and the columns
%           take the partial derivatives of the rows with respect to [x,y]
%           in that order.
%
%A derivation of the components of the Jacobian of a 3D u-v measurement is
%given in [1]. The results here are similar.
%
%EXAMPLE:
%Here, we validate the results using numerical differentiation.
% x=[1000;4];
% lRx=[-40;2];
% M=eye(2,2);
% J=uGradient2D(x,lRx,M);
% f=@(x)getUDirection2D(x,lRx,M);
% N=3;
% epsVal=[1;0.01];
% numJAug=numDiff(x,f,2,N,epsVal);
% JNumDiff=numJAug(2,:);
% RelErr=max(abs((J-JNumDiff)./JNumDiff))
%One will see that J and JNumDiff agree to more than 7 digits, indicating
%good agreeement.
%
%REFERENCES:
%[1] D. F. Crouse, "Basic tracking using nonlinear 3D monostatic and
%    bistatic measurements," IEEE Aerospace and Electronic Systems
%    Magazine, vol. 29, no. 8, Part II, pp. 4-53, Aug. 2014.
%
%April 2017 David F.Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

N=size(xTar,2);

if(nargin<4||isempty(includeV))
    includeV=false;
end

if(nargin<3||isempty(M))
   M=eye(2,2);
end

if(nargin<2||isempty(lRx))
   lRx=zeros(2,1);
end

J=zeros(1+includeV,2,N);
for curPoint=1:N
    %Convert the state into the local coordinate system.
    xLocal=M*(xTar(1:2,curPoint)-lRx(1:2));

    r3=norm(xLocal)^3;
    x=xLocal(1);
    y=xLocal(2);

    %The gradient of u in local coordinates.
    du=zeros(1,2);
    du(1)=y^2/r3;
    du(2)=-x*y/r3;

    J(1,:,curPoint)=du*M;

    if(includeV)
        dv=zeros(1,2);
        dv(1)=-x*y/r3;
        dv(2)=x^2/r3;
        %Now, the gradient vectors for the angular components must be
        %rotated back into the global coordinate system.
        J(2,:,curPoint)=dv*M;
    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.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [ hitsclean, polyidsclean, rayidsclean ] = invalidateintersections( pos_IMUinG, hits, polyids, rayids )
%INVALIDATEINTERSECTIONS Will handle that a single ray hits many planes
%   We will go through each rayid and from there we will check to see which
%   of the given intersections for that ray are closest to the given
%   position of the LIDAR. We know that the closest ray would intersect
%   first so that is the only valid ray intersection.

% Sort the arrays based on the rayid
[rayids, indices] = sort(rayids);
hits = hits(indices,:);
polyids = polyids(indices,:);

% Array to keep track of the max distance from our pos
maxdist = inf(1,max(rayids));
maxids = zeros(1,max(rayids));

% Now we know we are sorted, so lets loop through each one and find the
% shortest distance for the given ray ids
for ii=1:size(rayids,1)
    % compute the square difference between this point and the LIDAR pose
    diff = sqrt(sum((pos_IMUinG-hits(ii,1:3)).^2));
    % next lets see if it is the new min
    if diff < maxdist(rayids(ii,1))
        maxdist(rayids(ii,1)) = diff;
        maxids(rayids(ii,1)) = ii;
    end
end

% our output args
hitsclean = [];
polyidsclean = [];
rayidsclean = [];

% finally lets create the return values
cleanct = 0;
for ii=1:size(maxdist,2)
    % if it is not inf then we have the min
    if maxdist(1,ii) ~= Inf
        hitsclean = [hitsclean; hits(maxids(1,ii),:)];
        polyidsclean = [polyidsclean; polyids(maxids(1,ii),:)];
        rayidsclean = [rayidsclean; rayids(maxids(1,ii),:)];
        cleanct = cleanct + 1;
    end
end

% Debug info
fprintf('RAY: %d left after invalidation.\n',cleanct)

end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
no = 6;
ss = 1;
ns = 2;
sc = sqrt(2);
el = 2;
fb = fbCreate(no,ss,ns,sc,el);
for k = [32 64 128],
  tex = unitex(fb,k);
  tsim = textonsim(fb,tex);
  [tim,tperm] = visTextons(tex,fb);
  save(sprintf('unitex_%.2g_%.2g_%.2g_%.2g_%.2g_%d.mat',...
               no,ss,ns,sc,el,k),...
       'fb','tex','tsim','tim','tperm');
end

author: zouchuhang

repo: LayoutNet


%% set parameters
folderpath = '/home/gyeongsikmoon/workspace//Data/Hand_pose_estimation/HANDS2017/training/images/';
filepath = '/home/gyeongsikmoon/workspace/Data/Hand_pose_estimation/HANDS2017/training/Training_Annotation.txt';
frameNum = 957032;

%folderpath = '/home/gyeongsikmoon/workspace//Data/Hand_pose_estimation/HANDS2017/frame/images/';
%filepath = '/home/gyeongsikmoon/workspace/Data/Hand_pose_estimation/HANDS2017/frame/BoundingBox.txt';
%frameNum = 295510;

fp = fopen(filepath);
fid = 1;

tline = fgetl(fp);
while fid <= frameNum

    disp(fid);

    splitted = strsplit(tline);
    img_name = splitted{1};

    if exist(strcat(folderpath,img_name), 'file')
        img = imread(strcat(folderpath,img_name));

        fp_save = fopen(strcat(folderpath,img_name(1:size(img_name,2)-3),'bin'),'w');
        fwrite(fp_save,permute(img,[2,1,3]),'float');
        fclose(fp_save);

        %delete(strcat(folderpath,img_name));
    end

    tline = fgetl(fp);
    fid = fid + 1;
end

fclose(fp);


text:
function R = quat2rot(h)

% Uso: R = quat2rot(h)
%
% costruisce la matrice di rotazione R
% a partire dal quaternione h corrispondente
%
% From Luca Carlone's
% https://bitbucket.org/lucacarlone/pgo3d-duality-opencode/src/ebb6e1b8cebaad7f2aaf581b1d0c0bad737faebb/lib/quat2rot.m?at=master&fileviewer=file-view-default

s(4)=h(1);
s(1)=h(2);
s(2)=h(3);
s(3)=h(4);

R(1,1)=s(1)^2-s(2)^2-s(3)^2+s(4)^2;
R(1,2)=2*(s(1)*s(2)-s(3)*s(4));
R(1,3)=2*(s(1)*s(3)+s(2)*s(4));

R(2,1)=2*(s(1)*s(2)+s(3)*s(4));
R(2,2)=-s(1)^2+s(2)^2-s(3)^2+s(4)^2;
R(2,3)=2*(s(2)*s(3)-s(1)*s(4));

R(3,1)=2*(s(1)*s(3)-s(2)*s(4));
R(3,2)=2*(s(2)*s(3)+s(1)*s(4));
R(3,3)=-s(1)^2-s(2)^2+s(3)^2+s(4)^2;

author: uzh-rpg


text:
 function ob = fatrix2_vertcat(blocks, varargin)
%function ob = fatrix2_vertcat(blocks, [options])
%| vertcat: B = [A1; A2; ...]
%|
%| in
%|	blocks	{cell}	cell array of the (fatrix2) blocks
%|			all blocks must have same idim and imask
%|
%| option
%|	'dim_cat'	dimension along which to concatenate odims
%|			all other odims must match
%|			default is ndim+1 if *all* dimensions match
%|			and imask is not full 1D.
%|			otherwise default is 1
%| out
%|	ob = [A1; A2; ...]
%|
%| If all input blocks have the same odim and have a handle_forw_block
%| then output could also have a handle_forw_block (todo)

arg.dim_cat = []; % default is to infer from odims
arg = vararg_pair(arg, varargin);
arg.out_1d = false;

MM = numel(blocks);

fatrix2_vertcat_check_all_fatrix2(blocks)
fatrix2_vertcat_check_idim(blocks)

[odims_same, odims] = fatrix2_same_odims(blocks);
if ~odims_same % columnize
	warn 'mismatched odims so reverting to single column output'
	arg.out_1d = true;
	for mm=1:MM
		odims{mm} = prod(odims{mm});
	end
end

arg.blocks = blocks;

if isempty(arg.dim_cat)
	b1 = blocks{1};
	if numel(b1.idim) == 1 && size(b1,2) == b1.idim % full 1d input
		arg.dim_cat = 1;
	elseif odims_same
		arg.dim_cat = 'next';
	else
		arg.dim_cat = 1;
	end
end

if streq(arg.dim_cat, 'next')
	if ~odims_same
		fail 'dim_cat ''next'' works only if all odims same'
	end
	odim = [odims{1} MM];
	arg.dim_cat = numel(odim);
else
	id_check = 1:numel(odims{1});
	id_check(arg.dim_cat) = []; % odim must match except dim_cat
	for mm=1:MM
		if ~isequal(odims{mm}(id_check), odims{1}(id_check))
			fail('all blocks must have same odim (except dim_cat)')
		end
	end

	tmp = cell2mat(odims); % [MM ndim]
	tmp = sum(tmp(:,arg.dim_cat));
	odim = odims{1}; odim(arg.dim_cat) = tmp;
end

if fatrix2_vertcat_omask_all_empty(blocks)
	omask = []; % all empty omasks
else
	tmp = cell(MM,1);
	for mm=1:MM
		tmp{mm} = fatrix2_omask_array(blocks{mm});
		if arg.out_1d
			tmp{mm} = tmp{mm}(:);
		end
	end
	omask = cat(arg.dim_cat, tmp{:});
end

arg.mat2cell_arg = fatrix2_vertcat_mat2cell_arg(odims, arg.dim_cat, odim);

% build object
ob = fatrix2('arg', arg, 'caller', 'fatrix2_vertcat', ...
	'idim', blocks{1}.idim, 'imask', blocks{1}.imask, ...
	'odim', odim, 'omask', omask, ...
	'sparse', @fatrix2_vertcat_sparse, ...
	'abs', @fatrix2_vertcat_abs, 'power', @fatrix2_vertcat_power, ...
	'forw', @fatrix2_vertcat_forw, 'back', @fatrix2_vertcat_back);

% fatrix2_vertcat_check_all_fatrix2()
function fatrix2_vertcat_check_all_fatrix2(blocks)
if ~iscell(blocks)
	fail('blocks must be cell array, not "%s"', class(blocks))
end
for mm=1:numel(blocks);
	if ~isa(blocks{mm}, 'fatrix2')
		fail('vertcat requires all fatrix2; use Gmatrix if needed')
	end
end

% fatrix2_vertcat_check_idim()
% check compatibility of input dimensions
function fatrix2_vertcat_check_idim(blocks)
b1 = blocks{1};
for mm=2:numel(blocks);
	bm = blocks{mm};
	if size(bm,2) ~= size(b1,2)
		fail('all blocks must have same #cols')
	end
	if ~isequal(bm.idim, b1.idim)
		fail('all blocks must have same idim')
	end
	if ~isequal(bm.imask, b1.imask)
		fail('all blocks must have same imask')
	end
end

% fatrix2_vertcat_forw(): y = A * x
function y = fatrix2_vertcat_forw(arg, x)

MM = numel(arg.blocks);
yy = cell(MM,1);
for mm=1:MM
	bi = arg.blocks{mm};
	yy{mm} = fatrix2_do_forw(bi, x);
%	yy{mm} = bi.handle_forw(bi.arg, x);
%	if bi.scale ~= 1
%		yy{mm} = bi.scale * yy{mm};
%	end
	if arg.out_1d
		yy{mm} = yy{mm}(:);
	end
end

y = cat(arg.dim_cat, yy{:});

% fatrix2_vertcat_back(): x = A' * y
function x = fatrix2_vertcat_back(arg, y)

yy = mat2cell(y, arg.mat2cell_arg{:});

MM = length(arg.blocks);
for mm=1:MM
	bi = arg.blocks{mm};
	if arg.out_1d
		yy{mm} = reshape(yy{mm}, [bi.odim 1]);
	end
	tmp = fatrix2_do_back(bi, yy{mm});
%	tmp = bi.handle_back(bi.arg, yy{mm});
%	if bi.scale ~= 1
%		tmp = conj(bi.scale) * tmp;
%	end
	if mm == 1
		x = tmp;
	else
		x = x + tmp;
	end
end

% fatrix2_vertcat_abs(): abs(ob)
function ob = fatrix2_vertcat_abs(ob)

MM = numel(ob.arg.blocks);
for mm=1:MM
	ob.arg.blocks{mm} = abs(ob.arg.blocks{mm});
end

% fatrix2_vertcat_power(): ob.^sup
function ob = fatrix2_vertcat_power(ob, sup)

MM = numel(ob.arg.blocks);
for mm=1:MM
	ob.arg.blocks{mm} = (ob.arg.blocks{mm}) .^ sup;
end

% fatrix2_vertcat_mat2cell_arg()
function mat2cell_arg = fatrix2_vertcat_mat2cell_arg(odims, dim_cat, odim)

MM = numel(odims); % odims is cell(MM,1)

if dim_cat == numel(odims{1}) + 1
	for id = 1:numel(odims{1})
		mat2cell_arg{id} = odim(id);
	end
	mat2cell_arg{dim_cat} = ones(MM,1);
	tmp1 = 1:prod(odim);
	tmp1 = reshape(tmp1, odim);
	tmp2 = mat2cell(tmp1, mat2cell_arg{:});
else
	for mm=1:MM
		odim_m = odims{mm};
		if mm==1
			for id = 1:numel(odim_m)
				if id ~= dim_cat
					mat2cell_arg{id} = odim_m(id);
				end
			end
		end
		mat2cell_arg{dim_cat}(mm) = odim_m(dim_cat);
	end
end

% fatrix2_vertcat_omask_all_empty()
% see if all blocks have empty omask (often this will be the case)
function allempty = fatrix2_vertcat_omask_all_empty(blocks)
allempty = true;
for mm=1:numel(blocks)
	allempty = allempty & isempty(blocks{mm}.omask);
end

% fatrix2_vertcat_sparse(): S = sparse(A)
function S = fatrix2_vertcat_sparse(ob)
arg = ob.arg;
S = sparse([]);
MM = numel(arg.blocks);
S = cell(MM,1);
for mm=1:MM
	S{mm} = sparse(arg.blocks{mm});
end
S = cat(1, S{:});

%{
% fatrix2_vertcat_gram()
% todo: this would be useful if W is empty or (block) diagonal
% because it would save memory, but it is extra complexity
function [T, reuse] = fatrix2_vertcat_gram(ob, W, reuse, varargin)

blocks = ob.arg.blocks;
T = cell(size(blocks));
for mm=1:length(blocks)
	A = blocks{mm};
	if isnumeric(A)
		if isempty(W)
			T{mm} = A' * A;
		else
			warn 'todo: this may not work, need piece of W!'
			T{mm} = A' * W * A;
		end
	else
		if isempty(W)
			T{mm} = build_gram(A, [], reuse, varargin{:});
		else
			if isvar('W.arg.blocks{mm}')
				T{mm} = build_gram(A, W.arg.blocks{mm}, ...
					reuse, varargin{:});
			else
				fail('fatrix2_vertcat_gram needs block diag W')
			end
		end
	end
end
T = fatrix_block_sum(T{:});
%}

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
% [ wave, f, t, coh, phases, raw, coi, scale, period, scalef ] = getWavelet( x, Fs, fMin, fMax, nbins, graphics )
%
% x:
% vector - wavelet analysis is applied
% multi-column matrix - wavelet analysis is applied to each column
% separately
%
% two-column matrix - treated as two time series and the csd, coherence is
% computed
% 3D array, 3rd dimension is 2 - treated as multiple samples from two time
% series, coherence etc is computed
%
% Fs - sampling frequency
% fMin, fMax, nbins - parameters for spectral analysis - min/max freqnecy
% (Hz) and number of frequency bins (Morlet is always used so this is
% frequency per se and not just scale)
% graphcis - flag 0/1
%
% does: compute the CWT of each segment, then if two signals also compute
% the cross-spectrum/coherence/phase lag. then average over all segments
% and plot
%
% output:
% wave - the PSD/CSD
% f, t - vectors
% coh - only for 2-channel input
% phases - either for the 1-channel or the phase-difference (from the CSD,
%       not the smoothed estimate)
%
%
% call: wavelet, smoothwavelet, phaseplot, myjet, colormaps

% 08-oct-12 ES

% revisions
% 15-nov-12 (1) added single-channel phase; cross-spectrum, coherence, and
%               phase difference estimates
%           (2) local plotting of spectrogram + phasogram (single-channel)
%               or coherogram + phase differences (two-channel)

% to do: organize input / output handling better, plotting etc
%           also compute mean spectra/coherence/phase
% also - external coherence computation for multiple segments (i.e. average
%       across trials)

%function [ wave, f, t, phases, coi, scale, raw ] = getWavelet( x, Fs, fMin, fMax, nbins )
function [ wave, f, t, coh, phases, raw, coi, scale, period, scalef ] = getWavelet( x, Fs, fMin, fMax, nbins, scaling, graphics )

if isa( x, 'int16' ), x = single( x ); end

if nargin < 6 || isempty( scaling )
    scaling = 'var';
end
if nargin < 7 || isempty( graphics )
    if nargout == 0
        graphics = 1;
    else
        graphics = 0;
    end
end

dt = 1 / Fs;
s0 = 1 / fMax;
tMax = 1 / fMin;
dj = log2( tMax/s0 ) / nbins;

mother = 'MORLET';
%k0 = 6;
%fourier_factor = (4*pi)/(k0 + sqrt(2 + k0^2)); % scale->frequency

wave = [];
f = [];
t = [];
coh = [];
phases = [];
raw = [];
coi = [];
scale = [];
period = [];
scalef = [];

er = 0;
sx = size( x );
switch ndims( x )
    case 2
        switch min( sx )
            case 1
                nsignals = 1;
                nsegments = 1;
                nsamples = max( sx );
                x = x( : );
            case 2
                nsignals = 2;
                nsegments = 1;
                if size( x, 1 ) == 2
                    x = x';
                end
                y = x( :, 2 );
                x = x( :, 1 );
                nsamples = max( sx );
            otherwise
                nsignals = 1;
                nsegments = sx( 2 );
                nsamples = sx( 1 );
        end
    case 3
        if size( x, 3 ) ~= 2
            er = 1;
        else
            nsignals = 2;
            nsegments = sx( 2 );
            nsamples = sx( 1 );
            y = x( :, :, 2 );
            x = x( :, :, 1 );
        end
    otherwise
        er = 1;
end
if er
    error( 'not supported' )
end

% parameters
if isa( scaling, 'char' ) && strcmp( scaling, 'var' )
    switch nsignals
        case 1
            scalef = var( x ); % can be a different number for each segment
        case 2
            scalef = std( x ) .* std( y );
    end
elseif isa( scaling, 'double' ) && size( scaling, 1 ) == ( nbins + 1 )
    scalef = scaling; % can be a different number for each frequency
    scaling = 'z';
else
    scalef = ones( 1, nsegments );
    scaling = 'none';
end
flipidx = ( nbins + 1 ) : -1 : 1;
t = ( 1 : nsamples )' / Fs;

% actually comptue
%[wave,period,scale,coi] = wavelet(Y,dt,pad,dj,s0,J1,mother,param);
switch nsignals
    case 1
        xw = zeros( nbins + 1, nsamples, nsegments );
        for i = 1 : nsegments
            [ xw( :, :, i ), period, scale, coi ] = wavelet( x( :, i ), dt, 1, dj, s0, nbins, mother );
        end
        xw = xw( flipidx, :, : ); % freq, time, segments
        wave = abs( xw ) .^ 2;
        phases = angle( xw );
        raw = xw;
    case 2
        xw = zeros( nbins + 1, nsamples, nsegments );
        yw = xw;
        coh = xw;
        for i = 1 : nsegments
            [ xw( :, :, i ), period, scale, coi ] = wavelet( x( :, i ), dt, 1, dj, s0, nbins, mother );
            [ yw( :, :, i ) ] = wavelet( y( :, i ), dt, 1, dj, s0, nbins, mother );
            % coherence (copied as is from wtc.m):
            sinv=1./(scale');
            X = xw( :, :, i );
            Y = yw( :, :, i );
            wxy = X .* conj( Y ); % complex, single trial
            sX=smoothwavelet(sinv(:,ones(1,nsamples)).*(abs(X).^2),dt,period,dj,scale);
            sY=smoothwavelet(sinv(:,ones(1,nsamples)).*(abs(Y).^2),dt,period,dj,scale);
            sWxy=smoothwavelet(sinv(:,ones(1,nsamples)).*wxy,dt,period,dj,scale);
            Rsq=abs(sWxy).^2./(sX.*sY);
            %phases( :, :, i ) = angle( sWxy );
            coh( :, :, i ) = Rsq( flipidx, : );
            %coh = abs( yo( :, 1, 2 ) .^ 2 ) ./ ( yo( :, 1, 1 ) .* yo( :, 2, 2 ) );

%             subplot( 4, 2, 1 ), [ c h ] = contourf( t, f, log2( abs( flipud( X ) ).^2 ), 100 ); set( h, 'linestyle', 'none' );
%             subplot( 4, 2, 3 ), [ c h ] = contourf( t, f, log2( abs( flipud( Y ) ).^2 ), 100 ); set( h, 'linestyle', 'none' );
%             subplot( 4, 2, 5 ), [ c h ] = contourf( t, f, log2( abs( flipud( wxy ) ).^2 ), 100 ); set( h, 'linestyle', 'none' );
%             subplot( 4, 2, 2 ), [ c h ] = contourf( t, f, log2( flipud( sX ) ), 100 ); set( h, 'linestyle', 'none' );
%             subplot( 4, 2, 4 ), [ c h ] = contourf( t, f, log2( flipud( sY ) ), 100 ); set( h, 'linestyle', 'none' );
%             subplot( 4, 2, 6 ), [ c h ] = contourf( t, f, log2( flipud( sWxy ) ), 100 ); set( h, 'linestyle', 'none' );
%             subplot( 4, 2, 7 ), [ c h ] = contourf( t, f, flipud( abs(wxy).^2./(X.*Y) ), 100 ); set( h, 'linestyle', 'none' );
%             subplot( 4, 2, 8 ), [ c h ] = contourf( t, f, flipud( abs(sWxy).^2./(sX.*sY) ), 100 ); set( h, 'linestyle', 'none' );
%
%             x0 = (abs(X).^2);
%             fmat = sinv(:,ones(1,nsamples));
%             sX=smoothwavelet(fmat.*x0,dt,period,dj,scale);

        end

        % individual channels
        xw = xw( flipidx, :, : );
        yw = yw( flipidx, :, : );
        raw( :, :, :, 1 ) = xw;
        raw( :, :, :, 2 ) = yw;
        %xwave = abs( xw ) .^ 2;
        %ywave = abs( yw ) .^ 2;
        %xphases = angle( xw );
        %yphases = angle( yw );
        % cross spectrum
        xyw = xw .* conj( yw );
        phases = angle( xyw ); % from the CSD
        wave = abs( xyw );
        % for multiple trials, one can also compute the trial-averaged coherence and phase lag by:
        %cohTA = mean( abs( xyw ) .^ 2, 3 ) ./ ( mean( abs( xw ) .^ 2, 3 ) .* mean( abs( yw ) .^ 2, 3 ) );
        %phasesTA = mod( atan2( mean( sin( phases ), 3 ), mean( cos( phases ), 3 ) ), 2 * pi );

end

%xw = flipud( xw );
f = 1 ./ period( flipidx );
coi = 1 ./ coi; % minimum freq to consider at each time point
%f = fliplr( 1 ./ period );
%f = fliplr( 1 ./ scale );
%t = ( 1 : length( x ) )' / Fs;

if graphics

    % here the scaling is by the signal variance

    figure

    if nsignals == 1
        nplots = 1; % PSD
    else
        nplots = 4; % [ PSD1 CSD; COH PSD2 ]
    end

    for np = 1 : nplots
    end

    % scale
    %scalef = mean( scalef );
    %scaleres = 0.25;
    scaleres = 10;
    scalename = '{\sigma}^2';
    pow = zeros( size( wave ) );
    switch scaling
        case { 'var', 'none' }
            for i = 1 : nsegments
                pow( :, :, i ) = log2( abs( wave( :, :, i ) / scalef( i ) ) );
            end
        case 'z'

            for i = 1 : length( f )
                pow( i, :, : ) = ( wave( i, :, : ) - scalef( i, 1 ) ) / scalef( i, 2 );
            end
    end
    pow = mean( pow, 3 ); % average over segments (1 signal)
    levels = min( pow(:) ) : scaleres : max( pow( : ) );
    mphases = mod( atan2( mean( sin( phases ), 3 ), mean( cos( phases ), 3 ) ), 2 * pi );
    %mp = []; for i = 1 : size( phases, 1 ), mp( :, i ) = circ_mean( squeeze( phases( i, :, : ) )' ); end

    % plot
    h1 = subplot( 1, 1, 1 ); %subplot( 2, 1, 1 );
    %[ c h ] = contourf( t, f, pow, levels );
    [ c h ] = contourf( t, f, pow, 100 );
    set( h, 'linestyle','none')
    %xlabel( 'Time (sec)' )
    ylabel( 'Frequency (Hz)' )
    title( sprintf( '%d segments, %d signals', nsegments, nsignals ) )

    % center color limits around log2(1)=0
    if strcmp( scaling, 'var' )
        clim=get(gca,'clim');
        clim=[-1 1]*max(clim(2),3);
        set(gca,'clim',clim)
    end

    % add the cone of influence
    line( t, coi, 'color', [ 0 0 0 ] );
    hold on
    tt=[t([1 1])-dt*.5;t;t([end end])+dt*.5];
    hcoi=fill(tt,1./[period([end 1]) 1./coi period([1 end])],'w');
    %hcoi=fill(tt,[f([end 1]) coi f([1 end])],'w');
    set(hcoi,'alphadatamapping','direct','facealpha',.5)
    hold off

    set( h1, 'box', 'off', 'tickdir', 'out' )

    % add phase arrows (copied as is from xwt.m)
    if nsignals == 2
        Args.ArrowDensity = [30 30];
        Args.ArrowSize = 1;
        Args.ArrowHeadSize = 1;
        ad=mean(Args.ArrowDensity);
        Args.ArrowSize=Args.ArrowSize*30*.03/ad;
        Args.ArrowHeadSize=Args.ArrowHeadSize*Args.ArrowSize*220;
        phs_dt=round(length(t)/Args.ArrowDensity(1));
        tidx=max(floor(phs_dt/2),1):phs_dt:length(t);
        phs_dp=round(length(period)/Args.ArrowDensity(2));
        pidx=fliplr( max(floor(phs_dp/2),1):phs_dp:length(period) );
        phaseplot(t(tidx),f(pidx),2*pi-mphases(pidx,tidx),Args.ArrowSize,Args.ArrowHeadSize);
    end

    % add colorbar
    h = colorbar;
    subplot( h )
    barylbls=rats(2.^(get(h,'ytick')'));
    %barylbls([1 end],:)=' ';
    barylbls(:,all(barylbls==' ',1))=[];
    set(h,'yticklabel',barylbls);
    title( scalename )
    set( h, 'box', 'off', 'tickdir', 'out' )

    colormap( h1, myjet )

    if 1
        figure, %h2 = subplot( 2, 1, 2 );
        switch nsignals
            case 1
                % plot phases separately
                [ c h ] = contourf( t, f, mphases, 10 );
                set( h, 'linestyle','none')
                xlabel( 'Time (sec)' )
                ylabel( 'Frequency (Hz)' )
                %colormap( h2, colormaps( myjet ) )
                set( gca, 'clim', [ 0 2*pi ] )

                % add the cone of influence
                lh = line( t, coi, 'color', [ 0 0 0 ] );
                hold on
                tt=[t([1 1])-dt*.5;t;t([end end])+dt*.5];
                hcoi=fill(tt,1./[period([end 1]) 1./coi period([1 end])],'w');
                set(hcoi,'alphadatamapping','direct','facealpha',.5)
                hold off

                h = colorbar;
                subplot( h )
                title( 'Phase (rad)' )
                set( h, 'box', 'off', 'tickdir', 'out' )
                colormap( colormaps( myjet ) )

            case 2
                % also plot coherence (without phase arrows)
                % plot
                %h2 = subplot( 2, 1, 2 );
                [ c h ] = contourf( t, f, mean( coh, 3 ), 100 );
                %[ c h ] = contourf( t, f, cohTA, 100 );
                set( h, 'linestyle','none')
                xlabel( 'Time (sec)' )
                ylabel( 'Frequency (Hz)' )
                title( sprintf( '%d segments, %d signals', nsegments, nsignals ) )
                set( gca, 'clim', [ 0 1 ] )
                set( h1, 'box', 'off', 'tickdir', 'out' )

                % add the cone of influence
                lh = line( t, coi, 'color', [ 0 0 0 ] );
                hold on
                tt=[t([1 1])-dt*.5;t;t([end end])+dt*.5];
                hcoi=fill(tt,1./[period([end 1]) 1./coi period([1 end])],'w');
                set(hcoi,'alphadatamapping','direct','facealpha',.5)
                hold off

                % add phase plots (only for high coherence values inside the coi)
                aaa=2*pi-mphases;
                aaa(mean( coh, 3 )<.5)=NaN;
                aaa( bsxfun( @lt, f' * ones( 1, nsamples ), coi ) ) = NaN;
                phaseplot(t(tidx),f(pidx),aaa(pidx,tidx),Args.ArrowSize,Args.ArrowHeadSize);

                % add colorbar
                h = colorbar;
                subplot( h )
                title( 'Coherence' )
                set( h, 'box', 'off', 'tickdir', 'out' )
                colormap( myjet )
        end

    end

    %figure, xwt( [t x( :, 1 )],[t y( :, 1 ) ], 'Pad', 1, 'Dj', dj, 'S0', s0, 'J1', nbins );
    %figure, [ Rsq,aWxy ] = waveletCoherence( [t x( :, 1 )],[t y( :, 1 ) ], 'Pad', 1, 'Dj', dj, 'S0', s0, 'J1', nbins, 'mcc', 0, 'MakeFigure', 1 ); title( 'COH' )
    %figure, [ Wxy ] = xwt( [t x( :, i )],[t y( :, i ) ], 'Pad', 1, 'Dj', log2( fMax/fMin ) / nbins, 'S0', 1/fMax, 'MaxScale', 1/fMin, 'MakeFigure', 1 ); title( 'CSD' )

end

%t = [ 0 : 1 : length( x ) - 1 ]' / Fs;
%[ rawwave, period, scale, coi sig95 ] = wt( [ t x ], 'Pad', 1, 'dj', dj, 's0', s0, 'j1', nbins, 'mother', mother, 'MakeFigure', 1 );

return

% notes:
% (1) it is clear what the frequencies are (nbins, log-spaced between fMin and
% fMax), but the amp is unclear to me           DONE
% (2) should use a similar approach to filter the signal at various
% frequency ranges and calculate the spiking rate frequency/phase maps

% to do:
% (1) adjust the scale properly                 DONE
% (2) get the multi-segment version working; basically use the same call to
% wavelet.m, but concatenate the segments (with intervening portions), then
% call once with the first segment to get the proper coi.       DONE
% (3) get the plotting of phases working on this diagram        DONE
% (4) 2-signal version: compute the coherence as in Torrence and Compo
% (smoothing in time- and frequency-domains each of the spectra and the
% cross-spectrum) and add the phase                             DONE

% OK now (15nov12), but not super elegant:
% should partition into two function - getWavelet and plotWavelet,
% the first should only compute, the second should plot with options:
% single signal: just power,power + phase plots.
%       e.g. [ 1 1 ] would make two plots, whereas [ 1 0 ] just on
% two signals: power/phase for each signal; csd/coh/phase for the joint
%       i.e. if only power, plots 2x2 - [ s1 s12_csd; s12_coh; s2 ]
%       if also phase, cuts each plot into two and adds the phase

getWavelet( [ x0( :, 1 ) y0( :, 1 ) ], Fs, Fmin, Fmax, nBins );
getWavelet( [ x0( :, 1 ) ], Fs, Fmin, Fmax, nBins );
getWavelet( [ x0( :, 1 : 10 ) ], Fs, Fmin, Fmax, nBins );
xy = [];
xy( :, :, 1 ) = x0( :, 1 : 10 );
xy( :, :, 2 ) = y0( :, 1 : 10 );
getWavelet( xy, Fs, Fmin, Fmax, nBins );

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [C,S]=fitSpherHarmonics(M,points,vals,a,c)
%%FITSPHERHARMONICS Given a number of field values at a diversity
%               of points in spherical coordinates (either range, azimuth
%               and elevation, or just azimuth and elevation), fit fully
%               normalized spherical harmonic coefficients up to a
%               particular degree to the data. The data being fitted can be
%               real or complex. One can fit multiple sets of data at once
%               as long as they share the same set of points.
%
%INPUTS: M The maximum degree of the coefficients to use.
%   points The 2XN or 3XN set of points in spherical coordinates in the
%          form [r;azimuth;elevation] or just [azimuth;elevation]. Azimuth
%          is measured counterclockwise from the x-axis in the x-y plane.
%          Elevation is measured up from the x-y plane (towards the
%          z-axis). Angles must be given in radians. Note the system to be
%          observable, N>=(M+1)^2.
%     vals The NXnumSets sets of real or complex values at the given
%          points.
%        a The real or complex numerator in the (a/r)^n term in the
%          spherical harmonic sum. If this parameter is omitted or an empty
%          matrix is passed, then a=1 is used.
%        c The real or complex constant value by which the spherical
%          harmonic series is multiplied. If this parameter is omitted or
%          an empty matrix is passed, then c=1 is used.
%
%OUTPUTS: C A numCoeffsXnumSets collections of arrays, one for each set of
%           values passed, holding the coefficient terms that are
%           multiplied by cosines in the spherical harmonic expansion. If
%           given to a CountingClusterSet class, C(n+1,m+1) is the
%           coefficient of degree n and order m. This and S can be given to
%           the function spherHarmonicEval.
%         S A numCoeffsXnumSets collection of arrays, one for each set of
%           values passed, holding the coefficient terms that are
%           multiplied by sines in the spherical harmonic expansion. S has
%           the same structure as C.
%
%The spherical harmonic series being fit is assumed to be of the form
%V=(c/r)*sum_{n=0}^M\sum_{m=0}^n(a/r)^n*(C(n+1,m+1)*cos(m*lambda)+S(n+1,m+1)*sin(m*lambda))*\bar{P}_{nm}(sin(theta))
%where \bar{P}_{nm} is a fully-normalized associated Legendre polynomial.
%This function uses the LegendreCos function to evaluate the Legendre
%polynomials and it them constructs a linear system to solve for the values
%of C and S. Unobservable values of S are where m=0, so S(n+1,0+1)=0 and
%those components are not part of the estimation problem.
%
%Real spherical harmonic models are often used for gravitational, magnetic
%field and some terrain modeling, as discussed in part of [1], whereas
%complex spherical harmonic models are often used for electromagnetic
%modeling.
%
%The complex spherical harmonic model that is usually used in
%electromagnetics is
%V=(c/r)*sum_{n=0}^M\sum_{m=-n}^n(a/r)^n*A(n+1,m+1)*\bar{P}_{nm}(sin(theta))*exp(1j*m*lambda)
%Initially, it just looks like one has combined the S and C values into one
%complex coefficient A and then multiplied the sin term by 1j. However, if
%that were all that had been done, one wouldn't be able to match many
%models, because if A(n+1,m+1)=AR+1j*AI, expanding the term
%A(n+1,m+1)*exp(1j*m*lambda), one gets
%aR*cos(m*lambda)-aI*sin(m*lambda)+1j*(aI*cos(m*lambda)+aR*sin(m*lambda))
%Where one can see that it is not possible to independently set the real
%and the imaginary terms. However, note that the formulation of the inner
%sum used for complex values starts at -m and not at 0. This extension of
%the sum allows one to uniquely form real and complex terms. Though such a
%symmetry does not occur at m=0, that does not actually matter, since the
%sine terms are zero, so one can still independently set the real and the
%imaginary coefficients parts. Indeed, it can be shown that the resulting
%expression allows one to fit the real and complex parts independently if a
%and c are real.
%
%EXAMPLE:
%We take the coefficients for the World Magnetic Model, evaluate the
%magnetic field values at a number of points around the globe, use this
%function to reconstruct the spherical harmonic coefficients from the field
%values and compare the computed coefficients with the coefficients used to
%generate the field values.
% [C,S,a,c]=getWMMCoeffs();
% M=(1/2)*(sqrt(1+8*length(C))-1);
%
% %Evaluate the model on a grid.
% numPoints=21;
% totalGridPoints=numPoints*numPoints;
% deltaLat=180/numPoints;
% lat=-90+(0:(numPoints-1))*deltaLat;
% lat=lat*(pi/180);%Convert latitude to radians.
%
% deltaLon=360/numPoints;
% lon=-180+(0:(numPoints-1))*deltaLon;
% lon=lon*(pi/180);%Convert longitude to radians.
% [latGrid,lonGrid]=ndgrid(lat,lon);
% latLonEllipse=[latGrid(:)';lonGrid(:)'];
% %Convert from ellipsoidal latitudes to spherical latitudes (This assumes
% %that the points are on the surface of the reference ellipsoid).
% points=ellips2Sphere([latLonEllipse;zeros(1,totalGridPoints)]);
%
% %Evaluate the magnetic field values at the given points.
% V=spherHarmonicEval(C,S,points,a,c);
%
% %Reconstruct the coefficients from the magnetic field values.
% [CFit,SFit]=fitSpherHarmonics(M,points,V,a,c);
%
% %Check how good the fit is to the original set of coefficients.
% sel1=C(:)~=0;
% max(abs(CFit(sel1)-C(sel1))./C(sel1))
% sel1=S(:)~=0;
% max(abs(SFit(sel1)-S(sel1))./S(sel1))
%In both instances, the maximum relative error should have been less than
%8e-10. The fact that there is any error stems from finite precision
%effects.
%
%REFERENCES:
%[1] D. F. Crouse, "An overview of major terrestrial, celestial, and
%    temporal coordinate systems for target tracking," Naval Research
%    Laboratory, Washington, DC, Tech. Rep. NRL/FR/5344-16-10,279, 10 Aug.
%    2016.
%
%June 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<5||isempty(c))
    c=1;
end

if(nargin<4||isempty(a))
    a=1;
end

pointDims=size(points,1);

numEq=size(points,2);

%The total number of associated Legendre polynomials.
numPBar=(M+1)*(M+2)/2;

%The M+1 being subtracted is the number of S_{n,0} values,
%because these can just be set to zero as they are never used; they are
%always multiplied by sin(0).
numCCoeffs=numPBar;
numSCoeffs=numPBar-(M+1);

%The total number of coefficients if all values for C and S must be
%computed.
coeffsTotal=numCCoeffs+numSCoeffs;

if(coeffsTotal>numEq)
   warning('The value of M is too large for the number of values provided. The exact solution is unobservable.')
end

%The columns of Y corresponds to the C coefficients followed by the S
%coefficients.
Y=zeros(numEq,coeffsTotal);
for curEq=1:numEq
    curPoint=points(:,curEq);

    if(pointDims==3)
        r=curPoint(1);
        Az=curPoint(2);
        coLat=pi/2-curPoint(3);
    else
        r=1;
        Az=curPoint(1);
        coLat=pi/2-curPoint(2);
    end
    %We now have r, azimuth, and colatitude.

    %Get all of the associated Legendre polynomials.
    PBarVals=LegendreCos(coLat,M);
    PBarVector=PBarVals(:);

    outerCoeff=c/r;
    innerCoeff=1;
    for n=0:M
        for m=0:n
            %The index of PBarVals(n+1,m+1) in the vector PBarVector.
            idx=n*(n+1)/2+m+1;

            coeffVal=outerCoeff*innerCoeff*PBarVector(idx);

            Y(curEq,idx)=coeffVal*cos(m*Az);
            if(m>0)
                %The -(n+1) deals with the fact that there are no values of
                %S for m=0.
                Y(curEq,numCCoeffs+idx-(n+1))=coeffVal*sin(m*Az);
            end
        end
        innerCoeff=innerCoeff*(a/r);
    end
end

%If one were to use linsolve(Y,vals(:)), many systems would complain about
%Y having poor conditioning.
coeffs=lsqminnorm(Y,vals);

C=coeffs(1:numPBar,:);

%Fill in S, noting that the unused parts are left set to zero.
numSets=size(vals,2);
SData=zeros(numPBar,numSets);
curSIdx=numCCoeffs+1;
for n=1:M
    for m=1:n
        idx=n*(n+1)/2+m+1;
        SData(idx,:)=coeffs(curSIdx,:);
        curSIdx=curSIdx+1;
    end
end
S=SData;

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 net = pairwise(arg)

% PAIRWISE
%
% Construct a pairwise multi-class support vector classification network.
%
% Examples:
%
%    % default constructor (a 0-class pairwise network!)
%
%    net1 = pairwise;
%
%    % copy constructor
%
%    net2 = pairwise(net1);
%
%    % construct pairwise multi-class svc from a vector of two-class networks
%
%    net3 = pairwise(net)

%
% File        : @pairwise/pairwise.m
%
% Date        : Wednesday 13th September 2000
%
% Author      : Dr Gavin C. Cawley
%
% Description : Constructor for a class providing a framework for constructing
%               multi-class support vector classification networks from a set
%               of two-class networks using the pairwise rule.  Part of an
%               object-oriented implementation of Vapnik's Support Vector
%               Machine, as described in [1].
%
% References  : [1] V.N. Vapnik,
%                   "The Nature of Statistical Learning Theory",
%                   Springer-Verlag, New York, ISBN 0-387-94559-8,
%                   1995.
%
% History     : 13/09/2000 - v1.00
%
% Copyright   : (c) Dr Gavin C. Cawley, September 2000
%
%    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
%

if nargin == 0

   % this is the default constructor

   net.net = svc;
   net     = class(net, 'pairwise');

elseif nargin == 1

   if isa(arg, 'pairwise');

      % this is the copy constructor

      net = arg;

   end

elseif nargin > 1

   % there are no other constructors

   help pairwise

else

   net.net = arg;
   net     = class(net, 'pairwise');

end

% bye bye...

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%code for supplementary material
clear;
netStruct = load('../data/res152_batch32_Rankloss_2:1:0.1_margin1_both_drop0.75_shift_hard_256_152x152/net-epoch-80.mat');
net = dagnn.DagNN.loadobj(netStruct.net);
clear netStruct;
net.mode = 'test' ;
net.move('gpu') ;
net.removeLayer('RankLoss');
net.conserveMemory = true;
im_mean = reshape(net.meta.normalization.averageImage,1,1,3);
load('../url_data.mat');
p = imdb.images.data(imdb.images.set==3);
%%-------image feature
which_img = 621;
str = [ p{which_img}];
im = imread(['.',str]);
imshow(im);
oim = im;   % or oim = im;
f = getFeature2(net,oim,im_mean,'data','fc1_1bn');
f = sum(sum(f,1),2);
f2 = getFeature2(net,fliplr(oim),im_mean,'data','fc1_1bn');
f2 = sum(sum(f2,1),2);
f = f+f2;
size4 = size(f,4);
f = reshape(f,[],size4)';
f_img = norm_zzd(f);

% get raw text
load('../Flickr30k/rawdata-txt.mat');
raw_test = [];
load('../Flickr30k/train_val_test_split.mat');
test_set = find(set==3);
for i=1:1000
    tmp = test_set(i)*5-4 : test_set(i)*5;
    raw_test = cat(1,raw_test,raw_txt(tmp));
end
which_sentence = which_img*5-3;
title(raw_test{which_sentence});
load('../Flickr30k/dense_feature_word2.1.mat');
test_set = find(imdb.images.set==3);

%get text feature
content = wordcnn(:,test_set(which_img)*5-3);
len = sum(content>0);
load('/home/zzd/nlp/word2vector_matlab/flickr30k_dictionary.mat');
word_name = subset.names;
for k=1:len
    fprintf('%s ',word_name{content(k)});
end
fprintf('\n');

for i = 0:len
    content_tmp = content;
    if(i~=0) %start block words
        content_tmp(i)=0;
    end
    txtinput = zeros(len,20074,'single');
    kk = 1;
    for k=1:32
        if(content_tmp(k)==0)
            continue;
        end
        txtinput(kk,content_tmp(k))=1;
        kk = kk+1;
    end
    %transfer it to different location
    win = 33-len;
    input = zeros(32,20074,win,'single');
    for kk = 1:win
        input(kk:kk+len-1,:,kk) = txtinput;
    end

    input = reshape(input,1,32,20074,[]);
    f = getFeature2(net,input,[],'data2','fc6_2bn');
    f = sum(f,4);
    size4 = size(f,4);
    f = reshape(f,[],size4)';
    f_txt = norm_zzd(f);

    score = f_img * f_txt';
    if(i==0)
        s0 = score;
    else
      fprintf('%s,%.4f\n',word_name{content(i)},score-s0);
    end
end

author: layumi

repo: Image-Text-Embedding


text:
function [auc auh acc0 accm thrm thrs acc sens spec hull] = ...
    rocplot(scores, classes, plottype, Nthr)
%rocplot: plot a Receiver Operating Characteristic (ROC) curve
% ROC curves illustrate the performance on a binary classification problem
% where classification is based on simply thresholding a set of scores at
% varying levels. Low thresholds give high sensitivity but low specificity,
% high thresholds give high specificity but low sensitivity; the ROC curve
% plots this trade-off over a range of thresholds.
%
% Examples:
%   rocplot(scores, classes);
%   rocplot(scores, classes, plottype);
%   rocplot(scores, classes, plottype, Nthr);
%   [AUC AUH] = rocplot(scores, classes);
%   [AUC AUH acc0 accM thrM thr acc sens spec hull] = rocplot(...);
%
% classes is a Boolean vector of the same size as scores, which should
% be true where scores >= threshold should yield true, e.g. true for
% patients, and false for healthy control subjects, in medical diagnosis.
% If you have two vectors of scores, e.g. patient and control, first do:
%   scores  = [control(:); patient(:)];
%   classes = [false(numel(control), 1); true(numel(patient), 1)];
%
% plottype controls what is plotted, it defaults to 2, where:
%   0 gives no plot (useful to get AUC without creating a plot)
%   1 gives a standard ROC curve, with sensitivity vs (1 - specificity)
%   2 gives my preferred convention, with sensitivity vs specificity
% Nthr is an optional number of thresholds to consider (points in the ROC),
% if left off, all unique values in scores are considered, or, if there are
% more than 100 of those, then 100 values equally spaced from min to max of
% scores are used. Alternatively, Nthr can be a vector of thresholds.
%
% AUC is the area under the ROC curve, a measure of overall accuracy,
% which gives the probability that the classifier would rank a randomly
% chosen true instance (e.g. patient) higher than a random false one
% (e.g. control subject).
%
% AUH is the area under the convex hull of the ROC curve, which is of
% interest because it is theoretically possible to operate at any point on
% the convex hull of the points in an ROC curve (by using some proportional
% selection of classifiers that operate at two points defining the relevant
% section of the convex hull. This is just a more complex version of the
% logic that gives the null line for an ROC plot; if a threshold of -inf
% gives (spec=0,sens=1) and +inf gives (1,0), then using -inf for half of
% your data and +inf for the other half is expected to give (0.5,0.5) if
% you have equal numbers of true and false instances).
%
% Also recorded in the plot legend, and optionally returned, are acc0, the
% accuracy at a threshold of zero, which is of special importance
% in some algorithms, e.g. if your scores come from a linear classifier
% like a Support Vector Machine which can give scores(i) = w'*x(:, i) - b,
% as equivalent to testing for w'*x(:, i) > b, and accM, the max accuracy
% from the tested set of thresholds, which occurs at the threshold thrM.
%
% The function can also return a vector of all considered thresholds, along
% with the corresponding accuracies, sensitivities and specificities, in
% the variables thr, acc, sens, spec. The output hull contains the indices
% into sens and spec that give the convex hull. These arguments can then be
% used to plot multiple ROC curves and/or convex hulls on the same axis,
% e.g. after calling rocplot twice (with plottype 0), you could do:
%   plot(spec1,sens1,'b', spec2,sens2,'r', [0 1],[1 0],'g');
%
% See also: roc, prec_rec, get_concave_ROC, auroc, bookmaker
%
% These can be found in the MATLAB Central File Exchange:
% roc             http://www.mathworks.com/matlabcentral/fileexchange/19950
% roc             http://www.mathworks.com/matlabcentral/fileexchange/21318
% prec_rec        http://www.mathworks.com/matlabcentral/fileexchange/21528
% get_concave_ROC http://www.mathworks.com/matlabcentral/fileexchange/21382
% auroc           http://www.mathworks.com/matlabcentral/fileexchange/19468
% bookmaker       http://www.mathworks.com/matlabcentral/fileexchange/5648
% There is also a book "Ordinal Data Modeling", by Valen Johnson, with
% companion code  http://www.mathworks.com/matlabcentral/fileexchange/2264

% Copyright 2009 Ged Ridgway

scores = scores(:);
% need max thr greater than max(scores) so roc sens drops to zero
maxscoreplus = max(scores) + eps(max(scores));
classes = classes(:);
if any(double(logical(classes)) ~= double(classes))
    error('Classes vector must contain true (or 1) and false (or 0) only')
end

if ~exist('Nthr', 'var') || isempty(Nthr)
    if numel(scores) > 100
        Nthr = 100; % no need for more points in plot
    else
        thrs = unique([scores;maxscoreplus]);
        Nthr = numel(thrs);
    end
end
if ~exist('thrs', 'var')
    if ~isscalar(Nthr) % user-specified set of thrs to consider
        thrs = unique(Nthr); % (also sorts)
        Nthr = numel(thrs);
        % thrs(Nthr) = maxscoreplus; % commented out, assuming expert user.
    else
        thrs = linspace(min(scores), maxscoreplus, Nthr);
    end
end
thrs = thrs(:);

acc = zeros(Nthr, 1); sens = acc; spec = sens;
for n = 1:Nthr
    thr = thrs(n);
    result = scores >= thr;
    acc(n) = mean(result == classes);
    sens(n) = sum(result == 1 & classes == 1)/sum(classes == 1);
    spec(n) = sum(result == 0 & classes == 0)/sum(classes == 0);
end
% area under ROC curve
auc = trapz(spec, sens);
% test special case of 0 threshold, since training often assumes this
result0 = scores >= 0;
acc0 = mean(result0 == classes);
sens0 = sum(result0 == 1 & classes == 1)/sum(classes == 1);
spec0 = sum(result0 == 0 & classes == 0)/sum(classes == 0);
% threshold giving highest overall accuracy
[accm indm] = max(acc);
thrm = thrs(indm);
sensm = sens(indm);
specm = spec(indm);
% compute convex hull of ROC curve and area under hull
[hull auh] = convhull([spec;0], [sens;0]);
hull = unique(hull);
hull(end) = []; % don't close loop
% sanity check:
if abs(trapz(spec(hull), sens(hull)) - auh) > eps
    warning('rocplot:auh', 'Area under ROC convex hull may be incorrect');
end
if ~exist('plottype', 'var') || isempty(plottype)
    plottype = 2;
end
switch plottype
    case 0, return % no plot
    case 1 % standard TPR (= sens) vs FPR (= 1-spec)
        plot(1-spec,sens,'b', 1-spec(hull),sens(hull),'c', ...
            1-spec0,sens0,'ro', 1-specm,sensm,'ms', 1-[0 1],[1 0],'g');
        xlabel('False Positive Rate (1 - Specificity)');
        ylabel('True Positive Rate (Sensitivity)');
        legpos = 'SE';
        axis([0 1 0 1])
    case 2 % my convention, sens vs spec
        plot(spec,sens,'b', spec(hull),sens(hull),'c', ...
            spec0,sens0,'ro', specm,sensm,'ms', [0 1],[1 0],'g');
        xlabel('Specificity');
        ylabel('Sensitivity');
        legpos = 'SW';
        axis([0 1.05 0 1.05])
end
legend(sprintf('ROC Curve\nAUC = %.3g', auc), ...
    sprintf('Conv Hull\nAUC = %.3g', auh), ...
    sprintf('Acc(0) = %.3g', acc0), ...
    sprintf('Max Acc %.3g', accm), ... 'Chance', ...
    'Location', legpos)

if nargout == 0, clear auc, end % quieten rocplot(scores,classes) without ;

author: Sable

repo: mcbench-benchmarks


text:
% PRTBRV - PRT Bayesian Random Variable Class
%   Abstract Methods:
%       nDimensions

classdef prtBrv < prtAction
    methods (Abstract)
        self = estimateParameters(self, priorSelf, x)
        y = predictivePdf(self, x)
        y = predictiveLogPdf(self, x)

        d = getNumDimensions(self)

        self = initialize(self, x)
    end

    methods
        % Optional methods that would be nice to have
        % You don't have to implement these but you are welcome to some
        % types of processing for example VB with
        % calculateNegativeFreeEnergy = true, might require some of these
        function kld = conjugateKld(self, prior) %#ok<INUSD,STOUT>
            missingMethodError(self,'conjugateKld')
        end
        function x = posteriorMeanDraw(self, n, varargin) %#ok<INUSD,STOUT>
            missingMethodError(self,'posteriorMeanDraw')
        end
        function s = posteriorMeanStruct(self) %#ok<STOUT>
            missingMethodError(self,'posteriorMeanStruct')
        end
        function plotCollection(selfVec,colors) %#ok<INUSD>
           missingMethodError(selfVec(1),'plotCollection')
        end
        function val = plotLimits(self) %#ok<STOUT>
           missingMethodError(self,'plotLimits')
        end
        %function plot(self)
        %    plotCollection(self);
        %end
    end

    % All prtBrv's get the property nDimensions
    % For a particular subclass getNumDimensions must be implemented
    properties (Dependent)
        nDimensions
    end
    methods
        function val = get.nDimensions(self)
            val = zeros(size(self));
            for iRv = 1:numel(self)
                val(iRv) = self(iRv).getNumDimensions();
            end
        end
    end

    % Because we are a prtAction we also must implement something like this

    %%% properties (SetAccess=private)
    %%%     name = 'Maximum a Posteriori'   % Maximum a Posteriori
    %%%     nameAbbreviation = 'MAP'        % MAP
    %%%     isNativeMary = true;            % True
    %%% end

    % For prtAction, we also must implement trainAction() and runAction()
    methods (Access = protected, Hidden = true)
        function self = trainAction(self, ds)
            self = estimateParameters(self, ds);
        end

        function ds = runAction(self, ds)
            ds = ds.setObservations(self.predictiveLogPdf(ds.getObservations));
        end
    end

    methods (Access = 'protected', Hidden = true)
        function missingMethodError(self,methodName) %#ok<MANU>
            error('The method %s is not defined for this prtBrv object',methodName);
        end

        function self = constructorInputParse(self,varargin)

            nIn = length(varargin);

            % Quick Exit for the zero input constructor
            if nIn == 0
                return
            elseif mod(nIn,2)
                error('prt:prtBrv:constructorInputParse','Inputs must be supplied by as string/value pairs');
            end

            self = prtUtilAssignStringValuePairs(self, varargin{:});

        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
% Copyright (C) 2017-2018 Titus Cieslewski, RPG, University of Zurich,
%   Switzerland
%   You can contact the author at <titus at ifi dot uzh dot ch>
% Copyright (C) 2017-2018 Siddharth Choudhary, College of Computing,
%   Georgia Institute of Technology, Atlanta, GA, USA
% Copyright (C) 2017-2018 Davide Scaramuzza, RPG, University of Zurich,
%   Switzerland
%
% This file is part of dslam_open.
%
% dslam_open 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.
%
% dslam_open 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 dslam_open. If not, see <http://www.gnu.org/licenses/>.

function plotDataOverTime(data_increments, data_increment_times, vo_end_time)

opt_traffic = ...
    cumsum(cellfun(@(x) sum(sum(x(:, :, 1))), data_increments));
netvlad_traffic = ...
    cumsum(cellfun(@(x) sum(sum(x(:, :, 2))), data_increments));
gv_traffic = ...
    cumsum(cellfun(@(x) sum(sum(x(:, :, 3))), data_increments));

opt_filter = [true; opt_traffic(2:end) ~= opt_traffic(1:end-1)];

nv_filt = [true; netvlad_traffic(2:end) ~= netvlad_traffic(1:end-1)];
gv_filt = [true; gv_traffic(2:end) ~= gv_traffic(1:end-1)];

plot(data_increment_times(opt_filter), opt_traffic(opt_filter) / 1e6, '-x');
hold on;
plot(data_increment_times(nv_filt), netvlad_traffic(nv_filt) / 1e6);
plot(data_increment_times(gv_filt), gv_traffic(gv_filt) / 1e6);
y_lim = ylim;
plot([vo_end_time vo_end_time], ylim, 'k--');
set(gca, 'YLim', y_lim);
hold off;
legend('DOpt', 'DVPR', 'RelPose', 'VO end time', 'Location', 'NorthWest');
title('Total data transmission');
xlabel('time [s]');
ylabel('total transmitted [MB]');
grid on;

end


text:
%PARSEOPT Parse options.
%
%   [ GOT, VAL ] = PARSEOPT( COPTDEF, VARARGIN ) parses the property/value pairs
%   in VARARGIN, and adds the default options in the second column of
%   COPTDEF. The first column of COPTDEF contains strings of the field names.
%   Returns the got/val structs.

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function S = read_mclust_t(tfilelist)

% adapted from M-clust function LoadSpikes

%-------------------
% Check input type
%-------------------
if ~isa(tfilelist, 'cell')
   ft_error('LoadSpikes: tfilelist should be a cell-array.');
end

nFiles = length(tfilelist);

S = cell(nFiles, 1);
for iF = 1:nFiles
	tfn = tfilelist{iF};
	if ~isempty(tfn)
    try
 		  tfp = fopen_or_error(tfn, 'rb','b');
    catch err
			warning([ 'Could not open tfile ' tfn]);
      continue
		end

		ReadHeader(tfp);
		S{iF} = fread(tfp,inf,'uint64');
	  S{iF} = double(S{iF}*100);

		fclose(tfp);
	end 		% if tfn valid
end		% for all files
fprintf(2,'\n');

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function test_interp_nd_test03 ( m, n )

%*****************************************************************************80
%
%% TEST_INTERP_ND_TEST03 samples each derivative component in M dimensions, at N points.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 August 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of evaluation points.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST_INTERP_ND_TEST03\n' );
  fprintf ( 1, '  P00_F evaluates derivative components.\n' );
  fprintf ( 1, '  Here, we use spatial dimension M = %d\n', m );
  fprintf ( 1, '  The number of points is N = %d\n', n );

  seed = 123456789;
  [ x, seed ] = r8mat_uniform_01 ( m, n, seed );

  prob_num = p00_prob_num ( );

  for prob = 1 : prob_num

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Problem %d\n', prob );

    [ c, seed ] = p00_c ( prob, m, seed );
    r8vec_print ( m, c, '  C parameters:' );

    [ w, seed ] = p00_w ( prob, m, seed );
    r8vec_print ( m, w, '  W parameters:' );

    fprintf ( 1, '\n' );
    fprintf ( 1, '                        X(1)      X(2) ...\n' );
    fprintf ( 1, '      F(X)           dFdX(1)   dFdX(2) ...\n' );
    fprintf ( 1, '\n' );

    f = p00_f ( prob, m, c, w, n, x );

    d = zeros ( m, n );
    for id = 1 : m
      d(id,1:n) = p00_d ( prob, m, id, c, w, n, x );
    end

    for j = 1 : n
      fprintf ( 1, '                  ' );
      for i = 1 : m
        fprintf ( 1, '  %8.4f', x(i,j) );
      end
      fprintf ( 1, '\n' );

      fprintf ( 1, '  %14g  ', f(j) );
      for i = 1 : m
        fprintf ( 1, '  %8.4f', d(i,j) );
      end
      fprintf ( 1, '\n' );
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function points=regularNSimplexCoords(numDim,method)
%%REGULARNSIMPLEXCOORDS Find the coordinates of a regular n-dimensional
%               simplex. That is, a hypertetrahedron in n-dimensions.
%               All sides in the simplex have equal length, the centroid of
%               the simplex is zero, and the distance of all of the points
%               in the simplex from the origin is 1.
%
%INPUTS: numDim The number of dimensions of the simplex. numDim>=1.
%        method The simplex is not unique (it can be rotated). If the
%               method parameter is specified, it selects which method of
%               finding a simplex is used. Possible values are
%               0 (The default if omitted or an empty matrix is passed. Use
%                  the method of [1], which utilizes cosines, is used.
%               1 The method of [1] utilizing square roots instead of
%                 cosines is used.
%
%OUTPUTS: points numDimXnumPoints set of points that form a
%                numDim-simplex. numPoints is always equal to numDim+1.
%
%The two methods given are from Chapter 9 of [1]. The methods do not
%actually provide points that are a unit distance from the origin, though
%they are all equidistant from the origin. Thus, this scales the points to
%be a unit distance form the origin.
%
%REFERENCES:
%[1] A.H. Stroud, Approximate Calculation of Multiple Integrals. Cliffs,
%    NJ: Prentice-Hall, Inc., 1971.
%
%December 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<2||isempty(method))
    method=0;
end

switch(method)
    case 0%The method in [1] utilizing cosines
        points=zeros(numDim,numDim+1);

        r=1:fix(numDim/2);

        for k=0:numDim
            points(2*r-1,k+1)=cos(2*r*k*pi/(numDim+1));
            points(2*r,k+1)=sin(2*r*k*pi/(numDim+1));
            %For an odd number of dimensions
            if(mod(numDim,2)~=0)
                points(numDim,k+1)=(-1)^k/sqrt(2);
            end
        end
        points=points/sqrt(numDim/2);
    case 1%The method in [1] utilizing square roots.
        points=zeros(numDim,numDim+1);
        for k=1:(numDim+1)
            i=1:k;
            points(i,k)=-sqrt((numDim+1)./((numDim-i+2).*(numDim-i+1)));
            points(k,k)=sqrt((numDim+1)*(numDim-k+1)/(numDim-k+2));
            %Components for i>k are zero.
        end
        points=points/sqrt(numDim);
    otherwise
       error('Invalid method specified')
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 [X_resid, Y_resid, handles] = prplot_multilevel(Y, X, wh_col)
% Partial correlation plot for multi-level analysis
%
% :Usage:
% ::
%
%     [X_resid, Y_resid, handles] = prplot_multilevel(Y, X, wh_col)
%
% Uses unweighted estimates.
%
% do not enter intercept in X

N = length(X);

for i = 1:N

    X_other{i} = X{i};
    X_other{i}(:, wh_col) = [];

    X_part{i} = X{i}(:, wh_col);

end

wh = true(size(X{1}, 2) + 1, 1); % assume X is same size; k vars + intercept

wh(wh_col + 1) = false;  % omit from betas; consider intercept

stats = glmfit_multilevel(Y, X_other, [], 'noverbose');

for i = 1:N

    Y_resid{i} = Y{i} - [ones(size(X_other{i}, 1), 1) X_other{i}] * stats.first_level.beta(:, i);

end

stats = glmfit_multilevel(X_part, X_other, [], 'noverbose');

for i = 1:N

    X_resid{i} = X_part{i} - [ones(size(X_other{i}, 1), 1) X_other{i}] * stats.first_level.beta(:, i);

end

handles.overall = plot(cat(1, X_resid{:}), cat(1, Y_resid{:}), 'ko');

% partial fit plot
for i = 1:N

    handles.indiv = plot(X_resid{i}, Y_resid{i},'Color', rand(1, 3));

end

handles.refline = refline;
set(handles.refline, 'LineWidth', 3);

end % function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function im = UpsampleImage(img, line_factor, elem_factor)
% Returns an FFT upsampled image.
%
% Returns an image which is upsampled by the inverse FFTing, zero padding
% by the supplied amount (line_factor and elem_factor) and FFTing back.
%
% Usage:
%     im = UpsampleImage(img, line_factor, elem_factor)
%
% Example:
%     im = UpsampleImage(some_image data, 1.5, 1.5);
%
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

    [rows,cols] = size(img);

    % Inverse FFT
    sp = fftshift(ifft2(fftshift(img))) / (rows*cols);
    %%%sp = ifft2(img) / (rows*cols);

    padded_data = zeros( floor(rows*line_factor), floor(cols*elem_factor) );
    if (line_factor < 1  &&  elem_factor < 1)
        row_start  = floor((rows-(rows*line_factor))/2);
        row_length = size(padded_data,1);
        col_start  = floor((cols-(cols*elem_factor))/2);
        col_length = size(padded_data,2);
        padded_data = sp(row_start:(row_start+row_length), ...
                         col_start:(col_start+col_length));
    else
        % Zero pad around the edges
        off_line = floor((rows*line_factor - rows)/2.0);
        off_elem = floor((cols*elem_factor - cols)/2.0) ;
        padded_data(off_line:off_line+rows-1, off_elem:off_elem+cols-1) = sp;
    end

    % Then just FFT back to get the image.
    im = fftshift(fft2(fftshift(padded_data))) * (rows*cols);
    %%%im = fft2(padded_data) / (rows*cols);
end
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
% DPRIME  --  Signal-detection theory sensitivity measure.
%
%  d = dprime(pHit,pFA)
%  [d,beta] = dprime(pHit,pFA)
%
%  PHIT and PFA are numerical arrays of the same shape.
%  PHIT is the proportion of "Hits":        P(Yes|Signal)
%  PFA is the proportion of "False Alarms": P(Yes|Noise)
%  All numbers involved must be between 0 and 1.
%  The function calculates the d-prime measure for each <H,FA> pair.
%  The criterion value BETA can also be requested.
%  Requires MATLAB's Statistical Toolbox.
%
%  References:
%  * Green, D. M. & Swets, J. A. (1974). Signal Detection Theory and
%    Psychophysics (2nd Ed.). Huntington, NY: Robert Krieger Publ.Co.
%  * Macmillan, Neil A. & Creelman, C. Douglas (2005). Detection Theory:
%    A User's Guide (2nd Ed.). Lawrence Erlbaum Associates.
%
%  See also NORMINV, NORMPDF.

% Original coding by Alexander Petrov, Ohio State University.
% $Revision: 1.2 $  $Date: 2009-02-09 10:49:29 $
%
% Part of the utils toolbox version 1.1 for MATLAB version 5 and up.
% http://alexpetrov.com/softw/utils/
% Copyright (c) Alexander Petrov 1999-2006, http://alexpetrov.com
% Please read the LICENSE and NO WARRANTY statement:

% GNU Public License for the UTILS Toolbox
%
% ==============================================================================
%
% IN BRIEF
% ========
%
% This document refers to all Matlab scripts and documentation (referred
% to collectively here as "the software") contained in the "utils toolbox"
% by Alexander Petrov 1999-2006, http://alexpetrov.com/softw/utils/
%
% The software is freely available and freely redistributable, according
% to the conditions of the Gnu General Public License (below). You may not
% distribute the software, in whole or in part, in conjunction with
% proprietary code. That means you ONLY have my permission to distribute a
% program that uses my code IF you also make freely available (under the
% terms of the Gnu GPL) the source code for your whole project. You may
% not pass on the software to another party in its current form or any
% altered, embellished or reduced form, without acknowledging the author
% and including a copy of this license.
%
% The software is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
% Public License (reproduced below), and the Free Software Foundation
% website (http://www.fsf.org) for more details.
%
% Please notify the author, via the website, of any bugs, notes, comments
% or suggested changes, particularly of any useful changes you may have
% made to your own copy of the software.
%
% Alex Petrov, December 2006
%
% ==============================================================================
%
% GNU GENERAL PUBLIC LICENSE
% ==========================
%
% Version 2, June 1991
% Copyright (C) 1989, 1991 Free Software Foundation, Inc.
% 	59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
% Everyone is permitted to copy and distribute verbatim copies
% of this license document, but changing it is not allowed.
%
%
% Preamble
% ========
%
%   The licenses for most software are designed to take away your
% freedom to share and change it.  By contrast, the GNU General Public
% License is intended to guarantee your freedom to share and change free
% software--to make sure the software is free for all its users.  This
% General Public License applies to most of the Free Software
% Foundation's software and to any other program whose authors commit to
% using it.  (Some other Free Software Foundation software is covered by
% the GNU Library General Public License instead.)  You can apply it to
% your programs, too.
%
%   When we speak of free software, we are referring to freedom, not
% price.  Our General Public Licenses are designed to make sure that you
% have the freedom to distribute copies of free software (and charge for
% this service if you wish), that you receive source code or can get it
% if you want it, that you can change the software or use pieces of it
% in new free programs; and that you know you can do these things.
%
%   To protect your rights, we need to make restrictions that forbid
% anyone to deny you these rights or to ask you to surrender the rights.
% These restrictions translate to certain responsibilities for you if you
% distribute copies of the software, or if you modify it.
%
%   For example, if you distribute copies of such a program, whether
% gratis or for a fee, you must give the recipients all the rights that
% you have.  You must make sure that they, too, receive or can get the
% source code.  And you must show them these terms so they know their
% rights.
%
%   We protect your rights with two steps: (1) copyright the software, and
% (2) offer you this license which gives you legal permission to copy,
% distribute and/or modify the software.
%
%   Also, for each author's protection and ours, we want to make certain
% that everyone understands that there is no warranty for this free
% software.  If the software is modified by someone else and passed on, we
% want its recipients to know that what they have is not the original, so
% that any problems introduced by others will not reflect on the original
% authors' reputations.
%
%   Finally, any free program is threatened constantly by software
% patents.  We wish to avoid the danger that redistributors of a free
% program will individually obtain patent licenses, in effect making the
% program proprietary.  To prevent this, we have made it clear that any
% patent must be licensed for everyone's free use or not licensed at all.
%
%   The precise terms and conditions for copying, distribution and
% modification follow.
%
%
% GNU GENERAL PUBLIC LICENSE
% ==========================
% TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION
% ===============================================================
%
%   0. This License applies to any program or other work which contains
% a notice placed by the copyright holder saying it may be distributed
% under the terms of this General Public License.  The "Program", below,
% refers to any such program or work, and a "work based on the Program"
% means either the Program or any derivative work under copyright law:
% that is to say, a work containing the Program or a portion of it,
% either verbatim or with modifications and/or translated into another
% language.  (Hereinafter, translation is included without limitation in
% the term "modification".)  Each licensee is addressed as "you".
%
% Activities other than copying, distribution and modification are not
% covered by this License; they are outside its scope.  The act of
% running the Program is not restricted, and the output from the Program
% is covered only if its contents constitute a work based on the
% Program (independent of having been made by running the Program).
% Whether that is true depends on what the Program does.
%
%   1. You may copy and distribute verbatim copies of the Program's
% source code as you receive it, in any medium, provided that you
% conspicuously and appropriately publish on each copy an appropriate
% copyright notice and disclaimer of warranty; keep intact all the
% notices that refer to this License and to the absence of any warranty;
% and give any other recipients of the Program a copy of this License
% along with the Program.
%
% You may charge a fee for the physical act of transferring a copy, and
% you may at your option offer warranty protection in exchange for a fee.
%
%   2. You may modify your copy or copies of the Program or any portion
% of it, thus forming a work based on the Program, and copy and
% distribute such modifications or work under the terms of Section 1
% above, provided that you also meet all of these conditions:
%
%     a) You must cause the modified files to carry prominent notices
%     stating that you changed the files and the date of any change.
%
%     b) You must cause any work that you distribute or publish, that in
%     whole or in part contains or is derived from the Program or any
%     part thereof, to be licensed as a whole at no charge to all third
%     parties under the terms of this License.
%
%     c) If the modified program normally reads commands interactively
%     when run, you must cause it, when started running for such
%     interactive use in the most ordinary way, to print or display an
%     announcement including an appropriate copyright notice and a
%     notice that there is no warranty (or else, saying that you provide
%     a warranty) and that users may redistribute the program under
%     these conditions, and telling the user how to view a copy of this
%     License.  (Exception: if the Program itself is interactive but
%     does not normally print such an announcement, your work based on
%     the Program is not required to print an announcement.)
%
% These requirements apply to the modified work as a whole.  If
% identifiable sections of that work are not derived from the Program,
% and can be reasonably considered independent and separate works in
% themselves, then this License, and its terms, do not apply to those
% sections when you distribute them as separate works.  But when you
% distribute the same sections as part of a whole which is a work based
% on the Program, the distribution of the whole must be on the terms of
% this License, whose permissions for other licensees extend to the
% entire whole, and thus to each and every part regardless of who wrote it.
%
% Thus, it is not the intent of this section to claim rights or contest
% your rights to work written entirely by you; rather, the intent is to
% exercise the right to control the distribution of derivative or
% collective works based on the Program.
%
% In addition, mere aggregation of another work not based on the Program
% with the Program (or with a work based on the Program) on a volume of
% a storage or distribution medium does not bring the other work under
% the scope of this License.
%
%   3. You may copy and distribute the Program (or a work based on it,
% under Section 2) in object code or executable form under the terms of
% Sections 1 and 2 above provided that you also do one of the following:
%
%     a) Accompany it with the complete corresponding machine-readable
%     source code, which must be distributed under the terms of Sections
%     1 and 2 above on a medium customarily used for software interchange; or,
%
%     b) Accompany it with a written offer, valid for at least three
%     years, to give any third party, for a charge no more than your
%     cost of physically performing source distribution, a complete
%     machine-readable copy of the corresponding source code, to be
%     distributed under the terms of Sections 1 and 2 above on a medium
%     customarily used for software interchange; or,
%
%     c) Accompany it with the information you received as to the offer
%     to distribute corresponding source code.  (This alternative is
%     allowed only for noncommercial distribution and only if you
%     received the program in object code or executable form with such
%     an offer, in accord with Subsection b above.)
%
% The source code for a work means the preferred form of the work for
% making modifications to it.  For an executable work, complete source
% code means all the source code for all modules it contains, plus any
% associated interface definition files, plus the scripts used to
% control compilation and installation of the executable.  However, as a
% special exception, the source code distributed need not include
% anything that is normally distributed (in either source or binary
% form) with the major components (compiler, kernel, and so on) of the
% operating system on which the executable runs, unless that component
% itself accompanies the executable.
%
% If distribution of executable or object code is made by offering
% access to copy from a designated place, then offering equivalent
% access to copy the source code from the same place counts as
% distribution of the source code, even though third parties are not
% compelled to copy the source along with the object code.
%
%   4. You may not copy, modify, sublicense, or distribute the Program
% except as expressly provided under this License.  Any attempt
% otherwise to copy, modify, sublicense or distribute the Program is
% void, and will automatically terminate your rights under this License.
% However, parties who have received copies, or rights, from you under
% this License will not have their licenses terminated so long as such
% parties remain in full compliance.
%
%   5. You are not required to accept this License, since you have not
% signed it.  However, nothing else grants you permission to modify or
% distribute the Program or its derivative works.  These actions are
% prohibited by law if you do not accept this License.  Therefore, by
% modifying or distributing the Program (or any work based on the
% Program), you indicate your acceptance of this License to do so, and
% all its terms and conditions for copying, distributing or modifying
% the Program or works based on it.
%
%   6. Each time you redistribute the Program (or any work based on the
% Program), the recipient automatically receives a license from the
% original licensor to copy, distribute or modify the Program subject to
% these terms and conditions.  You may not impose any further
% restrictions on the recipients' exercise of the rights granted herein.
% You are not responsible for enforcing compliance by third parties to
% this License.
%
%   7. If, as a consequence of a court judgment or allegation of patent
% infringement or for any other reason (not limited to patent issues),
% conditions are imposed on you (whether by court order, agreement or
% otherwise) that contradict the conditions of this License, they do not
% excuse you from the conditions of this License.  If you cannot
% distribute so as to satisfy simultaneously your obligations under this
% License and any other pertinent obligations, then as a consequence you
% may not distribute the Program at all.  For example, if a patent
% license would not permit royalty-free redistribution of the Program by
% all those who receive copies directly or indirectly through you, then
% the only way you could satisfy both it and this License would be to
% refrain entirely from distribution of the Program.
%
% If any portion of this section is held invalid or unenforceable under
% any particular circumstance, the balance of the section is intended to
% apply and the section as a whole is intended to apply in other
% circumstances.
%
% It is not the purpose of this section to induce you to infringe any
% patents or other property right claims or to contest validity of any
% such claims; this section has the sole purpose of protecting the
% integrity of the free software distribution system, which is
% implemented by public license practices.  Many people have made
% generous contributions to the wide range of software distributed
% through that system in reliance on consistent application of that
% system; it is up to the author/donor to decide if he or she is willing
% to distribute software through any other system and a licensee cannot
% impose that choice.
%
% This section is intended to make thoroughly clear what is believed to
% be a consequence of the rest of this License.
%
%   8. If the distribution and/or use of the Program is restricted in
% certain countries either by patents or by copyrighted interfaces, the
% original copyright holder who places the Program under this License
% may add an explicit geographical distribution limitation excluding
% those countries, so that distribution is permitted only in or among
% countries not thus excluded.  In such case, this License incorporates
% the limitation as if written in the body of this License.
%
%   9. The Free Software Foundation may publish revised and/or new versions
% of the General Public License from time to time.  Such new versions will
% be similar in spirit to the present version, but may differ in detail to
% address new problems or concerns.
%
% Each version is given a distinguishing version number.  If the Program
% specifies a version number of this License which applies to it and "any
% later version", you have the option of following the terms and conditions
% either of that version or of any later version published by the Free
% Software Foundation.  If the Program does not specify a version number of
% this License, you may choose any version ever published by the Free Software
% Foundation.
%
%   10. If you wish to incorporate parts of the Program into other free
% programs whose distribution conditions are different, write to the author
% to ask for permission.  For software which is copyrighted by the Free
% Software Foundation, write to the Free Software Foundation; we sometimes
% make exceptions for this.  Our decision will be guided by the two goals
% of preserving the free status of all derivatives of our free software and
% of promoting the sharing and reuse of software generally.
%
% NO WARRANTY:
%
%   11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY
% FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW.  EXCEPT WHEN
% OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES
% PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED
% OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.  THE ENTIRE RISK AS
% TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU.  SHOULD THE
% PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING,
% REPAIR OR CORRECTION.
%
%   12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
% WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR
% REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES,
% INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING
% OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED
% TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY
% YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER
% PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGES.
%
% END OF TERMS AND CONDITIONS
%
%
% How to Apply These Terms to Your New Programs
% =============================================
%
%   If you develop a new program, and you want it to be of the greatest
% possible use to the public, the best way to achieve this is to make it
% free software which everyone can redistribute and change under these terms.
%
%   To do so, attach the following notices to the program.  It is safest
% to attach them to the start of each source file to most effectively
% convey the exclusion of warranty; and each file should have at least
% the "copyright" line and a pointer to where the full notice is found.
%
%     <one line to give the program's name and a brief idea of what it does.>
%     Copyright (C) <year>  <name of author>
%
%     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
%
%
% Also add information on how to contact you by electronic and paper mail.
%
% If the program is interactive, make it output a short notice like this
% when it starts in an interactive mode:
%
%     Gnomovision version 69, Copyright (C) year name of author
%     Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
%     This is free software, and you are welcome to redistribute it
%     under certain conditions; type `show c' for details.
%
% The hypothetical commands `show w' and `show c' should show the appropriate
% parts of the General Public License.  Of course, the commands you use may
% be called something other than `show w' and `show c'; they could even be
% mouse-clicks or menu items--whatever suits your program.
%
% You should also get your employer (if you work as a programmer) or your
% school, if any, to sign a "copyright disclaimer" for the program, if
% necessary.  Here is a sample; alter the names:
%
%   Yoyodyne, Inc., hereby disclaims all copyright interest in the program
%   `Gnomovision' (which makes passes at compilers) written by James Hacker.
%
%   <signature of Ty Coon>, 1 April 1989
%   Ty Coon, President of Vice
%
% This General Public License does not permit incorporating your program into
% proprietary programs.  If your program is a subroutine library, you may
% consider it more useful to permit linking proprietary applications with the
% library.  If this is what you want to do, use the GNU Library General
% Public License instead of this License.
%
% ==============================================================================

function [d,beta] = dprime(pHit,pFA)

%-- Convert to Z scores, no error checking
zHit = norminv(pHit) ;
zFA  = norminv(pFA) ;

%-- Calculate d-prime
d = zHit - zFA ;

%-- If requested, calculate BETA
if (nargout > 1)
  yHit = normpdf(zHit) ;
  yFA  = normpdf(zFA) ;
  beta = yHit ./ yFA ;
end

%%  Return DPRIME and possibly BETA
%%%%%% End of file DPRIME.M

author: PatternRecognition

repo: OpenBMI


text:
function ciu = consensus_und(d,tau,reps)
% CONSENSUS_UND      Consensus clustering
%
%   CIU = CONSENSUS(D,TAU,REPS) seeks a consensus partition of the
%   agreement matrix D. The algorithm used here is almost identical to the
%   one introduced in Lancichinetti & Fortunato (2012): The agreement
%   matrix D is thresholded at a level TAU to remove an weak elements. The
%   resulting matrix is then partitions REPS number of times using the
%   Louvain algorithm (in principle, any clustering algorithm that can
%   handle weighted matrixes is a suitable alternative to the Louvain
%   algorithm and can be substituted in its place). This clustering
%   produces a set of partitions from which a new agreement is built. If
%   the partitions have not converged to a single representative partition,
%   the above process repeats itself, starting with the newly built
%   agreement matrix.
%
%   NOTE: In this implementation, the elements of the agreement matrix must
%   be converted into probabilities.
%
%   NOTE: This implementation is slightly different from the original
%   algorithm proposed by Lanchichinetti & Fortunato. In its original
%   version, if the thresholding produces singleton communities, those
%   nodes are reconnected to the network. Here, we leave any singleton
%   communities disconnected.
%
%   Inputs:     D,      agreement matrix with entries between 0 and 1
%                       denoting the probability of finding node i in the
%                       same cluster as node j
%               TAU,    threshold which controls the resolution of the
%                       reclustering
%               REPS,   number of times that the clustering algorithm is
%                       reapplied
%
%   Outputs:    CIU,    consensus partition
%
%   References: Lancichinetti & Fortunato (2012). Consensus clustering in
%   complex networks. Scientific Reports 2, Article number: 336.
%
%   Richard Betzel, Indiana University, 2012
%
%   modified on 3/2014 to include "unique_partitions"

n = length(d); flg = 1;
while flg == 1

    flg = 0;
    dt = d.*(d >= tau).*~eye(n);
    if nnz(dt) == 0
        ciu = (1:n)';
    else
        ci = zeros(n,reps);
        for iter = 1:reps
            ci(:,iter) = community_louvain(dt);
        end
        ci = relabel_partitions(ci);
        ciu = unique_partitions(ci);
        nu = size(ciu,2);
        if nu > 1
            flg = 1;
            d = agreement(ci)./reps;
        end
    end

end

function cinew = relabel_partitions(ci)
[n,m] = size(ci);
cinew = zeros(n,m);
for i = 1:m
    c = ci(:,i);
    d = zeros(size(c));
    count = 0;
    while sum(d ~= 0) < n
        count = count + 1;
        ind = find(c,1,'first');
        tgt = c(ind);
        rep = c == tgt;
        d(rep) = count;
        c(rep) = 0;
    end
    cinew(:,i) = d;
end

function ciu = unique_partitions(ci)
ci = relabel_partitions(ci);
ciu = [];
count = 0;
c = 1:size(ci,2);
while ~isempty(ci)
    count = count + 1;
    tgt = ci(:,1);
    ciu = [ciu,tgt];                %#ok<AGROW>
    dff = sum(abs(bsxfun(@minus,ci,tgt))) == 0;
    ci(:,dff) = [];
    c(dff) = [];
end


text:
%% Examples of Numerical Solutions of PDEs
% We shall use a sequence of examples to introduce basic components of
% numerical solutions of PDEs.

%% Smooth solutions on regular domains
% We solve the Poisson equation $-\Delta u =1$ in the unit square $\Omega =
% (0,1)\times (0,1)$ with homogenous Dirichlet boundary condition
% $u|_{\partial \Omega}=0$ using linear finite element discrization on a
% sequence of uniform and structured grids.

close all; clear all
figure(1); set(gcf,'Units','normal'); set(gcf,'Position',[0,0,0.45,0.3]);
node = [0,0; 1,0; 1,1; 0,1];    % nodes
elem = [2,3,1; 4,1,3];          % elements
subplot(1,3,1); showmesh(node,elem)
[node,elem] = uniformrefine(node,elem);
subplot(1,3,2); showmesh(node,elem)
[node,elem] = uniformrefine(node,elem);
subplot(1,3,3); showmesh(node,elem)

%%
squarePoisson
% See squarePoisson for details of computation.

%% Smooth solutions on irregular domains
% We solve the Poisson equation in a Lake-type domain. This example shows
% the flexibility of FEM to complex domains.

load lakemesh
f = inline('ones(size(x,1),1)');    % right hand side
g_D = inline('zeros(size(x,1),1)'); % Dirichlet boundary condition
u = Poisson(node,elem,[],f,g_D,[]); % Poisson equation
% graphic
figure(1); set(gcf,'Units','normal'); set(gcf,'Position',[0,0,0.45,0.3]);
subplot(1,2,1); showmesh(node,elem); pause(0.05)
subplot(1,2,2); showsolution(node,elem,u,[0,90]);axis equal; colorbar;

%% Tutorial of iFEM

%% Features of iFEM
%
% * Simple data structures
% * Mesh adaptation in two- and three-dimensions
% * Fast solvers of algebraic equations
% * Easy to use, Easy to code, Easy to debug
% * Efficient programming

%% Basic Data Stucture
%
% *Mesh: node and elem.*
%
% Two matrices |node(1:N,1:d)| and |elem(1:NT,1:d+1)| are used to represent
% a d-dimensional triangulation embedded in $R^d$, where |N| is the number
% of vertices and |NT| is the number of elements.
%
% |node(k,1)| and |node(k,2)| are the x- and y-coordinates of the k-th node
% for points in 2-D. In 3-D, |node(k,3)| gives the additional z-coordinates
% of the k-th node.
%
% |elem(t,1:d+1)| are the global indices of |d+1| vertices which form the
% abstract $d$-simplex |t|. By convention, the vertices of a simplex is ordered
% such that the signed volume is positive. Therefore in 2-D, three vertices
% of a triangle is ordered counterclockwise and in 3-D, the ordering of
% vertices follows the right-hand rule.
%
% Related functions: <matlab:doc('fixorientation') fixorientation>,
% <matlab:doc('label') label>, <matlab:doc('label3') label3>
%
% Documentation: <matlab:ifemdoc('meshdoc') meshdoc>

clear all; close all;

%%
% *Example: L-shape domain in 2-D.*

node = [1,0; 1,1; 0,1; -1,1; -1,0; -1,-1; 0,-1; 0,0];
elem = [1,2,8; 3,8,2; 8,3,5; 4,5,3; 7,8,6; 5,6,8];
figure(1)
showmesh(node,elem)
axis on
findnode(node)
findelem(node,elem)

%%
% *Example: Cube in 3-D.*

node = [-1,-1,-1; 1,-1,-1; 1,1,-1; -1,1,-1; -1,-1,1; 1,-1,1; 1,1,1; -1,1,1];
elem = [1,2,3,7; 1,6,2,7; 1,5,6,7; 1,8,5,7; 1,4,8,7; 1,3,4,7];
figure(2)
showmesh3(node,elem,[],'FaceAlpha',0.25);
view([-53,8]);
axis on
findnode3(node)
findelem3(node,elem)

%%
% *Boundary: bdEdge or bdFace.*
%
% For 2-D triangulations, we use |bdEdge(1:NT,1:3)| to record the type of
% three edges of each element. Similarly in 3-D, we use |bdFace(1:NT,1:4)|
% to record the type of four faces of each element.
% The value is the type of boundary condition listed as follows.
%
% * 0: non-boundary, i.e., an interior edge or face.
% * 1: first type, i.e., a Dirichlet boundary edge or face.
% * 2: second type, i.e., a Neumann boundary edge or face.
% * 3: third type, i.e., a Robin boundary edge or face.
%
% We label three edges of a triangle such that |bdEdge(t,i)| is the edge
% opposite to the i-th vertex. Similarly |bdFace(t,i)| is the face opposite
% to the i-th vertex.
%
% Related functions: <matlab:doc('findboundary') findboundary>,
% <matlab:doc('setboundary') setboundary>, <matlab:doc('findboundary3') findboundary3>
% <matlab:doc('setboundary3') setboundary3>.
%
% Documentation: <matlab:ifemdoc('bddoc') bddoc>

node = [1,0; 1,1; 0,0];
elem = [1 2 3];
edge = [2 3; 1 3; 1 2];
subplot(1,2,1);
showmesh(node,elem);
findedge(node,edge);
findnode(node);
node = [1,1,1; 0,0,0; 1,1,0; 1,0,0];
elem = [1,2,3,4];
subplot(1,2,2);
showmesh3(node,elem,[],'FaceAlpha',0.35); view([-10,18]);
findnode3(node);
findelem3(node,[2,3,4; 1,3,4; 1,2,4; 1,2,3])

%% Finite Element Method
%
%
%%
% * Example 1: The Poisson equation in 2-D
% * Example 2: The Poisson equation: complex domains
% * Example 2: The Poisson equation in 3-D with multigrid solvers

%% Adaptive Finite Element Method
%
%
%%
% * Example 1: The Poisson equation on a 2-D L-shaped domain
% * Example 2: The Poisson equation in 3-D with jump diffusion coefficients

%% Time-dependent Problems
%%
% * Example 1: Heat equation
% * Example 2: Moving interface

repo: Aero-matlab


text:
% SIMPDEMO
% This program solves the simple two dimensional problem in Chapter 2,
% and makes figure 2.1.
%
tol=[1.d-6,1.d-6];
%
% Create the mesh for the contour plot of \| f \|.
%
vl=.1:.5:2; vr=2:4:40; v=[vl,vr];
v=.5:4:40;
v=[.25,.5:2:40];
xr=-5:.2:5; n=length(xr); z=zeros(n,n);
for i=1:n
for j=1:n
    w=[xr(i),xr(j)]';
    z(i,j)=norm(simple(w));
end
end
%
% Newton's method
%
params=[40, 1, 0,0];
%
% x0 is a good initial iterate.
%
x0=[2,.5]';
[sn, errsn, ierrn, x_hist]=nsold(x0, 'simple', tol, params);
%
% x1 is a poor initial iterate. The iteration from x1 will stagnate
%    at a point where $F'$ is singular.
%
x1=[3,5]';
[sn2, errsn2, ierrn2, x_hist2]=nsold(x1, 'simple', tol, params);
%
% Draw a contour plot of $\| f \|$.
%
figure(1)
contour(xr,xr,z,v)
hold
%
% Use the x_hist array to plot the iterations on the coutour plot.
%
plot(x_hist(1,:),x_hist(2,:),'-*', x_hist2(1,:),x_hist2(2,:),'-o');
legend('Convergence','Stagnation');
xlabel('x_1');
ylabel('x_2');
axis([0 5 -5 5])
hold
%
% Default
%
params=[40, 1000, .5, 0];
[sc, errsd, ierrd]=nsold(x0, 'simple', tol, params);
inewt=length(errsn(:,1)); cnewt=0:inewt-1;
idefault=length(errsd(:,1)); cdefault=0:idefault-1;
figure(2)
semilogy(cnewt,errsn(:,1),'-',cdefault,errsd(:,1),'--')
legend('Newton','default')


text:
function i4col_sorted_singleton_count_test ( )

%*****************************************************************************80
%
%% I4COL_SORTED_SINGLETON_COUNT_TEST tests I4COL_SORTED_SINGLETON_COUNT;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 November 2014
%
%  Author:
%
%    John Burkardt
%
  m = 3;
  n = 10;
  test_num = 2;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4COL_SORTED_SINGLETON_COUNT_TEST\n' );
  fprintf ( 1, '  I4COL_SORTED_SINGLETON_COUNT counts singletons\n' );
  fprintf ( 1, '  in a sorted I4COL;\n' );

  seed = 123456789;

  for test = 1 : test_num

    b = 0;
    c = 3;

    [ a, seed ] = i4mat_uniform_ab ( m, n, b, c, seed );

    a = i4col_sort_a ( m, n, a );

    i4mat_print ( m, n, a, '  Ascending sorted I4COL:' );

    singleton_num = i4col_sorted_singleton_count ( m, n, a );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Number of singletons = %d\n', singleton_num );

  end

  return
end


text:
function [bbox3M, planC] = crop_for_larynx(planC,paramS,varargin)
% Custom crop function for larynx segmentation model
% AI 10/04/19

indexS = planC{end};
strC = {planC{indexS.structures}.structureName};
bboxIdx = [];

%% Return bounding box if it exists
if isfield(paramS,'saveStrToPlanCFlag') && paramS.saveStrToPlanCFlag
    if isfield(paramS,'outStrName')
        outStrName = paramS.outStrName;
    else
        outStrName = 'crop_for_larynx';
    end
    bboxIdx = getMatchingIndex(outStrName,strC,'EXACT');
end
if ~isempty(bboxIdx)
    bbox3M = getStrMask(bboxIdx,planC);
else
    %% Create union of masseters (left & right) and medial pterygoids (left & right)
    unionStrListC = paramS.structureName.chewingStructures;

    for k = 1:length(unionStrListC)
        idx = getMatchingIndex(unionStrListC{k},strC,'EXACT');
        if k==1
            scanIdx = getStructureAssociatedScan(idx,planC);
            sizV = size(getScanArray(scanIdx,planC));
            mask3M = false(sizV);
        end
        [strMask3M, planC] = getStrMask(idx,planC);
        mask3M = mask3M | strMask3M;
    end

    %Get limits of bbox around union str
    %mask3M = getStrMask(endStr,planC);
    if sum(mask3M(:))>1
        [minr,~,minc,maxc,mins,~] = compute_boundingbox(mask3M);
    else
        minr = 1;
        minc = 1;
        maxc = size(mask3M,2);
        mins = 1;
    end

    %% Get limits of bounding box around cropped pt outline
    outStrName = paramS.structureName.cropStructure;
    idx = getMatchingIndex(outStrName,strC,'EXACT');
    cropMask3M = getStrMask(idx,planC);
    [~,maxr,~,~,~,maxs] = compute_boundingbox(cropMask3M);

    %Get final bounding box
    bbox3M = false(size(mask3M));
    bbox3M(minr:maxr,minc:maxc,mins:maxs) = true;
end

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR