text:
function diff = getDiffuseness_SV(i_vecs)
%GETDIFFUSENESS_SV Compute the DoA spherical variance diffuseness
%
%   This routine estimates diffuseness based on the assumption that DoA
%   estimates should have a small variance around the source DoA, if the
%   direct-to-diffuse ratio is high, and be completely uniformly
%   distributed if the sound is fully diffuse. This variation can be
%   captured on the spherical variance of the DoA estimates. The method can
%   utilize intensity vectors, or any other vectors that indicate DoAs. The
%   statistics can be computed from multiple temporal observations or
%   across different frequency bands. This estimator has been used e.g. in
%
%       Politis, A., Delikaris-Manias, S. and Pulkki, V., 2015.
%       Direction-of-arrival and diffuseness estimation above spatial aliasing for symmetrical directional microphone arrays.
%       In 2015 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP).
%
%   Inputs:
%       i_vecs:  Kx3 matrix of DoA vectors, for K temporal or frequency
%           observations
%
%   Outputs:
%       diff:       diffuseness value
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% GETDIFFUSENESS_SV.M - 5/10/2016
% Archontis Politis, archontis.politis@aalto.fi
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% unit DoA vectors
doa_vecs = i_vecs ./ (sqrt(sum(i_vecs.^2,2))*ones(1,3));
% mean DoA vector
mean_doa_vec = mean(doa_vecs,1);
% spherical variance
diff = 1 - sqrt(sum(mean_doa_vec.^2, 2));

end

author: polarch

repo: Spherical-Array-Processing


text:
%% mat2strIntDouble
% Below is a demonstration of the features of the |mat2strIntDouble| function

%%
clear; close all; clc;

%% Syntax
% |[t]=mat2strIntDouble(varargin);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
clear all
close all
%Setting variables to parameters
img=imread('synthetic1.png');
if length(size(img))==3             %generalized to rgb images
	img=rgb2gray(img);
end
img= double(img);

p=0;
phi=1;  % for 1.5 octave bandwidth
K = sqrt(2*log(2))*((2^phi+1)/(2^phi-1)); % for phi denoting spatial freq bandwidth
%{ admissible phi for simple and complex cells found 2 b in [0.5 2.5] }

w0_deg=[8 16 32 64 128];
%w0_deg=[8 16 32];
theta_deg=[0 22.5 45 67.5 90 112.5 135 157.5];
%theta_deg=[0 22.5];

w0=[8 16 32 64 128].*(pi/180);
%w0=[8 16 32].*(pi/180);
theta=[0 22.5 45 67.5 90 112.5 135 157.5].*(pi/180);
%theta=[0 22.5].*(pi/180);

len_w0=length(w0);
len_theta=length(theta);

k=1;
for i=1:1:len_w0
	for j=1:1:len_theta
		grid_size=2^(i-1)
		[gabor_even, gabor_odd]=gabor_filter_new(img,K,w0(i),theta(j),p,grid_size);

		sze = size(gabor_even,3);
		%             gabor_even_store(:,:,k:1:k+sze-1)=gabor_even;
		%             gabor_odd_store(:,:,k:1:k+sze-1)=gabor_odd;

		str_title=strcat('angle:',num2str(theta_deg(j)),' freq:', num2str(w0_deg(i)));

		for n=1:1:sze
			img_filtd_even(:,:,1)=normalize(conv2(gabor_even(:,:,n),img,'same'));
			img_filtd_odd(:,:,1)= normalize(conv2(gabor_odd(:,:,n),img,'same'));
			img_filtd_even_coeff(1)=sum(sum(img_filtd_even(:,:,1).^2));
			img_filtd_odd_coeff(1)=sum(sum(img_filtd_odd(:,:,1).^2));
			img_filtd_coeff(k+n-1)=sqrt(img_filtd_even_coeff(1) + img_filtd_odd_coeff(1));
		end
		%             clear img_filtd_even;
		%             clear img_filtd_odd;
		%             clear img_filtd_even_coeff;
		%             clear img_filtd_odd_coeff;
		%             figure;
		%             subplot(1,2,1)
		%             imshow(gabor_even(:,:,n),[])
		%             title(str_title);
		%             subplot(1,2,2)
		%             imshow(img_filtd_even(:,:,k+n-1),[])
		%             title('even filtered img')
		%
		%
		%             figure;
		%             subplot(1,2,1)
		%             imshow(gabor_odd(:,:,n),[])
		%             title(str_title);
		%             subplot(1,2,2)
		%             imshow(img_filtd_odd(:,:,k+n-1),[])
		%             title('odd filtered img')
		k=k+sze;
	end
end

figure;
plot(0:1:(k-2),img_filtd_coeff,'ro')
hold on
plot(0:1:(k-2),img_filtd_coeff,'g')
xlabel('w0:theta'); ylabel('image contrast energy');grid on;
title('Image Contrast Energy computed Quadrature pair-wise')

% contrast energy of filtered image computed quadrature pair-wise for
% every orientation and freq.

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function test_bug2336

% WALLTIME 00:30:00
% MEM 4gb
% DEPENDENCY

% script to test the integration of FEM output files created by BESA MRI.
%
% Copyright (C) 2013, BESA GmbH
%
% File name: test_bug2336.m
%
% Author: Robert Spangler
% Created: 2013-11-27

% ensure that fieldtrip/extrenal/besa is on the path
ft_hastoolbox('besa', 1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% The first part of the test script pertains to the low-level reading functions

% Specify path to data folder
strPath_Data = dccnpath('/home/common/matlab/fieldtrip/data/test/bug2336/');

% Filenames of BESA MRI output files
strPath_Leadfield         = sprintf('%sBESA_MRI_Electrode_Config_34_PB_FEM_DATA.lft', strPath_Data);
strPath_SourceSpaceNodes  = sprintf('%sBESA_MRI_Electrode_Config_34_PB_FEM_DATA.loc', strPath_Data);
strPath__ElectrodeConfig  = sprintf('%sElectrode_Config_34.sfh', strPath_Data);
strPath__SkullSurface     = sprintf('%sBESA_MRI_MRI_PB_FEM_SKULL_ACPC.srf', strPath_Data);
strPath__BrainSurface     = sprintf('%sBESA_MRI_PB_FEM_BRAIN_ACPC.srf', strPath_Data);
strPath__CSFSurface       = sprintf('%sBESA_MRI_PB_FEM_CSF_ACPC.srf', strPath_Data);
strPath__ScalpSurface     = sprintf('%sBESA_MRI_PB_FEM_SCALP_ACPC.srf', strPath_Data);

% Read leadfield
[dim lf] = readBESAlft(strPath_Leadfield);
fprintf('Number of channels: %i\n', dim(1));
fprintf('Number of nodes: %i\n', dim(2));
fprintf('Number of orientations: %i\n', dim(3));

% Read source space grid nodes
[ssg IdxNeighbour] = readBESAloc(strPath_SourceSpaceNodes);
fprintf('Number of nodes: %i\n', size(ssg, 1));
fprintf('Number of neighbours/node: %i\n', size(IdxNeighbour, 1));

% Read electrode config &  fiducials
cfg_eeg = readBESAsfh(strPath__ElectrodeConfig, 'eeg');
cfg_fid = readBESAsfh(strPath__ElectrodeConfig, 'fiducial');

% Read scalp, skull, csf & brain surface
cfg_scalp = readBESAsrf(strPath__ScalpSurface);
cfg_skull = readBESAsrf(strPath__SkullSurface);
cfg_csf   = readBESAsrf(strPath__CSFSurface);
cfg_brain = readBESAsrf(strPath__BrainSurface);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% The second part of the test script pertains to the high-level FieldTrip implementation

cfg = [];
cfg.method = 'besa';
cfg.elecfile = strPath__ElectrodeConfig;
cfg.hdmfile = strPath_Leadfield;
cfg.outputfile = tempname;
cfg.smooth = 'no';
vol = ft_prepare_headmodel(cfg);

pos = ssg(1,:);
sens1 = ft_read_sens(strPath__ElectrodeConfig);
sens2 = sens1;
sens2.elecpos = sens2.elecpos+randn(size(sens2.elecpos));

[vol1, sens1] = ft_prepare_vol_sens(vol, sens1);
lf1 = ft_compute_leadfield(pos, sens1, vol1);

[vol2, sens2] = ft_prepare_vol_sens(vol, sens2);
lf2 = ft_compute_leadfield(pos, sens2, vol2);

figure
subplot(2,2,1)
ft_plot_topo3d(sens1.elecpos, lf(:,1))
subplot(2,2,2)
ft_plot_topo3d(sens1.elecpos, lf(:,2))
subplot(2,2,3)
ft_plot_topo3d(sens1.elecpos, lf(:,3))

% show 3D stuff
figure
ft_plot_dipole(pos, [0 0 1]', 'diameter', 10, 'length', 20);
ft_plot_topo3d(sens1.elecpos, lf(:,3))

[X, Y, Z] = ndgrid(1:vol.dim(1), 1:vol.dim(2), 1:vol.dim(3));
pos = ft_warp_apply(vol.transform, [X(:) Y(:) Z(:)]);
ft_plot_mesh(pos);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function test_suite = test_parametrize
%TESTPARAMETRIZE  One-line description here, please.
%
%   output = testParametrize(input)
%
%   Example
%   testParametrize
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2011-05-27,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

test_suite = functiontests(localfunctions);

function testSquare(testCase) %#ok<*DEFNU>
% Tests with a square of perimeter 40
p1 = [10 10];
p2 = [20 10];
p3 = [20 20];
p4 = [10 20];
square = [p1;p2;p3;p4];
exp = [0;10;20;30];

testCase.assertEqual(exp, parametrize(square));

function testSquare3d(testCase)
% Tests with a square of perimeter 40
p1 = [10 10 0];
p2 = [20 10 0];
p3 = [20 20 0];
p4 = [10 20 0];
square = [p1;p2;p3;p4];
exp = [0;10;20;30];

testCase.assertEqual(exp, parametrize(square));

function testSquareNormaized(testCase) %#ok<*DEFNU>
% Tests with a square of perimeter 40
p1 = [10 10];
p2 = [20 10];
p3 = [20 20];
p4 = [10 20];
square = [p1;p2;p3;p4;p1];
exp = [0;.25;.5;.75;1];

testCase.assertEqual(exp, parametrize(square, 'normalize', true));

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
%Copyright (c) October,15 2008 by Varsha Hedau, UIUC.  All rights reserved.
function [vp Pout]=ordervp(vp,h,w,Pin)
%Orders vanishing points as vertical, farther horizontal and closer
%horizontal vanishing points.
%Note this takes VP=[x1 y1; x2 y2; x3 y3 ]as input
%if P is given P is also shuffled
if nargin < 4
Pin=[];
Pout=[];
end

vptemp=vp;
dists = ((vp(:,1)-w/2).^2 + (vp(:,2)-h/2).^2).^0.5;
[vv,ii] = sort(dists,'descend');
vp = vp(ii,:);
dot1 = dot(vp(1,:)-[w/2,h/2],[1 0])/norm(vp(1,:)-[w/2,h/2]);
dot2 = dot(vp(2,:)-[w/2,h/2],[1 0])/norm(vp(2,:)-[w/2,h/2]);
if abs(dot1)>abs(dot2)
    tempvar = vp(1,:);
    vp(1,:) = vp(2,:);
    vp(2,:) = tempvar;
end

if numel(Pin)>0

ind=find(vptemp(:,1)==vp(1,1) & vptemp(:,2)==vp(1,2));
Pout=Pin(:,ind);
ind=find(vptemp(:,1)==vp(2,1) & vptemp(:,2)==vp(2,2));
Pout=[Pout Pin(:,ind)];
ind=find(vptemp(:,1)==vp(3,1) & vptemp(:,2)==vp(3,2));
Pout=[Pout Pin(:,ind)];
Pout=[Pout Pin(:,4)];% 4th is outlier
end

 return

author: zouchuhang

repo: LayoutNet


text:
function [errCode, errMsgs] = compareEEG(EEG1, EEG2, dataFlg)
tol = 1E-10;
errCode = 0;
errMsgs = {};
%% Check data
if size(EEG1.data, 1) ~= size(EEG2.data, 1)
    errCode = errCode + 1;
    errMsgs{end + 1} = 'data channel numbers do not agree';
elseif size(EEG1.data, 2) ~= size(EEG2.data, 2)
    errCode = errCode + 1;
    errMsgs{end + 1} = 'data frame numbers do not agree';
elseif dataFlg && sum(sum(abs(EEG1.data - EEG2.data))) > tol
    errCode = errCode + 1;
    errMsgs{end + 1} = 'data values do not agree';
end

%% Check times
if size(EEG1.times, 1) ~= size(EEG2.times, 1)
    errCode = errCode + 1;
    errMsgs{end + 1} = 'times sizes do not agree';
elseif size(EEG1.times, 2) ~= size(EEG2.times, 2)
    errCode = errCode + 1;
    errMsgs{end + 1} = 'times sizes not agree';
elseif sum(sum(abs(EEG1.times - EEG2.times))) > 0
    errCode = errCode + 1;
    errMsgs{end + 1} = 'times values do not agree';
end

%% Check timestamps
if (isfield(EEG1.etc, 'timestamp'))
    if size(EEG1.etc.timestamp, 1) ~= size(EEG2.etc.timestamp, 1)
        errCode = errCode + 1;
        errMsgs{end + 1} = 'timestamp sizes do not agree';
    elseif size(EEG1.etc.timestamp, 2) ~= size(EEG2.etc.timestamp, 2)
        errCode = errCode + 1;
        errMsgs{end + 1} = 'timestamp sizes not agree';
    elseif sum(sum(abs(EEG1.etc.timestamp - EEG2.etc.timestamp))) > tol
        errCode = errCode + 1;
        errMsgs{end + 1} = 'timestamp values do not agree';
    end
end

%% Check events
if length(EEG1.event) ~= length(EEG2.event)
    errCode = errCode + 1;
    errMsgs{end + 1} = 'event sizes do not agree';
else
    if sum(~strcmpi(strtrim({EEG1.event.type}), strtrim({EEG2.event.type}))) > 0
        errCode = errCode + 1;
        errMsgs{end + 1} = 'event type do not agree';
    end
    if sum(sum(abs(cell2mat({EEG1.event.latency}) - cell2mat({EEG2.event.latency})))) > 0
        errCode = errCode + 1;
        errMsgs{end + 1} = 'event latencies do not agree';
    end
    if isfield(EEG1.event, 'usertags') && ...
        sum(~strcmpi(strtrim({EEG1.event.usertags}), strtrim({EEG2.event.usertags}))) > 0
           errCode = errCode + 1;
          errMsgs{end + 1} = 'event usertags do not agree';
    end
end
%% check chanlocs
if length(EEG1.chanlocs) ~= length(EEG2.chanlocs)
    errCode = errCode + 1;
    errMsgs{end + 1} = 'chanlocs sizes do not agree';
else
    if sum(~strcmpi(strtrim({EEG1.chanlocs.type}), strtrim({EEG2.chanlocs.type}))) > 0
        errCode = errCode + 1;
        errMsgs{end + 1} = 'chanlocs type do not agree';
    end
    if sum(~strcmpi(strtrim({EEG1.chanlocs.labels}), strtrim({EEG2.chanlocs.labels}))) > 0
        errCode = errCode + 1;
        errMsgs{end + 1} = 'chanlocs labels do not agree';
    end
    if sum(sum(abs(cell2mat({EEG1.chanlocs.X}) - cell2mat({EEG1.chanlocs.X})))) > tol
        errCode = errCode + 1;
        errMsgs{end + 1} = 'chanlocs X do not agree';
    end
    if sum(sum(abs(cell2mat({EEG1.chanlocs.Y}) - cell2mat({EEG1.chanlocs.Y})))) > tol
        errCode = errCode + 1;
        errMsgs{end + 1} = 'chanlocs Y do not agree';
    end
    if sum(sum(abs(cell2mat({EEG1.chanlocs.Z}) - cell2mat({EEG1.chanlocs.Z})))) > tol
        errCode = errCode + 1;
        errMsgs{end + 1} = 'chanlocs Z do not agree';
    end
    if sum(sum(abs(cell2mat({EEG1.chanlocs.theta}) - cell2mat({EEG1.chanlocs.theta})))) > tol
        errCode = errCode + 1;
        errMsgs{end + 1} = 'chanlocs theta do not agree';
    end
    if sum(sum(abs(cell2mat({EEG1.chanlocs.radius}) - cell2mat({EEG1.chanlocs.radius})))) > tol
        errCode = errCode + 1;
        errMsgs{end + 1} = 'chanlocs radius do not agree';
    end
    if sum(sum(abs(cell2mat({EEG1.chanlocs.sph_theta}) - cell2mat({EEG1.chanlocs.sph_theta})))) > tol
        errCode = errCode + 1;
        errMsgs{end + 1} = 'chanlocs sph_theta do not agree';
    end
    if sum(sum(abs(cell2mat({EEG1.chanlocs.sph_phi}) - cell2mat({EEG1.chanlocs.sph_phi})))) > 0
        errCode = errCode + 1;
        errMsgs{end + 1} = 'chanlocs sph_phi do not agree';
    end
end

author: VisLab

repo: EEG-Clean-Tools


text:
function [ x, seed ] = r8vec_normal_01 ( n, seed )

%*****************************************************************************80
%
%% R8VEC_NORMAL_01 returns a unit pseudonormal R8VEC.
%
%  Discussion:
%
%    The standard normal probability distribution function (PDF) has
%    mean 0 and standard deviation 1.
%
%    This routine can generate a vector of values on one call.  It
%    has the feature that it should provide the same results
%    in the same order no matter how we break up the task.
%
%    Before calling this routine, the user may call RANDOM_SEED
%    in order to set the seed of the random number generator.
%
%    The Box-Muller method is used, which is efficient, but
%    generates an even number of values each time.  On any call
%    to this routine, an even number of new values are generated.
%    Depending on the situation, one value may be left over.
%    In that case, it is saved for the next call.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 July 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of values desired.  If N is negative,
%    then the code will flush its internal memory; in particular,
%    if there is a saved value to be used on the next call, it is
%    instead discarded.  This is useful if the user has reset the
%    random number seed, for instance.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real X(N), a sample of the standard normal PDF.
%
%    Output, integer SEED, an updated seed for the random number generator.
%
%  Local parameters:
%
%    Local, integer MADE, records the number of values that have
%    been computed.  On input with negative N, this value overwrites
%    the return value of N, so the user can get an accounting of
%    how much work has been done.
%
%    Local, real R(N+1), is used to store some uniform random values.
%    Its dimension is N+1, but really it is only needed to be the
%    smallest even number greater than or equal to N.
%
%    Local, integer SAVED, is 0 or 1 depending on whether there is a
%    single saved value left over from the previous call.
%
%    Local, integer X_LO_INDEX, X_HI_INDEX, records the range of entries of
%    X that we need to compute.  This starts off as 1:N, but is adjusted
%    if we have a saved value that can be immediately stored in X(1),
%    and so on.
%
%    Local, real Y, the value saved from the previous call, if
%    SAVED is 1.
%
  persistent made;
  persistent saved;
  persistent y;
%
%  I'd like to allow the user to reset the internal data.
%  But this won't work properly if we have a saved value Y.
%  I'm making a crock option that allows the user to signal
%  explicitly that any internal memory should be flushed,
%  by passing in a negative value for N.
%
  if ( n < 0 )
    made = 0;
    saved = 0;
    y = 0.0;
    x = [];
    return
  elseif ( n == 0 )
    x = [];
    return
  end
%
%  Record the range of X we need to fill in.
%
  x_lo_index = 1;
  x_hi_index = n;
%
%  Use up the old value, if we have it.
%
  if ( saved == 1 )
    x(1) = y;
    saved = 0;
    x_lo_index = 2;
  end
%
%  Maybe we don't need any more values.
%
  if ( x_hi_index - x_lo_index + 1 == 0 )
%
%  If we need just one new value, do that here to avoid null arrays.
%
  elseif ( x_hi_index - x_lo_index + 1 == 1 )

    [ r(1), seed ] = r8_uniform_01 ( seed );

    if ( r(1) == 0.0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'R8VEC_NORMAL_01 - Fatal error!\n' );
      fprintf ( 1, '  R8_UNIFORM_01 returned a value of 0.\n' );
      error ( 'R8VEC_NORMAL_01 - Fatal error!' );
    end

    [ r(2), seed ] = r8_uniform_01 ( seed );

    x(x_hi_index) = ...
             sqrt ( - 2.0 * log ( r(1) ) ) * cos ( 2.0 * pi * r(2) );
    y =      sqrt ( - 2.0 * log ( r(1) ) ) * sin ( 2.0 * pi * r(2) );

    saved = 1;

    made = made + 2;
%
%  If we require an even number of values, that's easy.
%
  elseif ( mod ( x_hi_index - x_lo_index + 1, 2 ) == 0 )

    m = floor ( ( x_hi_index - x_lo_index + 1 ) / 2 );

    [ r, seed ] = r8vec_uniform_01 ( 2*m, seed );

    x(x_lo_index:2:x_hi_index-1) = ...
      sqrt ( -2.0 * log ( r(1:2:2*m-1) ) ) ...
      .* cos ( 2.0 * pi * r(2:2:2*m) );

    x(x_lo_index+1:2:x_hi_index) = ...
      sqrt ( -2.0 * log ( r(1:2:2*m-1) ) ) ...
      .* sin ( 2.0 * pi * r(2:2:2*m) );

    made = made + x_hi_index - x_lo_index + 1;
%
%  If we require an odd number of values, we generate an even number,
%  and handle the last pair specially, storing one in X(N), and
%  saving the other for later.
%
  else

    x_hi_index = x_hi_index - 1;

    m = floor ( ( x_hi_index - x_lo_index + 1 ) / 2 ) + 1;

    [ r, seed ] = r8vec_uniform_01 ( 2*m, seed );

    x(x_lo_index:2:x_hi_index-1) = ...
      sqrt ( -2.0 * log ( r(1:2:2*m-3) ) ) ...
      .* cos ( 2.0 * pi * r(2:2:2*m-2) );

    x(x_lo_index+1:2:x_hi_index) = ...
      sqrt ( -2.0 * log ( r(1:2:2*m-3) ) ) ...
      .* sin ( 2.0 * pi * r(2:2:2*m-2) );

    x(n) = sqrt ( -2.0 * log ( r(2*m-1) ) ) ...
      * cos ( 2.0 * pi * r(2*m) );

    y = sqrt ( -2.0 * log ( r(2*m-1) ) ) ...
      * sin ( 2.0 * pi * r(2*m) );

    saved = 1;

    made = made + x_hi_index - x_lo_index + 2;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%--- help for arima/simulate ---
%
% SIMULATE Simulate ARIMA model responses and conditional variances
%
%  Syntax:
%
%    [Y,E,V] = simulate(Mdl,numObs)
%    [Y,E,V] = simulate(Mdl,numObs,param1,val1,param2,val2,...)
%
%  Description:
%
%    Simulate sample paths of responses, innovations, and conditional
%    variances of a univariate ARIMA process.
%
%  Input Arguments:
%
%    Mdl - ARIMA model specification object, as produced by the ARIMA
%      constructor or ARIMA/ESTIMATE method.
%
%    numObs - Positive integer indicating the number of observations (rows)
%      generated for each path of the outputs Y, E, and V.
%
%  Optional Input Parameter Name/Value Pairs:
%
%    'NumPaths'   Positive integer indicating the number of sample paths
%                 (columns) generated for all simulated outputs. The default
%                 is 1.
%
%    'Y0'         Presample response data, providing initial values for the
%                 model. Y0 is a column vector or a matrix. If Y0 is a
%                 column vector, then it is applied to each simulated path.
%                 If Y0 is a matrix, then it must have at least NumPaths
%                 columns. Y0 may have any number of rows, provided at least
%                 Mdl.P observations exist to initialize the model. If the
%                 number of rows exceeds Mdl.P, then only the most recent
%                 Mdl.P observations are used. If the number of columns
%                 exceeds NumPaths, then only the first NumPaths columns are
%                 used. If Y0 is unspecified, any necessary presample
%                 observations are set to the unconditional mean for stationary
%                 AR processes, and to zero if the process is non-stationary
%                 or contains a regression component. The last row contains
%                 the most recent observation.
%
%    'E0'         Mean-zero presample innovations, providing initial values
%                 for the model. E0 is a column vector or a matrix. If E0 is
%                 a column vector, then it is applied to each simulated path.
%                 If E0 is a matrix, then it must have at least NumPaths
%                 columns. E0 may have any number of rows, provided sufficient
%                 observations exist to initialize the ARIMA model as well
%                 as any conditional variance model (the number of observations
%                 required is at least Mdl.Q, but may be more if a conditional
%                 variance model is included). If the number of rows
%                 exceeds the number necessary, then only the most recent
%                 observations are used. If the number of columns exceeds
%                 NumPaths, then only the first NumPaths columns are used. If
%                 no presample data is specified, any necessary observations
%                 are set to zero. The last row contains the most recent
%                 observation.
%
%    'V0'         Positive presample conditional variances, providing initial
%                 values for any conditional variance model; if the variance
%                 of the model is constant, then V0 is unnecessary. V0 is a
%                 column vector or a matrix. If V0 is a column vector, then
%                 it is applied to each simulated path. If V0 is a matrix,
%                 then it must have at least NumPaths columns. V0 may have
%                 any number of rows, provided sufficient observations exist
%                 to initialize any conditional variance model. If the number
%                 of rows exceeds the number necessary, then only the most
%                 recent observations are used. If the number of columns
%                 exceeds NumPaths, then only the first NumPaths columns are
%                 used. If no presample variance data is specified, any
%                 necessary observations are set to the unconditional variance
%                 of the conditional variance process. The last row contains
%                 the most recent observation.
%
%    'X'          Matrix of predictor data used to include a regression
%                 component in the conditional mean. Each column of X is a
%                 separate time series, and the last row of each contains
%                 the most recent observation of each series. The number of
%                 observations in X must equal or exceed numObs. When the
%                 number of observations in X exceeds numObs, only the most
%                 recent observations are used. If missing, the conditional
%                 mean will have no regression component regardless of the
%                 presence of any regression coefficients found in the model.
%
%  Output Arguments:
%
%    Y - numObs-by-NumPaths matrix of simulated response data.
%
%    E - numObs-by-NumPaths matrix of simulated mean-zero innovations.
%
%    V - numObs-by-NumPaths matrix of conditional variances of the
%      innovations in E.
%
%  Notes:
%
%    o Missing data values, indicated by NaNs, are removed from X by listwise
%      deletion (i.e., any row in X with at least one NaN is removed), reducing
%      the effective sample size. The presample data Y0, E0, and V0 are merged
%      into a composite series, and any row of the combined series with at
%      least one NaN is also removed by listwise deletion. The presample data
%      is also synchronized such that the last (most recent) observation of
%      each series occurs at the same time.
%
%   o  Regression models included in the conditional mean are based on the
%      presence of the predictor matrix X. Although each column of the output
%      time series represents a different path of the corresponding univariate
%      stochastic process, the regression matrix X represents as a single
%      path of a (possibly) multivariate time series matrix in which each
%      column is a different time series. When the conditional mean has a
%      regression component, the entire predictor matrix X is applied to
%      every column of the output time series.
%
%  References:
%
%    [1] Box, G. E. P., G. M. Jenkins, and G. C. Reinsel. Time Series
%        Analysis: Forecasting and Control. 3rd edition. Upper Saddle River,
%        NJ: Prentice-Hall, 1994.
%
%    [2] Enders, W. Applied Econometric Time Series. Hoboken, NJ: John Wiley
%        & Sons, 1995.
%
%    [3] Hamilton, J. D. Time Series Analysis. Princeton, NJ: Princeton
%        University Press, 1994.
%
%  See also ARIMA, FORECAST, ESTIMATE, INFER, FILTER.
%
%    Reference page in Doc Center
%       doc arima/simulate
%
%    Other functions named simulate
%
%       conjugateblm/simulate    generic/simulate
%       customblm/simulate       gjr/simulate
%       diffuseblm/simulate      regARIMA/simulate
%       dsge/simulate            sde/simulate
%       dtmc/simulate            semiconjugateblm/simulate
%       egarch/simulate          ssm/simulate
%       empiricalblm/simulate    varm/simulate
%       garch/simulate           vecm/simulate
%

author: jmaih


text:
function [h]= identify_classifier(h)
% load classifier files

rootS2p = which('run_pipeline');
if isempty(rootS2p)
    warndlg('Suite2p location not in PATH! where is run_pipeline.m?')
    [filename1,rootS2p]=uigetfile(root, 'Select Data File');

end
rootS2p = fileparts(rootS2p);
rootS2p = fullfile(rootS2p, 'configFiles');

flag=0
if isfield(h.dat.ops, 'classifier') && ~isempty(h.dat.ops.classifier)
    if ~exist(fullfile(rootS2p, h.dat.ops.classifier), 'file')
         warndlg('specified classifier database not found, reverting  to last used');
    else
        def_file = h.dat.ops.classifier;
        flag = 1;
        h.dat.cl.fpath  = fullfile(rootS2p, def_file);

        hload = load(h.dat.cl.fpath);
        if ~isfield(hload, 'st') || ~isfield(hload, 'statLabels') || ~isfield(hload, 'prior')
            error('found a non-classifier file in configFiles, called %s. \n Please remove and try again!', def_file)
        end

        h.st        = hload.st;
        h.prior     = hload.prior;
        h.statLabels = hload.statLabels;
    end
else
%     warning('no specified classifier database, reverting  to last used');
end
if (flag==0)
    fs = dir(fullfile(rootS2p, '*.mat'));
    if isempty(fs)
        warndlg('no classifier found, please make a new one!')
        h           = make_classifier(h);
    else
        [~, imax]    = max([fs.datenum]);

        def_file = fs(imax).name;

        h.dat.cl.fpath  = fullfile(rootS2p, def_file);

        hload = load(h.dat.cl.fpath);
        if ~isfield(hload, 'st') || ~isfield(hload, 'statLabels') || ~isfield(hload, 'prior')
            error('found a non-classifier file in configFiles, called %s. \n Please remove and try again!', def_file)
        end

        h.st        = hload.st;
        h.prior     = hload.prior;
        h.statLabels = hload.statLabels;
        h.dat.ops.classifier = def_file;
    end
end

end

function h = make_classifier(h)
% new classifier button
rootS2p = which('run_pipeline');
if isempty(rootS2p)
    warndlg('Suite2p location not in PATH! where is run_pipeline.m?')
    [filename1,rootS2p]=uigetfile(root, 'Select run_pipeline.m');

end
rootS2p = fileparts(rootS2p);
rootS2p = fullfile(rootS2p, 'configFiles');

def_name = fullfile(rootS2p, 'cl_new.mat');
[FileName,PathName] = uiputfile('*.mat', 'Create new classifier', def_name);

if FileName
    load(fullfile(rootS2p, 'priors', 'prior_default.mat'));
    st = [];
    save(fullfile(PathName, FileName), 'st', 'prior', 'statLabels')

    h.st = st;
    h.prior = prior;
    h.statLabels = statLabels;

    h.dat.cl.fpath          = fullfile(PathName, FileName);
    h                       = classROI(h);

%     h = splitROIleftright(h);
%     h = buildLambdaValue(h);
%     redraw_figure(h);

    hload = load(h.dat.cl.fpath);
    h.st        = hload.st;
    h.statLabels = hload.statLabels;

    h.prior     = hload.prior;
else
    error('you must make a new classifier first!')
end

end

author: cortex-lab


text:
function [hpb] = kW2hpb(kW)
% Convert power from kilowatts to boiler horsepower.
% Chad A. Greene 2012
hpb = kW*0.101941995;

author: Sable

repo: mcbench-benchmarks


text:
function len = strlen(str)
%  len = strlen(str)
% compute the # of characters in str (ignoring 0s at the end)

%
% strlen.m
%
% Original Author: Bruce Fischl
%
% Copyright © 2011 The General Hospital Corporation (Boston, MA) "MGH"
%
% Terms and conditions for use, reproduction, distribution and contribution
% are found in the 'FreeSurfer Software License Agreement' contained
% in the file 'LICENSE' found in the FreeSurfer distribution, and here:
%
% https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense
%
% Reporting: freesurfer@nmr.mgh.harvard.edu
%

len = length(str) ;
for i=length(str):-1:1
	if (str(i) ~= 0)
		break ;
	end

	len = len-1;
end


text:
function [y, x] = max(f, flag, dim)
%MAX   Maximum value of a CHEBFUN.
%   MAX(F) and MAX(F, 'global') return the maximum value of the CHEBFUN F.
%
%   [Y, X] = MAX(F) returns also a point X such that F(X) = Y.
%
%   [Y, X] = MAX(F, 'local') returns not just the global maximum value but all
%   of the local maxima.
%
%   If F is complex-valued, absolute values are taken to determine the maxima,
%   but the resulting values correspond to those of the original function.
%
%   If F is array-valued, then the columns of X and Y correspond to the columns
%   of F. NaNs are used to pad Y and X when the 'local' flag is used and the
%   columns are not of the same length.
%
%   H = MAX(F, G), where F and G are CHEBFUNs defined on the same domain,
%   returns a CHEBFUN H such that H(x) = max(F(x), G(x)) for all x in the
%   domain of F and G. Alternatively, either F or G may be a scalar.
%
%   MAX(F, [], DIM) computes the maximum of the CHEBFUN F in the dimension DIM.
%   If DIM = 1 and F is a column CHEBFUN or DIM = 2 and F is a row CHEBFUN, this
%   is equivalent to MAX(F). Otherwise, MAX(F, [], DIM) returns a CHEBFUN which
%   is the maximum across the discrete dimension of F. For example, if F is a
%   quasimatrix with two columns, MAX(F, [], 2) = MAX(F(:,1), F(:,2)).
%
% See also MIN, MINANDMAX, ROOTS.

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

% Trivial empty case.
if ( isempty(f) )
    x = [];
    y = [];
    return
end

if ( nargin == 3 )
    if ( ~any(dim == [1 2]) )
        error('CHEBFUN:CHEBFUN:max:badDim', ...
            'DIM input to CHEBFUN MAX must be 1 or 2.');
    end

    if ( dim ~= 1 + f(1).isTransposed )
        % Take max across discrete dimension of a quasimatrix:
        f = cheb2cell(f);
        y = -realmax;
        for k = 1:numel(f)
            y = maxOfTwoChebfuns(y, f{k});
        end
        y = merge(y);
        return
    end
end

if ( (nargin > 2) && isempty(flag) )
    % MAX(F, [], 1).
    flag = 'global';
end

if ( (nargin == 1) || strcmp(flag, 'global') )
    % MAX(F) or MAX(F, 'global')
    [y, x] = globalMax(f);

elseif ( isa(flag, 'chebfun') || isnumeric(flag) )
    % MAX(F, G)
    y = maxOfTwoChebfuns(f, flag);

elseif ( strcmp(flag, 'local') )
    % MAX(F, 'local')
    [y, x] = localMax(f);

else
    error('CHEBFUN:CHEBFUN:max:flag', 'Unrecognized flag.');

end

end

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%% GLOBAL MAXIMUM %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [y, x] = globalMax(f)

% Call MINANDMAX():
[y, x] = minandmax(f);

% Extract the maximum:
y = y(2,:);
x = x(2,:);

end

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%% LOCAL MAXIMA %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [y, x] = localMax(f)

% Call MINANDMAX():
[y, x] = minandmax(f, 'local');

% Determine which are maxima.

ends = f(1).domain([1, end]).'; % Endpoints of the domain are special.
f = mat2cell(f); % Convert f into a cell of scalar-valued CHEBFUNs.

% Loop over the FUNs:
for k = 1:numel(f)
    % Compute 1st and 2nd derivatives.
    dfk = diff(f{k});
    dfk2 = diff(dfk);

    % For interior extrema, look at 2nd derivative:
    maximaLoc = feval(diff(f{k}, 2), x(:,k)) < 0;

    % For end-points, look at 1st derivative:
    dfk_ends = feval(dfk, ends);
    endptMaxLoc = dfk_ends.*[1, -1]' < 0;

    % If 1st derivative is small at an endpoint, assume it's zero and try to
    % use 2nd derivative to determine if it's a minimum:
    %
    % [TODO]:  What if the 2nd derivative is zero, so that rounding error
    % precludes us from accurately determining the sign?
    smallEndDer = abs(dfk_ends) < 1e3*vscale(dfk)*eps;
    endptMaxLoc(smallEndDer) = feval(dfk2, ends(smallEndDer)) < 0;

    maximaLoc(1) = endptMaxLoc(1);
    maximaLoc(x(:,k) == ends(2)) = endptMaxLoc(2);

    % Set points corresponding to local maxima to NaN:
    y(~maximaLoc,k) = NaN;
    x(~maximaLoc,k) = NaN;

    % Sort the result
    [x(:,k), maximaLoc] = sort(x(:,k));
    y(:,k) = y(maximaLoc,k);
end

% Remove any rows which contain only NaNs.
x(all(isnan(x), 2),:) = [];
y(all(isnan(y), 2),:) = [];

end

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%% MAX(F, G) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function h = maxOfTwoChebfuns(f, g)
% Return the function h(x) = max(f(x), g(x)) for all x.

% If one is complex, use abs(f) and abs(g) to determine which function values to
% keep. (experimental feature)
if ( isreal(f) && isreal(g) && (nargin < 3) )
	S = sign(f - g);
else
	S = sign(abs(f) - abs(g));
end

% Heaviside function (0 where f > g, 1 where f < g);
H = 0.5*(S + 1);
notH = 0.5*(1 - S); % ~H.

% Combine for output:
h = H.*f + notH.*g;

% [TODO]: Enforce continuity?

% TODO: Why simplify?
% Simplify:
h = simplify(h);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [CT] = calculate_CT(Seconds, Throttle, RPM, b, h, gramsMeas, varargin)
% This function calculates the parameter for a zero intercept linear fit
% between RPM^2 and thrust of the form:
%   Thrust = (RPM^2)*CT
% Output is given in the form of returned CT value and plots.
%
% Argument definitions:
%   b = length of test rig from pin joint to scale contact (same units as h)
%   h = length of test rig from pin joint to center motor axis (same units as b)
%   gramsMeas = user measured values direct from scale, MUST BE IN GRAMS!
%   These values must correspond to ALL unique throttle values, even those
%   outside of "LowCutoff" and "HighCutoff."
% (Note: values of gramsMeas corresponding to throttle setting outside
% "LowCutoff" and "HighCutoff" are not used in calculation, therefore if
% measurements are only available for the range of "LowCutoff" to
% "HighCutoff", the vector can be padded with zeros in the appropriate
% locations.
%
% See also calculate_CR_B, calculate_CQ.

% Assign variable length input arguments
if ((nargin-5)==0) % if no input arguments are given for cutoff range:
LowCutoff = min(Throttle); % Use full range of input throttle data
HighCutoff = max(Throttle);
else
LowCutoff = varargin{1}; % Ignore data below this throttle value
HighCutoff = varargin{2}; % Ignore data above this throttle value
end

% Find unique Throttle Setting values (Not counting 0)
ThrottleU = unique(Throttle(Throttle~=0));

% Check length of inputs match
lengthGM = length(gramsMeas);
lengthTU = length(ThrottleU);
if lengthGM ~= lengthTU
    msgbox({['Length of measured grams values vector (' num2str(lengthGM) ')'];
        [' and number of unique Throttle values (' num2str(lengthTU) ')']; ...
       ' do not match! Check input data!'},'Argument Error','warn');
    return
end

% Find unique/filtered Throttle Setting values between low/high cutoff limits
ThrottleUF = ThrottleU(ThrottleU>=LowCutoff & ThrottleU<=HighCutoff);

% Find coresponding gramsMeasurments
gramsMeasU = gramsMeas(ThrottleU>=LowCutoff & ThrottleU<=HighCutoff); % Unique scale values

% Calculate arms ratio
armRatio = b/h;
g = 9.81; % Gravity (m/s^2)

% Calculate thrust generated by motors
newtonsThrust = g*armRatio*gramsMeasU/1000; % Vector of thrusts in Newtons

% Find an average RPM value for each unique/filtered Throttle Setting
RPMAverages100 = zeros(length(ThrottleUF),1); % preallocate
for i = 1:length(ThrottleUF) % For all of the unique and filtered throttle settings:
    RPMUF = RPM(Throttle==ThrottleUF(i));  % get RPM values cooresponding to unique/filtered Throttle value
    if length(RPMUF)>=120 % check to make sure there are enough data points for filtering
        RPMfilt100 = RPMUF(100:end); % RPMUF minus the first 100 data points (to remove transient RPM fluctuation as steady state value is reached)
    else RPMfilt100 = RPMUF((fix(length(RPMUF)/2)):end); % If there's a limited number of points, just take the second half of the values as steady state
    end
    RPMAverages100(i,1) = mean(RPMfilt100); % Take the average RPM value from these filtered data points
end

RPMsq = RPMAverages100.^2;
CT = RPMsq\newtonsThrust; % This is a technique for forcing a zero intercept using matrix multiplication techniques
CTp = [CT; 0];

% Plot the linear relationship and display calculated CT value
figure
plot(RPMsq,newtonsThrust,'*',RPMsq,polyval(CTp,RPMsq))
xlabel('RPM^2 (rev^2/min^2)','FontSize',12)
ylabel('Thrust (N)','FontSize',12)
title({['Linear Fit of Thrust vs. RPM^2']; ...
    [' '];['CT: ' num2str(CT) ' (N/RPM^2)']},'FontSize',14)
xlims = get(gca,'xlim');
ylims = get(gca,'ylim');
grid on

% Plot Residuals for assesing general quality of fit
predictedThrust = polyval(CTp, RPMsq); % predict thrust for each throttle setting of interest
residuals = newtonsThrust - predictedThrust; % find absolute error between prediction and actual
figure
plot(predictedThrust,residuals,'r*',[min(predictedThrust) max(predictedThrust)],[0 0])
title('CT Residuals Plot','FontSize',14)
xlabel('Predicted Thrust (N)','FontSize',12)
ylabel('Prediction Error (N)','FontSize',12)
grid on
axis tight
end

repo: Quad-Sim

path: github-repos/MATLAB/dch33-Quad-Sim/Quad-Sim-961bc69d4939c8d0661eeb9820de668994262f65/Quadcopter Dynamic Modeling and Simulation/Data Acquisition & Analysis/MATLAB Data Analysis/calculate_CT.m


text:
function a = dec2bin(d)
i=1;
a=zeros(1,65535);
while d >= 2
    r=rem(d,2);
    if r==1
        a(i)=1;
    else
        a(i)=0;
    end
    i=i+1;
    d=floor(d/2);
end
if d == 2
    a(i) = 0;
else
    a(i) = 1;
end
x=[a(16) a(15) a(14) a(13) a(12) a(11) a(10) a(9) a(8) a(7) a(6) a(5) a(4) a(3) a(2) a(1)];


text:
function y = ssrfunc(x)

% sun-synchronous repeating ground track
% orbit objective function

% required by ssrepeat.m

% input

%  x = array of current dependent variables

% output

%  y = function value array evaluated at x

% Orbital Mechanics with MATLAB

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global j2 mu req omega xldot

global ecc xqp

y = zeros(2, 1);

% "unload" current solution

sma = x(1);
inc = x(2);

% determine current values of nonlinear system

mm = sqrt(mu / sma ^ 3);

slr = sma * (1.0 - ecc * ecc);

cinc = cos(inc);

% compute perturbations

raandot = -1.5 * j2 * mm * (req / slr) ^ 2 * cinc;

apdot = 0.75 * j2 * mm * (req / slr) ^ 2 * (5.0 * cinc ^ 2 - 1);

madot = mm + 0.75 * j2 * mm * (req / sma) ^ 2 * (3.0 * cinc ^ 2 - 1) ...
    / (1.0 - ecc * ecc) ^ 1.5;

% sun-synchronous inclination equation

y(1) = cinc + 2.0 * xldot * sma ^ 3.5 * (1 - ecc * ecc) ^ 2 ...
    / (3 * j2 * req * req * sqrt(mu));

% repeating groundtrack equation

y(2) = (1.0 / xqp) / (omega - raandot) - (1.0 / (apdot + madot));


text:
% a=0.005, x^2
% a=0.0002, abs(x)^3
% a=0.0001, x^4

a = 0.005;
b = 0.00025;
xs = 0:50;
ys2 = zeros(length(xs), 1);
ys3 = zeros(length(xs), 1);
for i = 1:length(xs)
    ys2(i) = exp(-a*abs(xs(i))^2);
    ys3(i) = exp(-b*abs(xs(i))^3);
end
plot(xs, ys2);
hold on;
plot(xs, ys3);
hold off;

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
% add activation of all previous layers, which are required to have the
% same size.
%
classdef AddNode < GraphNode

    methods
        function obj = AddNode(dimOut)
            obj = obj@GraphNode('Add',dimOut);
        end

        function obj = forward(obj,prev_layers)
            obj = obj.preprocessingForward(prev_layers);
            obj.a = prev_layers{1}.a;
            for i=2:length(prev_layers)
                obj.a = obj.a + prev_layers{i}.a;
            end
            obj = forward@GraphNode(obj, prev_layers);
        end

        function obj = backward(obj,prev_layers, future_layers)
            if obj.skipGrad || obj.skipBP
                return;
            end

            future_grad = obj.GetFutureGrad(future_layers);
            for i=1:length(prev_layers)
                obj.grad{i} = future_grad;
            end

            obj = backward@GraphNode(obj, prev_layers, future_layers);
        end

    end

end

author: singaxiong

repo: SignalGraph


text:
function [uOutput] = scecdcimp(nFunction, sFilename)
% import UW catalog

% Filter function switchyard
if nFunction == FilterOp.getDescription
    uOutput = 'UW catalog ';
elseif nFunction == FilterOp.importCatalog
    % Read formated data
    mData = textread(sFilename, '%s', 'delimiter', '\n', 'whitespace', '');
    % Create empty catalog
    uOutput = zeros(length(mData), 9);
    % Loop thru all lines of catalog and convert them
    for i = 1:length(mData)
        if rem(i,100) == 0 ; disp([ num2str(i) ' of ' num2str(length(mData)) ' events processed ']); end
        try
            l = find(mData{i} == ' ');
            mData{i}(l) = '0';
            uOutput(i,1) = -str2num(mData{i}(30:32)) - str2num(mData{i}(34:37))/6000 ;  % Longitude
            if mData{i}(33) == 'E' ; uOutput(i,1)  = -1*uOutput(i,1); end
            uOutput(i,2) =  str2num(mData{i}(22:23)) + str2num(mData{i}(25:28))/6000 ;  % Latitude
            if mData{i}(24) == 'S' ; uOutput(i,1)  = -1*uOutput(i,1); end
            uOutput(i,3) = str2num(mData{i}(3:6));    % Year
            uOutput(i,4) = str2num(mData{i}(7:8));    % Month
            uOutput(i,5) = str2num(mData{i}(9:10));   % Day
            uOutput(i,6) = str2num(mData{i}(45:48));  % Magnitude
            uOutput(i,7) = str2num(mData{i}(38:43));  % Depth
            uOutput(i,8) = str2num(mData{i}(11:12));  % Hour
            uOutput(i,9) = str2num(mData{i}(13:14));  % Minute

            %Create decimal year
            uOutput(i,3) = decyear([uOutput(i,3) uOutput(i,4) uOutput(i,5) uOutput(i,8) uOutput(i,9)]);

            if mData{i}(2) == 'P'  ||  mData{i}(2) == 'X'  ||  mData{i}(2) == 'L'
                uOutput(i,:)=nan;
            end

        catch
            msg.dbfprintf('Import: Problem in line %d of %s. Line ignored.\n',i, sFilename);  % eliminate explosions
            uOutput(i,:)=nan;
        end
    end
    l = isnan(uOutput(:,1));
    uOutput(l,:) = [];
end

author: CelsoReyes


text:
function[]=jlab_addpath
%JLAB_ADDPATH  Adds JLAB and JDATA subdirectories to your Matlab search path.
%
%   This script adds JLAB subdirectories to your Matlab search path, as
%   well as the JDATA directory if it is found to exist.
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2015--2019 J.M. Lilly --- type 'help jlab_license' for details

%First, find the path to JLAB
fullpath=which('jlab_license');
jlab_path=fullpath(1:end-15);
bool=real(jlab_path)==real('\');
jlab_path(bool)='/';

%jlab_path

for k=1:2
    if k==1
        disp('Adding JLAB subdirectories to your Matlab search path.')
        if exist(jlab_path)
            jlab_addpath_subdirectories(jlab_path);
        end
    end
    if k==2
        jlab_path=[jlab_path(1:end-5) '/jdata'];
        if exist(jlab_path)
            disp('Looks like you have JDATA installed; adding this to your path also.')
            addpath(jlab_path)
            jlab_addpath_subdirectories(jlab_path);
        end
    end
end
clear i j k dirlist dirlisti jlab_path fullpath bool

function[]=jlab_addpath_subdirectories(path)
dirlist=dir(path);
for i=1:length(dirlist)
    if dirlist(i).isdir
        if ~strcmpi(dirlist(i).name(1),'.')
            %dirlist(i).name
            addpath([path '/' dirlist(i).name])
            dirlisti=dir([path '/' dirlist(i).name]);
            for j=1:length(dirlisti)
                if dirlisti(j).isdir
                    if ~strcmpi(dirlisti(j).name(1),'.')
                        %[path '/' dirlist(i).name '/' dirlisti(j).name]
                        addpath([path '/' dirlist(i).name '/' dirlisti(j).name]);
                    end
                end
            end
        end
    end
end

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function [ratio, ratiominus, ratioplus, proba_post, yt] = classification_evaluation(bnet, BDT, class)
% Computes the classification ratio of a bnet structure on a test dataset BDT
% [ratio ratiominus ratioplus] = classification_evaluation(bnet, BDT, class)
%
% [ratiominus rationplus] is the 95 percent confident interval.
% results are in percentage [0 100].
%
%  francois.olivier.c.h@gmail.com
%

    [proba_post,engine] = inference(bnet, mat_to_bnt(BDT), class);
               [tmp yt] = max(proba_post, [],2);
                  [N L] = size(BDT);
                  count = length(find(BDT(class,:)==yt'));
                  ratio = 100*count/L;
[ratiominus, ratioplus] = confiance(ratio,L);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [I, J] = confiance(t, N)
% Compute the 95 percent confident interval
%
% see Y. Bennani and F. Bossaert,
%     Predictive neural networks for traffic disturbance detection in the telephone network
%     In Proceedings of IMACS-CESA 1996, Lille, France.

Z    = 1.96; % this value for the 95 percent confident interval
T    = t/100;
tmp  = (Z*Z)/N;
D    = 1+tmp;
N1   = T+tmp/2;
tmp2 = T*(1-T)/N + tmp/(4*N);
N2   = Z*sqrt(tmp2);
I    = 100*(N1-N2)/D;
J    = 100*(N1+N2)/D;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function ring_render(ring_param,x_list,y_list)

%Extract ring parameters:
%------------------------
Ri=ring_param{2}; %inner radius
Ro=ring_param{3}; %outer modulation size
H=ring_param{4}; %height
Rfreq=ring_param{1}; %modulation frequency
face_color=ring_param{5}/255; %color

delta_theta=5; %angle spacing (smaller values will increase render time)
theta=0:delta_theta:360;

%Complete the list of cross-section nodes:
%-----------------------------------------
x=x_list;
y=y_list;
y_new=[x(1:length(x)-1) x(1)];
z_new=[y(1:length(y)-1) y(1)];
x_new=zeros(size(y));
x=x_new;
y=y_new+Ri;
z=z_new;
clear x_new y_new z_new

%Render ring:
%------------
for nt=1:length(theta)-1 %loop through each face angle

    th1=theta(nt);
    th2=theta(nt+1);

    for np=1:length(x)-1 %generate the edges of the face

        xtemp=x(np:np+1);
        ytemp=y(np:np+1);
        ztemp=z(np:np+1);

        %apply modulation to outer cross-section coordinates:
        r_temp1=(ytemp>Ri)*Ro*((sind(th1*Rfreq)+1)/2);
        r_temp2=(ytemp>Ri)*Ro*((sind(th2*Rfreq)+1)/2);

        %apply rotation transformation to cross-section coordinate:
        xtemp_rot1=xtemp*cosd(th1)+(ytemp+r_temp1)*sind(th1);
        xtemp_rot2=xtemp*cosd(th2)+(ytemp+r_temp2)*sind(th2);
        ytemp_rot1=xtemp*-sind(th1)+(ytemp+r_temp1)*cosd(th1);
        ytemp_rot2=xtemp*-sind(th2)+(ytemp+r_temp2)*cosd(th2);

        %edge lists:
        xtemp_face=[xtemp_rot1(1) xtemp_rot2(1) xtemp_rot2(2) xtemp_rot1(2) xtemp_rot1(1)];
        ytemp_face=[ytemp_rot1(1) ytemp_rot2(1) ytemp_rot2(2) ytemp_rot1(2) ytemp_rot1(1)];
        ztemp_face=[ztemp(1) ztemp(1) ztemp(2) ztemp(2) ztemp(1)];

        %render face:
        hf=fill3(xtemp_face,ytemp_face,ztemp_face,face_color);
        set(hf,'EdgeAlpha',0,'FaceLighting','phong','BackFaceLighting','lit','SpecularStrength',1,'FaceAlpha',.5)
        material shiny
        hold on
    end
end
axis equal
axis([-1 1 -1 1 -1 1]*(Ri+H))
axis off
set(gcf,'Color','k')
light('Position',[1 1 1]*25,'Style','local')
light('Position',[-1 1 -1]*25,'Style','local')
light('Position',[1 -1 -1]*25,'Style','local')
light('Position',[-1 -1 1]*25,'Style','local')


text:
function sT = som_table_modify(sT,action,arg1,arg2,arg3)

%SOM_TABLE_MODIFY Modify table: add or remove columns or rows.
%
% sTable = som_table_modify(sTable,action,arg1,[arg2],[arg3])
%
%  Input and output arguments ([]'s are optional):
%   sTable      (struct) table struct
%   action      (string) action id (see below).
%   arg1        (varies) Depending on action, 1 to 3 arguments
%   [arg2]      (varies) are needed. See below.
%   [arg3]      (varies)
%
%   sTable      (struct) the modified table struct
%
%  Actions and their arguments:
%   'addcol'    Add one or several new columns.
%               arg1 (cell array) new values
%                    (char)       new values (a single column can be given)
%                    (matrix)     new values
%               arg2 (cell array) new headers
%               arg3 (scalar)     at which position the new columns
%                                 should be inserted (at the end by default)
%   'addrow'    Add one or several new rows.
%               arg1 (cell array) new values
%                    (char)       new values (a single row can be given)
%                    (matrix)     new values
%               arg2 (scalar)     at which position the new rows
%                                 should be inserted (at the end by default)
%   'removecol' Remove one or several columns.
%               arg1 (vector)     indeces of columns to be removed
%   'removerow' Remove one or several rows.
%               arg1 (vector)     indeces of rows to be removed
%
% See also  SOM_TABLE_STRUCT, SOM_TABLE_PRINT.

% Contributed to SOM Toolbox 2.0, January 4th, 2002 by Juha Vesanto
% Copyright (c) by Juha Vesanto
% http://www.cis.hut.fi/projects/somtoolbox/

% Version 2.0beta juuso 040102

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[nrT,ncT] = size(sT.values);

switch action,
 case 'addcol',
  values = arg1;
  if ischar(values), values = cellstr(values); end
  if isnumeric(values), values = num2cell(values); end
  spans  = ones([size(values) 2]);
  [nr,nc] = size(values);
  if nargin<4, header = cell(1,nc); header(:) = {''}; else header = arg2; end
  if ischar(header), header = cellstr(header); end
  if nargin<5, where = ncT+1; else  where  = arg3; end
  if nrT ~= nr,
    error('Mismatch between sizes of given table and additional columns')
  else
    sT.headers = [sT.headers(:,1:where-1), header, sT.headers(:,where:end)];
    sT.values  = [sT.values(:,1:where-1), values, sT.values(:,where:end)];
    sT.span    = [sT.span(:,1:where-1,:), spans, sT.span(:,where:end,:)];
  end
 case 'addrow',
  values  = arg1;
  if ischar(values), values = cellstr(values); end
  if isnumeric(values), values = num2cell(values); end
  [nr,nc] = size(values);
  spans   = ones([size(values) 2]);
  if nargin<4, where = nrT+1; else where  = arg2; end
  if ncT ~= nc,
    error('Mismatch between sizes of given table and additional rows')
  else
    sT.values = [sT.values(1:where-1,:); values; sT.values(where:end,:)];
    sT.span   = [sT.span(1:where-1,:,:); spans; sT.span(where:end,:,:)];
  end
 case 'removecol',
  where      = setdiff(1:ncT,arg1);
  sT.values  = sT.values(:,where);
  sT.headers = sT.headers(:,where);
  sT.span    = sT.span(:,where,:);
 case 'removerow',
  where      = setdiff(1:nrT,arg1);
  sT.values  = sT.values(where,:);
  sT.headers = sT.headers(where,:);
  sT.span    = sT.span(where,:,:);
end

return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: ilarinieminen

repo: SOM-Toolbox


text:
function [g1, g2] = simwhiteXrbfwhiteKernGradient(simKern, rbfKern, t1, varargin)

% SIMWHITEXRBFWHITEKERNGRADIENT Compute a cross gradient between a SIM-WHITE
% and an RBF-WHITE kernels.
% FORMAT
% DESC computes cross gradient of parameters of a cross kernel between a
% SIM-WHITE and an RBF-WHITE kernels for the multiple output kernel.
% ARG simKern : the kernel structure associated with the SIM-WHITE kernel.
% ARG rbfKern : the kernel structure associated with the RBF-WHITE kernel.
% ARG t1 : inputs for which kernel is to be computed.
% ARG covGrad : gradient of the objective function with respect to
% the elements of the cross kernel matrix.
% RETURN g1 : gradient of the parameters of the SIM-WHITE kernel, for
% ordering see simwhiteKernExtractParam.
% RETURN g2 : gradient of the parameters of the RBF-WHITE kernel, for
% ordering see rbfwhiteKernExtractParam.
%
% FORMAT
% DESC computes cross kernel terms between a SIM-WHITE and an RBF-WHITE
% kernels for the multiple output kernel.
% ARG simKern : the kernel structure associated with the SIM-WHITE kernel.
% ARG rbfKern : the kernel structure associated with the RBF-WHITE kernel.
% ARG t1 : row inputs for which kernel is to be computed.
% ARG t2 : column inputs for which kernel is to be computed.
% ARG covGrad : gradient of the objective function with respect to
% the elements of the cross kernel matrix.
% RETURN g1 : gradient of the parameters of the SIM-WHITE kernel, for
% ordering see simwhiteKernExtractParam.
% RETURN g2 : gradient of the parameters of the RBF-WHITE kernel, for
% ordering see rbfwhiteKernExtractParam.
%
% SEEALSO : multiKernParamInit, multiKernCompute, simwhiteKernParamInit,
% rbfwhiteKernParamInit, simwhiteKernExtractParam, rbfwhiteKernExtractParam
%
% COPYRIGHT : David Luengo, 2009

% KERN

if nargin < 5
    t2 = t1;
else
    t2 = varargin{1};
end
covGrad = varargin{end};

if size(t1, 2) > 1 | size(t2, 2) > 1
  error('Input can only have one column');
end
if simKern.variance ~= rbfKern.variance
  error('Kernels cannot be cross combined if they have different variances.')
end
if simKern.isStationary ~= rbfKern.isStationary
  error('Stationary and non-stationary kernels cannot be cross combined.')
end

g1 = zeros(1, simKern.nParams);
g2 = zeros(1, rbfKern.nParams);

T1 = repmat(t1, 1, size(t2, 1));
T2 = repmat(t2.', size(t1, 1), 1);
deltaT = T1 - T2;
indT = double(deltaT >= 0);

% Parameters required for further computations
isStationary = simKern.isStationary;
variance = simKern.variance;
decay = simKern.decay;
sensitivity = simKern.sensitivity;
invWidth = rbfKern.inverseWidth;

% Computing a normalised (i.e. variance = 1 and sensitivity = 1) kernel
simKern.variance = 1;
simKern.sensitivity = 1;
rbfKern.variance = 1;
K = simwhiteXrbfwhiteKernCompute(simKern, rbfKern, t1, t2);

% Gradient w.r.t. the decay (simKern)
g1(1) = variance * sensitivity * sum(sum( ((-deltaT+decay/invWidth) .* K  ...
    + 1/sqrt(2*pi*invWidth) ...
        * exp(-decay*deltaT + 0.5*(decay^2)/invWidth) ...
        .* (exp(-0.5*invWidth*((T2+decay/invWidth).^2)) ...
        - exp(-0.5*invWidth*((abs(deltaT).*(1-indT)+decay/invWidth).^2)))) ...
    .* covGrad));

% Gradient w.r.t. the inverse width (rbfKern)
g2(1) = variance * sensitivity * sum(sum( (-0.5*((decay/invWidth)^2) .* K  ...
    + 1/sqrt(8*pi*invWidth) ...
        * exp(-decay*deltaT + 0.5*(decay^2)/invWidth) ...
        .* ((T2-decay/invWidth).*exp(-0.5*invWidth*((T2+decay/invWidth).^2)) ...
        - (abs(deltaT).*(1-indT)-decay/invWidth) ...
            .* exp(-0.5*invWidth*((abs(deltaT).*(1-indT)+decay/invWidth).^2)))) ...
    .* covGrad));

% Gradient w.r.t. sigma_r^2
g1(2) = sensitivity * sum(sum(K .* covGrad));
g2(2) = 0; % Otherwise it is counted twice

% Gradient w.r.t. sensitivity (only simKern)
g1(3) = variance * sum(sum(K .* covGrad));

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
%[2015]-"Multi-verse optimizer: A nature-inspired algorithm for global
%optimization"

% (9/12/2020)

function MVO = jMultiVerseOptimizer(feat,label,opts)
% Parameters
lb    = 0;
ub    = 1;
thres = 0.5;
p     = 6;      % control TDR
Wmax  = 1;      % maximum WEP
Wmin  = 0.2;    % minimum WEP
type  = 1;

if isfield(opts,'T'), max_Iter = opts.T; end
if isfield(opts,'N'), N = opts.N; end
if isfield(opts,'p'), p = opts.p; end
if isfield(opts,'Wmin'), Wmin = opts.Wmin; end
if isfield(opts,'Wmax'), Wmax = opts.Wmax; end
if isfield(opts,'ty'), type = opts.ty; end
if isfield(opts,'thres'), thres = opts.thres; end

% Objective function
fun = @jFitnessFunction;
% Number of dimensions
dim = size(feat,2);
% Initial
X   = zeros(N,dim);
for i = 1:N
  for d = 1:dim
    X(i,d) = lb + (ub - lb) * rand();
  end
end
% Pre
fit  = zeros(1,N);
fitG = inf;

curve = inf;
t = 1;
% Iterations
while t <= max_Iter
  % Calculate inflation rate
  for i = 1:N
    fit(i) = fun(feat,label,(X(i,:) > thres),opts);
    % Best universe
    if fit(i) < fitG
      fitG = fit(i);
      Xgb  = X(i,:);
    end
  end
  % Sort universe from best to worst
  [fitSU, idx] = sort(fit,'ascend');
  X_SU         = X(idx,:);
  % Elitism (first 1 is elite)
  X(1,:) = X_SU(1,:);
  % Either 1-norm or 2-norm
  if type == 1
    % Normalize inflation rate using 2-norm
    NI = fitSU ./ sqrt(sum(fitSU .^ 2));
  elseif type == 2
    % Normalize inflation rate using 1-norm
    NI = fitSU / sum(fitSU);
  end
  % Normalize inverse inflation rate using 1-norm
  inv_fitSU = 1 ./ (1 + fitSU);
  inv_NI    = inv_fitSU / sum(inv_fitSU);
  % Wormhole Existence probability (3.3), increases from 0.2 to 1
  WEP = Wmin + t * ((Wmax - Wmin) / max_Iter);
  % Travelling disrance rate (3.4), descreases from 0.6 to 0
  TDR = 1 - ((t ^ (1 / p)) / (max_Iter ^ (1 / p)));
  % Start with 2 since first is elite
  for i = 2:N
    % Define black hole
    idx_BH = i;
    for d = 1:dim
      % White/black hole tunnels & exchange object of universes (3.1)
      r1 = rand();
      if r1 < NI(i)
        % Random select k with roulette wheel
        idx_WH       = jRouletteWheelSelection(inv_NI);
        % Position update
        X(idx_BH, d) = X_SU(idx_WH, d);
      end
      % Local changes for universes (3.2)
      r2 = rand();
      if r2 < WEP
        r3 = rand();
        r4 = rand();
        if r3 < 0.5
          X(i,d) = Xgb(d) + TDR * ((ub - lb) * r4 + lb);
        else
          X(i,d) = Xgb(d) - TDR * ((ub - lb) * r4 + lb);
        end
      else
        X(i,d) = X(i,d);
      end
    end
    % Boundary
    XB = X(i,:); XB(XB > ub) = ub; XB(XB < lb) = lb;
    X(i,:) = XB;
  end
  curve(t) = fitG;
  fprintf('\nIteration %d Best (MVO)= %f',t,curve(t))
  t = t + 1;
end
% Select features
Pos   = 1:dim;
Sf    = Pos((Xgb > thres) == 1);
sFeat = feat(:,Sf);
% Store results
MVO.sf = Sf;
MVO.ff = sFeat;
MVO.nf = length(Sf);
MVO.c  = curve;
MVO.f  = feat;
MVO.l  = label;
end

%// Roulette Wheel Selection //
function Index = jRouletteWheelSelection(prob)
% Cummulative summation
C = cumsum(prob);
% Random one value, most probability value [0~1]
P = rand();
% Route wheel
for i = 1:length(C)
	if C(i) > P
    Index = i;
    break;
  end
end
end

author: JingweiToo

repo: Wrapper-Feature-Selection-Toolbox


text:
function obj = updateHeatmapAnnotation(obj,anIndex)

%-------X/YLABEL FIELDS--------%
% title...[DONE]
% titlefont.size...[DONE]
% titlefont.family...[DONE]
% titlefont.color...[DONE]

%------ANNOTATION FIELDS-------%

% x: ...[DONE]
% y: ...[DONE]
% xref: ...[DONE]
% yref: ...[DONE]
% text: ...[DONE]
% showarrow: ...[HANDLED BY CALL TO ANNOTATION];
% font: ...[DONE]
% xanchor: ...[DONE]
% yanchor: ...[DONE]
% align: ...[DONE]
% arrowhead: ...[HANDLED BY CALL FROM ANNOTATION];
% arrowsize: ...[HANDLED BY CALL FROM ANNOTATION];
% arrowwidth: ...[HANDLED BY CALL FROM ANNOTATION];
% arrowcolor: ...[HANDLED BY CALL FROM ANNOTATION];
% ax: ...[HANDLED BY CALL FROM ANNOTATION];
% ay: ...[HANDLED BY CALL FROM ANNOTATION];
% textangle: ...[DONE]
% bordercolor: ...[DONE]
% borderwidth: ...[DONE]
% borderpad: ...[DONE]
% bgcolor: ...[DONE]
% opacity: ...[NOT SUPPORTED IN MATLAB]

%-AXIS INDEX-%
nanns = length(obj.layout.annotations);
axIndex = nanns + obj.getAxisIndex(obj.State.Text(anIndex).AssociatedAxis);

%-CHECK FOR MULTIPLE AXES-%
[xsource, ysource] = findSourceAxis(obj,anIndex);

%-get heatmap title name-%
title_name = obj.State.Text(anIndex).Handle;

%-show arrow-%
obj.layout.annotations{axIndex}.showarrow = false;

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

%-anchor title to paper-%
if obj.State.Text(anIndex).Title
    %-xref-%
    obj.layout.annotations{axIndex}.xref = 'paper';
    %-yref-%
    obj.layout.annotations{axIndex}.yref = 'paper';
else
    %-xref-%
    obj.layout.annotations{axIndex}.xref = ['x' num2str(xsource)];
    %-yref-%
    obj.layout.annotations{axIndex}.yref = ['y' num2str(ysource)];
end

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

%-xanchor-%
obj.layout.annotations{axIndex}.xanchor = 'middle';

%-align-%
obj.layout.annotations{axIndex}.align = 'middle';

%-xanchor-%
obj.layout.annotations{axIndex}.yanchor = 'top';

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

%-text-%
obj.layout.annotations{axIndex}.text = sprintf('<b>%s</b>', title_name);
obj.layout.annotations{axIndex}.font.size = 14;

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

if obj.State.Text(anIndex).Title

    %-AXIS DATA-%
    eval(['xaxis = obj.layout.xaxis' num2str(xsource) ';']);
    eval(['yaxis = obj.layout.yaxis' num2str(ysource) ';']);

    %-x position-%
    obj.layout.annotations{axIndex}.x = mean(xaxis.domain);
    %-y position-%
    obj.layout.annotations{axIndex}.y = (yaxis.domain(2) + 0.04);
else
    %-x position-%
    obj.layout.annotations{axIndex}.x = text_data.Position(1);
    %-y position-%
    obj.layout.annotations{axIndex}.y = text_data.Position(2);
end

%-------------------------------------------------------------------------%
end

author: plotly


text:
function data_num = p04_data_num ( )

%*****************************************************************************80
%
%% P04_DATA_NUM returns the dimension of the data vector for problem 4.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 August 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer DATA_NUM, the dimension of the data vector.
%
  data_num = 13;

  return
end


text:
function result = sphere01_triangle_qua02 ( v1, v2, v3, f )

%*****************************************************************************80
%
%% SPHERE01_TRIANGLE_QUAD02 quadrature over a triangle on the unit sphere.
%
%  Discussion:
%
%    The integral is approximated by the average of the vertex values,
%    multiplied by the area.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 April 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real V1(3), V2(3), V3(3), the XYZ coordinates of
%    the vertices of the triangle.
%
%    Input, function v = F ( x ), evaluates the integrand at X.
%
%    Output, real RESULT, the approximate integral.
%
  area = sphere01_triangle_vertices_to_area ( v1, v2, v3 );

  quad = ( f ( v1 ) + f ( v2 ) + f ( v3 ) ) / 3.0;

  result = quad * area;

  return
end


text:
function area = polyArea(v,varargin)
% area of a flat polygon given by vertices v1, v2, ..., v_n
%
% Input
%  v - @vector3d
%
% Output
%  area - area the polygon v1, v2, ..., v_n
%

nv = length(v);
if nv<3, area = 0; return; end

% formula from http://geomalgorithms.com/a01-_area.html
v1 = v.subSet(1);
N = normalize(mean(cross(v.subSet(2)-v1,v.subSet(3:nv-1)-v1)));

area = dot(N, sum(cross(v.subSet(1:nv-1),v.subSet(2:nv))))./2;

author: mtex-toolbox

repo: mtex

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


text:
function imdisp(I)
% balance the gray level range to display the image
x = (0:255)./255;
grey = [x;x;x]';
minI = min(min(I));
maxI = max(max(I));
I = (I-minI)/(maxI-minI)*255;
image(I);
axis('square','off');
colormap(grey);


text:
function [data,units] = compute_max_absdwing_length(trx,n)

flies = trx.exp2flies{n};
nflies = numel(flies);
data = cell(1,nflies);
for i = 1:nflies,
  fly = flies(i);

  data{i} = max(abs(diff(trx(fly).wing_lengthl_mm,1,2)),abs(diff(trx(fly).wing_lengthr_mm,1,2))) ./ trx(fly).dt;

end
units = parseunits('mm/s');

author: kristinbranson

repo: JAABA


text:
%Finite Element Method 101-2
%National Taiwan University
%2D Elasticity problem

%%clear memory
close all; clc; clear all;
format long;
%% Load Mesh for Lab 10
[nodeCoordinates,elementNodes]=Mesh_Lab10('T3');
% node coordinates are given in mm
NodePerElement=3;
numberNodes=length(nodeCoordinates);
numberElements=length(elementNodes);

%% Import BCs

% Essential BC's
GDof = 2*numberNodes;
prescribedDof = [1,2,3,4];

% Natural BC's
force = zeros(GDof,1);
force(end) = -10; % 10 [kN]

%% Import material and section properties
E = 3E7; % [GPa]
poisson = 0.3; %[-]
thickness = 1; % [mm]

%% Evalute force vector
%force=formForceVectorT3(GDof,naturalBCs,surfaceOrientation,...
%    elementNodes,nodeCoordinates,P,thickness);

%% Construct Stiffness matrix for T3 element
D=E/(1-poisson^2)*[1 poisson 0;poisson 1 0;0 0 (1-poisson)/2];

stiffness=formStiffness2D(GDof,numberElements,...
    elementNodes,numberNodes,nodeCoordinates,D,thickness);

%% solution
displacements=solution(GDof,prescribedDof,stiffness,force);

%% output displacements
outputDisplacements(displacements, numberNodes, GDof);
scaleFactor=1.E5;
drawingMesh(nodeCoordinates+scaleFactor*[displacements(1:2:2*numberNodes) ...
    displacements(2:2:2*numberNodes)],elementNodes,'T3','r--');

% Computes elements stresses
for e=1:numberElements
  numNodePerElement = length(elementNodes(e,:));
  numEDOF = 2*numNodePerElement;
  elementDof=zeros(1,numEDOF);
  for i = 1:numNodePerElement
      elementDof(2*i-1)=2*elementNodes(e,i)-1;
      elementDof(2*i)=2*elementNodes(e,i);
  end

  %  B matrix
  x1 = nodeCoordinates(elementNodes(e,1),1);
  y1 = nodeCoordinates(elementNodes(e,1),2);
  x2 = nodeCoordinates(elementNodes(e,2),1);
  y2 = nodeCoordinates(elementNodes(e,2),2);
  x3 = nodeCoordinates(elementNodes(e,3),1);
  y3 = nodeCoordinates(elementNodes(e,3),2);
  A = 1/2*det([1 x1 y1; 1 x2 y2; 1 x3 y3]);
  B = 1/(2*A).*[y2-y3 0 y3-y1 0 y1-y2 0;
                        0 x3-x2 0 x1-x3 0 x2-x1;
                        x3-x2 y2-y3 x1-x3 y3-y1 x2-x1 y1-y2];

  stress=D*B*displacements(elementDof);
  vonmises=sqrt(0.5*((stress(1)-(stress(2)))^2+(stress(2))^2+(stress(1))^2+6*(stress(3))^2));
  fprintf('\n Stress in element % u \n',e)
  fprintf('Sigma_xx : %0 .6f \n',stress(1))
  fprintf('Sigma_yy : %0 .6f \n',stress(2))
  fprintf('Sigma_xy : %0 .6f \n',stress(3))
  fprintf('Vonmises : %0 .6f \n',vonmises)
end

repo: Aero-matlab


text:
function value = parallelipiped_volume_3d ( x, y, z )

%*****************************************************************************80
%
%% PARALLELIPIPED_VOLUME_3D returns the volume of a parallelipiped in 3D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 November 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X(4), Y(4), Z(4), the coordinates of one corner
%    of the parallelipiped, and its 3 immediate neighbors.
%
%    Output, real PARALLELIPIPED_VOLUME_3D, the volume of
%    the parallelipiped.
%
  value = abs ( ...
    ( z(2) - z(1) ) * ( y(4) * x(3) - y(3) * x(4) ) + ...
    ( z(3) - z(1) ) * ( x(4) * y(2) - x(2) * y(4) ) + ...
    ( z(4) - z(1) ) * ( x(2) * y(3) - x(3) * y(2) ) + ...
    ( z(3) - z(2) ) * ( y(4) * x(1) - y(1) * x(4) ) + ...
    ( z(4) - z(2) ) * ( x(3) * y(1) - x(1) * y(3) ) + ...
    ( z(4) - z(3) ) * ( x(1) * y(2) - x(2) * y(1) ) );

  return
end


text:
function [xi,w]=fifthOrderCrossPolyCubPoints(numDim,algorithm)
%%FIFTHORDERCROSSPOLYCUBPOINTS Generate third order cubature points for
%               integration over the region sum(abs(x))<=1 with weighting
%               function w(x)=1. In two dimensions, this is a square, in
%               three an octohedron, in n a cross polytope (or
%               n-dimensional octahedron).
%
%INPUTS: numDim An integer specifying the dimensionality of the points to
%               be generated.
%     algorithm A value indicating which algorithm should be used.
%               Possible values are
%               0 (The default if omitted or an empty matrix is passed)
%                 Formula G_n 5-2 in [1], pg. 305, 4*numDim^2-4*numDim+1
%                 points, 2<=numDim<=7.
%               1 G_n 5-3 in [1], pg. 305, 2^numDim+2*numDim points,
%                 numDim=3,4,5, with a correction that the first term with
%                 (r;0;0...0) be put through fullSymPerms.
%               2 G_n 5-4 in [1], pg. 306, 2^{numDim+1)-1 points, n<4.
%               3 G_3 5-1 in [1], pg. 306, 13 points, numDim=3.
%
%OUTPUTS: xi A numDim X numCubaturePoints matrix containing the cubature
%              points. (Each "point" is a vector)
%            w A numCubaturePoints X 1 vector of the weights associated
%              with the cubature points.
%
%REFERENCES:
%[1] A.H. Stroud, Approximate Calculation of Multiple Integrals. Cliffs,
%    NJ: Prentice-Hall, Inc., 1971.
%
%November 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release

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

switch(algorithm)
    case 0%G_n 5-2 in [1], pg. 305, 4*numDim^2-4*numDim+1 points,
          %2<=numDim<=7.
        if(numDim<2||numDim>7)
           error('This formula requires 2<=numDim<=7.')
        end
        V=exp(numDim*log(2)-gammaln(numDim+1));
        n=numDim;

        r=sqrt((n+5-sqrt((n+5)*(7-n)))/((n+3)*(n+4)));
        s=sqrt((n+5+sqrt((n+5)*(7-n)))/((n+3)*(n+4)));

        B0=V*(n^2+5*n+10)/((n+1)*(n+2)*(n+5));
        B1=V*(n+3)*(n+4)/(4*(n-1)*(n+1)*(n+2)*(n+5));

        xi=[zeros(numDim,1),fullSymPerms([r;s;zeros(numDim-2,1)])];
        w=[B0;B1*ones(4*numDim^2-4*numDim,1)];
    case 1%G_n 5-3 in [1], pg. 305, 2^numDim+2*numDim points, numDim=3,4,5,
          %with a correction that the first term with (r;0;0...0) must be
          %put through fullSymPerms.
        if(numDim~=3 &&numDim~=4&&numDim~=5)
           error('This formula requires numDim=3, 4, or 5.')
        end
        V=exp(numDim*log(2)-gammaln(numDim+1));
        n=numDim;

        r=sqrt((5*(n+3)*(n+4)+sqrt(5*(n+3)*(n+4)*(n^2+5*n+10)))/((n+3)*(n+4)*(2*n+5)));
        s=sqrt((2*n*(n+3)*(n+4)-2*sqrt(5*(n+3)*(n+4)*(n^2+5*n+10)))/((n-2)*(n+3)*(n+4)*(n^2+4*n+5)));

        B1=10*V/((n+1)*(n+2)*(n+3)*(n+4)*r^4);
        B2=2^(2-n)*V/((n+1)*(n+2)*(n+3)*(n+4)*s^4);

        xi=[fullSymPerms([r;zeros(n-1,1)]),PMCombos(s*ones(numDim,1))];
        w=[B1*ones(2*n,1);B2*ones(2^n,1)];
    case 2%G_n 5-4 in [1], pg. 306, 2^{numDim+1)-1 points, n<4.
        if(numDim>=4)
           error('This formula requires numDim<4')
        end
        V=exp(numDim*log(2)-gammaln(numDim+1));
        n=numDim;

        s=sqrt(2/((n+3)*(n+4)));
        B0=V*(n^2+5*n+10)/((n+1)*(n+2)*(n+5));

        numPoints=2^(n+1)-1;
        xi=zeros(numDim,numPoints);
        w=zeros(numPoints,1);
        w(1)=B0;

        curStart=2;
        vec=s*ones(numDim,1);
        for k=1:numDim
            r=sqrt(2*(k+5)/((n+3)*(n+4)));
            B=V*2^(k-n-1)*5*(n+3)*(n+4)/((n+1)*(n+2)*(k+4)*(k+5));

            vec(k)=r;
            xiCur=PMCombos(vec(k:end));
            num2Add=size(xiCur,2);
            xi(k:end,curStart:(curStart+num2Add-1))=xiCur;
            w(curStart:(curStart+num2Add-1))=B;

            curStart=curStart+num2Add;
        end
    case 3%G_3 5-1 in [1], pg. 306, 13 points, numDim=3.
        if(numDim~=3)
           error('This formula requires numDim=3')
        end
        V=exp(numDim*log(2)-gammaln(numDim+1));

        r=sqrt((4+2*sqrt(2))/21);
        s=sqrt((4-2*sqrt(2))/21);
        B0=V*68/320;
        B1=V*21/320;

        xi=[[0;0;0],PMCombos([r;s;0]),PMCombos([0;r;s]),PMCombos([s;0;r])];
        w=[B0;B1*ones(12,1)];
    otherwise
        error('Unknown algorithm specified');
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:
%% Uninstall
%
% Copyright (c) 2018-present, Mahmoud Afifi
% York University, Canada
% mafifi@eecs.yorku.ca | m.3afifi@gmail.com
%
% This source code is licensed under the license found in the
% LICENSE file in the root directory of this source tree.
% All rights reserved.
%
% Please cite the following work if this program is used:
% Mahmoud Afifi, Brian Price, Scott Cohen, and Michael S. Brown,
% "When color constancy goes wrong: Correcting improperly white-balanced
% images", CVPR 2019.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%

disp('Uninstalling WB for sRGB images...')
current = pwd;
rmpath(fullfile(current,'classes'));
savepath
disp('Done!');

author: mahmoudnafifi


text:
function [Pro,FreeNodec] = transferP1red(elemc,elemf,FreeNode)
%% TRANSFERP2RED
%
%

if ~exist('FreeNode','var'), FreeNode = []; end

%% Data structure
N = max(elemc(:));
if exist('elemf','var')
    [elem,HB] = uniformcoarsenred(elemf);
else
    [elem2dof,edge] = dofP2(elemc);
    NE = size(edge,1);
    HB(:,[1 2 3]) = [(N+1:N+NE)', edge(:,1:2)];
end

%% Prolongation matrix
coarseNode = (1:N)';

repo: ifem


% save the tracking results to txt required by the MOTChallenge
function save_results(seq_name, file_name, bboxes, initialization_threshold, opt)

num_target = max(bboxes.id);
len_traj = zeros(num_target, 1);
for i = 1:num_target
    len_traj(i) = numel(find(bboxes.id == i & (bboxes.state == opt.STATE_TRACKED | bboxes.state == opt.STATE_ACTIVATED)));
end

fid = fopen(file_name, 'w');
num_bbox = numel(bboxes.x);

for i = 1:num_bbox
    % <frame>, <id>, <bb_left>, <bb_top>, <bb_width>, <bb_height>, <conf>, <x>, <y>, <z>
    if len_traj(bboxes.id(i)) > initialization_threshold && (bboxes.state(i) == opt.STATE_TRACKED || bboxes.state(i) == opt.STATE_ACTIVATED)
        fprintf(fid, '%d,%d,%f,%f,%f,%f,%f,%f,%f,%f\n', ...
            bboxes.fr(i), bboxes.id(i), bboxes.x(i), bboxes.y(i), bboxes.w(i), bboxes.h(i), -1, -1, -1, -1);
    end
end

fclose(fid);


text:
function [trans, rot] = mrDoAlignVol(inpts, volpts, scaleFac, inpSize, ...
	sagSize, volume, numSlices, retwin, volwin, obwin);
%
%[trans, rot] = mrDoAlignVol(inpts, volpts, scaleFac, inpSize, ...
%	sagSize, volume, numSlices, retwin, volwin, obwin);
%
%	returns alignment matrix given inpts and volpts as corresponding points
%	rot rotates inpts into volpts coordinate frame.
%	scaleFac is a vector containing scalings of the x,y,and z axes
%		such that inpts*scalefac is at the same scale as volpts

global volslimin1 volslimax1;

nuinpts = inpts ./ (ones(length(inpts),1)*scaleFac(1,:));
nuvolpts = volpts ./ (ones(length(volpts),1)*scaleFac(2,:));
nuvolpts = nuvolpts - (mean(nuvolpts)'*ones(1,length(nuvolpts)))';
nuinpts = nuinpts - (mean(nuinpts)'*ones(1,length(nuinpts)))';

H = zeros(3,3);
for i = 1:length(nuvolpts)
	H = H + (nuinpts(i,:)')*(nuvolpts(i,:));
end
[U,S,V] = svd(H);

mirrorFixer = [ 1 0 0;0 1 0; 0 0 det(U*V);];
rot = [V*mirrorFixer*(U')];

if det(rot) == -1
	disp('Warning: rotation matrix has -1 determinant');
end

alinpts = (rot*(inpts'./(ones(length(inpts),1)*scaleFac(1,:))'))';
nuvolpts = volpts ./ (ones(length(volpts),1)*scaleFac(2,:));
trans = mean(nuvolpts) - mean(alinpts);

% Check that vol = rot*inp + trans;
%alinpts = alinpts+(trans'*ones(1,length(alinpts)))';
%figure(retwin);
%hold on
%plot(inpts(:,1),inpts(:,2),'r-');
%hold off
%figure(volwin);
%hold on
%plot(alinpts(:,1),alinpts(:,2),'r-');
%hold off

% Check it out by displaying anatomy rotated to inplanes first plane

img = mrCheckAlignVol(rot,trans,scaleFac,inpSize,inpts(size(inpts,1),3), ...
			volume,sagSize,numSlices,obwin);

author: vistalab

repo: vistasoft


text:
function [ w, cvx_optval ] = fdla( A ) %#ok

% Computes the fastest distributed linear averaging (FDLA) edge weights
%
% [W,S] = FDLA(A) gives a vector of the fastest distributed linear averaging
% edge weights for a graph described by the incidence matrix A (n x m).
% Here n is the number of nodes and m is the number of edges in the graph;
% each column of A has exactly one +1 and one -1.
%
% The FDLA edge weights are given by the SDP:
%
%   minimize    s
%   subject to  -s*I <= I - L - (1/n)11' <= s*I
%
% where the variables are edge weights w in R^m and s in R.
% Here L is the weighted Laplacian defined by L = A*diag(w)*A'.
% The optimal value is s, and is returned in the second output.
%
% For more details see the references:
% "Fast linear iterations for distributed averaging" by L. Xiao and S. Boyd
% "Convex Optimization of Graph Laplacian Eigenvalues" by S. Boyd
%
% Written for CVX by Almir Mutapcic 08/29/06

[n,m] = size(A); %#ok
I = eye(n,n);
J = I - (1/n) * ones(n,n);
cvx_begin sdp
    variable w(m,1)   % edge weights
    variable s        % epigraph variable
    variable L(n,n) symmetric
    minimize( s )
    subject to
        L == A * diag(w) * A'; %#ok
        -s * I <= J - L <= s * I; %#ok
cvx_end

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function bbs = acfDetect( I, detector, fileName )
% Run aggregate channel features object detector on given image(s).
%
% The input 'I' can either be a single image (or filename) or a cell array
% of images (or filenames). In the first case, the return is a set of bbs
% where each row has the format [x y w h score] and score is the confidence
% of detection. If the input is a cell array, the output is a cell array
% where each element is a set of bbs in the form above (in this case a
% parfor loop is used to speed execution). If 'fileName' is specified, the
% bbs are saved to a comma separated text file and the output is set to
% bbs=1. If saving detections for multiple images the output is stored in
% the format [imgId x y w h score] and imgId is a one-indexed image id.
%
% A cell of detectors trained with the same channels can be specified,
% detected bbs from each detector are concatenated. If using multiple
% detectors and opts.pNms.separate=1 then each bb has a sixth element
% bbType=j, where j is the j-th detector, see bbNms.m for details.
%
% USAGE
%  bbs = acfDetect( I, detector, [fileName] )
%
% INPUTS
%  I          - input image(s) of filename(s) of input image(s)
%  detector   - detector(s) trained via acfTrain
%  fileName   - [] target filename (if specified return is 1)
%
% OUTPUTS
%  bbs        - [nx5] array of bounding boxes or cell array of bbs
%
% EXAMPLE
%
% See also acfTrain, acfModify, bbGt>loadAll, bbNms
%
% Piotr's Image&Video Toolbox      Version 3.20
% Copyright 2013 Piotr Dollar & Ron Appel.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Simplified BSD License [see external/bsd.txt]

% run detector on every image
if(nargin<3), fileName=''; end; multiple=iscell(I);
if(~isempty(fileName) && exist(fileName,'file')), bbs=1; return; end
if(~multiple), bbs=acfDetectImg(I,detector); else
  n=length(I); bbs=cell(n,1);
  parfor i=1:n, bbs{i}=acfDetectImg(I{i},detector); end
end

% write results to disk if fileName specified
if(isempty(fileName)), return; end
d=fileparts(fileName); if(~isempty(d)&&~exist(d,'dir')), mkdir(d); end
if( multiple ) % add image index to each bb and flatten result
  for i=1:n, bbs{i}=[ones(size(bbs{i},1),1)*i bbs{i}]; end
  bbs=cell2mat(bbs);
end
dlmwrite(fileName,bbs); bbs=1;

end

function bbs = acfDetectImg( I, detector )
% Run trained sliding-window object detector on given image.
Ds=detector; if(~iscell(Ds)), Ds={Ds}; end; nDs=length(Ds);
opts=Ds{1}.opts; pPyramid=opts.pPyramid; pNms=opts.pNms;
imreadf=opts.imreadf; imreadp=opts.imreadp;
shrink=pPyramid.pChns.shrink; pad=pPyramid.pad;
separate=nDs>1 && isfield(pNms,'separate') && pNms.separate;
% perform actual computations
if(all(ischar(I))), I=feval(imreadf,I,imreadp{:}); end
P = chnsPyramid(I,pPyramid); bbs = cell(P.nScales,nDs);
for i=1:P.nScales
  for j=1:nDs, opts=Ds{j}.opts;
    modelDsPad=opts.modelDsPad; modelDs=opts.modelDs;
    bb = acfDetect1(P.data{i},Ds{j}.clf,shrink,...
      modelDsPad(1),modelDsPad(2),opts.stride,opts.cascThr);
    shift=(modelDsPad-modelDs)/2-pad;
    bb(:,1)=(bb(:,1)+shift(2))/P.scaleshw(i,2);
    bb(:,2)=(bb(:,2)+shift(1))/P.scaleshw(i,1);
    bb(:,3)=modelDs(2)/P.scales(i);
    bb(:,4)=modelDs(1)/P.scales(i);
    if(separate), bb(:,6)=j; end; bbs{i,j}=bb;
  end
end; bbs=cat(1,bbs{:});
if(~isempty(pNms)), bbs=bbNms(bbs,pNms); end
end

author: zouchuhang

repo: LayoutNet


text:
%compute dtailinflang
function [data,units]=compute_dtailinflang(trx,n)

larvae=trx.exp2flies{n};
numlarvae=numel(larvae);
%absdtailheadang=cell(1,numlarvae);
dtailinflang=cell(1,numlarvae);
for i=1:numlarvae
    larva=larvae(i);
    dtailinflang{1,i}=modrange(diff(trx(larva).tailinflang),-pi,pi)./trx(larva).dt;
    % KB: this angle distance didn't make sense to me, changed to one that
    % made sense to me
    % dtailinflang{1,i}=(mod1(trx(larva).tailinflang(2:end)-trx(larva).tailinflang(1:end-1),pi)-pi)./trx(larva).dt;
end

units=parseunits('rad/s');
data=dtailinflang;


text:
function res = vl_ffdnet_concise(net, x)

global sigmas;
n = numel(net.layers);
res = struct('x', cell(1,n+1));
res(1).x = x ;
cudnn = {'CuDNN'} ;
%cudnn = {'NoCuDNN'} ;

for i=1:n
    l = net.layers{i} ;
    switch l.type
        case 'conv'
            res(i+1).x = vl_nnconv(res(i).x, l.weights{1}, l.weights{2}, ...
                'pad', l.pad, ...
                'stride', l.stride, ...
                'dilate', l.dilate, ...
                l.opts{:}, ...
                cudnn{:}) ;

        case 'concat'
            if size(sigmas,1)~=size(res(i).x,1)
                sigmaMap   = bsxfun(@times,ones(size(res(i).x,1),size(res(i).x,2),1,size(res(i).x,4),'single'),permute(sigmas,[3 4 1 2]));
                res(i+1).x = cat(3,res(i).x,sigmaMap);
            else
                res(i+1).x = cat(3,res(i).x,sigmaMap);
            end

        case 'SubP'
            res(i+1).x = vl_nnSubP(res(i).x, [],'scale',l.scale);

        case 'relu'
            res(i+1).x = max(res(i).x,0) ;
    end
        res(i).x = [] ;
end

author: cszn

repo: DnCNN

save_path: github-repos/MATLAB/cszn-DnCNN


text:
function d = p01_fd ( p )

%*****************************************************************************80
%
%% P01_FD is a signed distance function for problem 1.
%
%  Licensing:
%
%    (C) 2004 Per-Olof Persson.
%    See COPYRIGHT.TXT for details.
%
%  Modified:
%
%    06 February 2006
%
%  Parameters:
%
%    Input, real P, one or more points.
%
%    Output, real D, the signed distance of each point to the boundary of the region.
%
  d = sqrt ( sum ( p.^2, 2 ) ) - 1.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function scale = p24_scale ( neqn )

%*****************************************************************************80
%
%% P24_SCALE returns scale factors for problem p24.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Output, real SCALE(NEQN), the scaling factors.
%
  scale(1:neqn,1) = [ 1.0; 1.0 ];

  return
end