text:
function [] = showChannelTrend(EEG, EEGChannels, displayChannels, seconds)
%% Displays the trend line for EEG over an interval
%
% Parameters
%    EEG              EEGLAB EEG structure
%    EEGChannels      vector numbers of channels corresponding to EEG
%    displayChannels  vector containing list of channels to visualize
%    seconds          interval (in seconds) to analyze (empty = all)
%
% This function shows a mumber of visualizations both in time and
% frequency space
%

%% Set up the data to be visualized
EEGOrig = EEG;
t = 0:(size(EEGOrig.data, 2) - 1);
t = t/EEGOrig.srate;
if ~isempty(seconds)
    tIndex = (t >= seconds(1) & t <= seconds(2));
    t = t(tIndex);
    EEGOrig.data = EEGOrig.data(:, tIndex);
    EEGOrig.pnts = sum(tIndex);
    EEGOrig.times = EEGOrig.times(tIndex);
end

%% Display channels must be EEG channels
displayChannels = intersect(EEGChannels, displayChannels);

%% Calculate the correlations and the models

myData = EEGOrig.data;
myPval = zeros(size(myData, 1), 2);
myCorr = zeros(size(myData, 1), 1);
for k = 1:size(myData, 1)
  myPval(k, :) = polyfit(t, myData(k, :), 1);
  myCorr(k) = corr(t', (myData(k, :))');
end

%% Display the correlation values and the polynomial coefficients
cTitle = 'Summary of trends for EEG channels';
figure('Color', [1, 1, 1], 'Name', cTitle)
hold on
line([-1, 1], [0, 0], 'Color', [.85, .85, .85], 'LineWidth', 3);
plot(myCorr(EEGChannels), myPval(EEGChannels, 1), 'xk', 'MarkerSize', 12, 'LineWidth', 2)
xlabel('Channel-trend correlation')
ylabel('Slope of trend')
title(cTitle)
box on
hold off

%%
for k = 1:length(displayChannels)
    chan = displayChannels(k);
    myTitle = ['Channel ' num2str(chan) ' (corr=' num2str(myCorr(chan)) ')'];
    figure('Color', [1, 1, 1], 'Name', myTitle)
    plot(t, myData(chan, :), '-k')
    t1 = get(gca, 'XLim');
    hold on
    line(t1, myPval(chan, 1)*t1 + myPval(chan, 2), 'Color', [1, 0, 0], 'LineWidth', 2);
    xlabel('Seconds')
    ylabel('Voltage (uV)')
    hold off
    box on
end

%% Display the correlation values and the polynomial coefficients
%% Parameters that must be preset
chans = displayChannels;
EEGSmall = EEGOrig;
EEGSmall.nbchan = length(chans);
EEGSmall.data = EEGOrig.data(chans, :);
params.detrendType = 'high pass';
params.detrendCutoff = 1;
[EEGLine, lineNoiseOut] = cleanLineNoise(EEGSmall, params);
EEGLineFilt = removeTrend(EEGLine, params);
EEGFilt = removeTrend(EEGSmall, params);
EEGFiltLine = cleanLineNoise(EEGFilt, params);
%% Remove from model
t = 0:(size(EEGSmall.data, 2) - 1);
t = t/EEGSmall.srate;
EEGModel = EEGSmall;
for k = 1:length(chans)
  ch = chans(k);
  myModel = polyval(myPval(ch, :), t);
  EEGModel.data(k, :) = EEGModel.data(k, :) - myModel;
end
%%
EEGModelLine = cleanLineNoise(EEGModel, params);
EEGModelLineFilt = removeTrend(EEGModelLine, params);

%% Compute the filtering for models and EEG
SOrig = cell(length(chans), 1);
fOrig = cell(length(chans), 1);
SFilt = cell(length(chans), 1);
fFilt = cell(length(chans), 1);
SFiltLine = cell(length(chans), 1);
fFiltLine = cell(length(chans), 1);
SLine = cell(length(chans), 1);
fLine = cell(length(chans), 1);
SLineFilt = cell(length(chans), 1);
fLineFilt = cell(length(chans), 1);
SModel = cell(length(chans), 1);
fModel = cell(length(chans), 1);
SModelLine = cell(length(chans), 1);
fModelLine = cell(length(chans), 1);
SModelLineFilt = cell(length(chans), 1);
fModelLineFilt = cell(length(chans), 1);
winSize = lineNoiseOut.taperWindowSize;
sParams = lineNoiseOut;
for k = 1:length(chans)
   [SOrig{k},fOrig{k}] = mtspectrumsegc(EEGSmall.data(k, :), winSize, sParams);
   [SLine{k},fLine{k}] = mtspectrumsegc(EEGLine.data(k, :), winSize, sParams);
   [SLineFilt{k},fLineFilt{k}] = mtspectrumsegc(EEGLineFilt.data(k, :), winSize, sParams);
   [SFilt{k},fFilt{k}] = mtspectrumsegc(EEGFilt.data(k, :), winSize, sParams);
   [SFiltLine{k},fFiltLine{k}] = mtspectrumsegc(EEGFiltLine.data(k, :), winSize, sParams);
   [SModel{k},fModel{k}] = mtspectrumsegc(EEGModel.data(k, :), winSize, sParams);
   [SModelLine{k},fModelLine{k}] = mtspectrumsegc(EEGModelLine.data(k, :), winSize, sParams);
   [SModelLineFilt{k},fModelLineFilt{k}] = mtspectrumsegc(EEGModelLineFilt.data(k, :), winSize, sParams);
end
%% Figures with the entire spectral range and all filtering variations
thisName = EEG.setname;
colors = jet(10);
for k = 1:length(chans)
    tString = [thisName ': Channel ' num2str(chans(k)) ' (all versions)'];
    figure('Name', tString, 'Color', [1, 1, 1])
    hold on
    plot(fOrig{k}, 10*log10(SOrig{k}), 'Color', [0.85, 0.85, 0.85], 'LineWidth', 3', 'LineStyle', '-')
    plot(fLine{k}, 10*log10(SLine{k}), 'Color', [0, 0, 0])
    plot(fLineFilt{k}, 10*log10(SLineFilt{k}), 'Color', colors(1, :))
    plot(fFilt{k}, 10*log10(SFilt{k}), 'Color', colors(2, :))
    plot(fFiltLine{k}, 10*log10(SFiltLine{k}), 'Color', colors(4, :))
    plot(fModel{k}, 10*log10(SModel{k}), 'Color', colors(6, :));
    plot(fModelLine{k}, 10*log10(SModelLine{k}), 'Color', colors(7, :))
    plot(fModelLineFilt{k}, 10*log10(SModelLineFilt{k}), 'Color', colors(10, :)')
    legend('Orig', 'Line', 'L-Filt', 'Filt', 'Filt-L', 'Model', 'Model-L', 'Model-LF');
    hold off
    xlabel('Hz')
    ylabel('Power')
    title(tString, 'Interpreter', 'none')
    box on
end
%% Figures with just interaction of filtering and line noise -- limited range
for k = 1:length(chans)
    tString = [thisName ': Channel ' num2str(chans(k)) ' (limited spectrum)'];
    figure('Name', tString, 'Color', [1, 1, 1])
    hold on
    plot(fOrig{k}, 10*log10(SOrig{k}), 'Color', [0.7, 0.7, 0.7], 'LineWidth', 3, 'LineStyle', '-')
    plot(fFiltLine{k}, 10*log10(SFiltLine{k}), 'Color', [0.85, 0.85, 0.85], 'LineWidth', 6,  'LineStyle', '-')
    plot(fLineFilt{k}, 10*log10(SLineFilt{k}), 'Color', [0, 0, .8])
    plot(fModelLineFilt{k}, 10*log10(SModelLineFilt{k}), 'Color', [1, 0.2, 0.2])
    legend('Orig', 'Filt-L', 'L-Filt',  'Model-LF');
    hold off
    xlabel('Hz')
    ylabel('Power')
    set(gca, 'XLim', [0, 110], 'XLimMode', 'manual')
    title(tString, 'Interpreter', 'none')
    box on
end

author: VisLab

repo: EEG-Clean-Tools


text:
function [y,coef,window,Cx,Ff] = lanczosfilter(x,varargin)
%LANCZOSFILTER   Filters a time series via Lanczos method (cosine filter).
%   [Y,coef,window,Cx,Ff] = LANCZOSFILTER(X,dT,Cf,M,pass) Filters the time
%   series via the Lanczos filter in the frequency space (FFT), where
%
%   INPUTS:
%      X    - Time series
%      dT   - Sampling interval       (default: 1)
%      Cf   - Cut-off frequency       (default: half Nyquist)
%      M    - Number of coefficients  (default: 100)
%      pass - Low or high-pass filter (default: 'low')
%
%   OUTPUTS:
%      Y      - Filtered time series
%      coef   - Coefficients of the time window (cosine)
%      window - Frequency window (aprox. ones for Ff lower(greater) than Fc
%               if low(high)-pass filter and ceros otherwise)
%      Cx     - Complex Fourier Transform of X for Ff frequencies
%      Ff     - Fourier frequencies, from 0 to the Nyquist frequency.
%
%   The program removes from the time series the frequencies greater than
%   the cut off frequency if "pass" is 'low', i.e., low-pass filter .
%   Otherwise, if pass is 'high', frequencies from zero to Cf are removed,
%   i.e., a high-pass filter. Units of the cut-off frequency, [Cf], must be
%   [dT]^{-1}.
%
%   In consequence, when used as a low-pass the time series is smoothed
%   like a cosine filter in time space with M coefficients where greater is
%   better (see the reference).
%
%   If any option is empty, defaults are used.
%
%   Note: NaN's elements are replaced by the mean of the time series and
%         ignored. If you have a better idea, just let me know.
%
%   Reference:
%   Emery, W. J. and R. E. Thomson. "Data Analysis Methods in Physical
%     Oceanography". Elsevier, 2d ed., 2004. On pages 533-539.
%
%   Example:
%      dT = 30; % min
%      N = 7*24*60/dT; t = (0:N-1)*dT; % data for 7 days
%      pnoise = 0.30;
%      T1 = 12.4*60; T2 = 24*60; T3 = 15*24*60; Tc = 10*60; % min
%      xn = 5 + 3*cos(2*pi*t/T1) + 2*cos(2*pi*t/T2) + 1*cos(2*pi*t/T3);
%      xn = xn + pnoise*max(xn-mean(xn))*(0.5 - rand(size(xn)));
%      [xs,c,h,Cx,f] = lanczosfilter(xn,dT,1/Tc,[],'low');
%      subplot(211), plot(t,xn,t,xs), legend('noisy','smooth'), axis tight
%      subplot(212), plot(f,h,f,abs(Cx)/max(abs(Cx)),...
%         [1 1]/Tc,[min(h) max(h)],'-.',...
%         [1/T1 1/T2 1/T3],([1/T1 1/T2 1/T3]<=1/Tc),'o'), axis tight
%
%   See also FILTER, FFT, IFFT

%   Written by
%   Lic. on Physics Carlos Adrián Vargas Aguilera
%   Physical Oceanography MS candidate
%   UNIVERSIDAD DE GUADALAJARA
%   Mexico, 2004
%
%   nubeobscura@hotmail.com

% Check arguments:
if nargin<1 || nargin>5
 error('Lanczosfilter:ArgumentNumber','Incorrect number of arguments.')
elseif ~isvector(x) || ~isreal(x)
 error('Lanczosfilter:ArgumentType','Incorrect time series.')
end
if nargin<2 || isempty(varargin{1})
  dT = 1;
elseif ~(numel(varargin{1})==1) || ~isreal(varargin{1})
 error('Lanczosfilter:ArgumentType','Incorrect time interval.')
else
 dT = varargin{1};
end
Nf = 1/(2*dT); % Nyquist frequency
if nargin<3 || isempty(varargin{2})
 Cf = Nf/2;
elseif ~(numel(varargin{2})==1) || ~isreal(varargin{2}) || varargin{2}<=0 || varargin{2}>Nf
 error('Lanczosfilter:ArgumentType','Incorrect cut-off frequency.')
else
 Cf = varargin{2};
end
if nargin<4 || isempty(varargin{3})
 M = 100;
elseif ~(numel(varargin{3})==1) || ~isreal(varargin{3}) || (varargin{3}==round(varargin{3}))
 error('Lanczosfilter:ArgumentType','Incorrect Number of coeffients.')
else
 M = varargin{3};
end
if nargin<5 || isempty(varargin{4})
 LoH = 'l';
elseif ~ischar(varargin{4}) || isempty(strfind('lh',lower(varargin{4}(1))))
 error('Lanczosfilter:ArgumentType','Incorrect filter pass type.')
else
 LoH = varargin{4};
end
if strcmpi(LoH(1),'h')
 LoH = 2;
else
 LoH = 1;
end

% Normalize the cut off frequency with the Nyquist frequency:
Cf = Cf/Nf;

% Lanczos cosine coeficients:
coef = lanczos_filter_coef(Cf,M); coef = coef(:,LoH);

% Filter in frequency space:
[window,Ff] = spectral_window(coef,length(x)); Ff = Ff*Nf;

% Replace NaN's with the mean (ideas?):
inan = isnan(x);
xmean = mean(x(~inan));
x(inan) = xmean;

% Filtering:
[y,Cx] = spectral_filtering(x,window);

function coef = lanczos_filter_coef(Cf,M)
% Positive coeficients of Lanczos [low high]-pass.
hkcs = lowpass_cosine_filter_coef(Cf,M);
sigma = [1 sin(pi*(1:M)/M)./(pi*(1:M)/M)];
hkB = hkcs.*sigma;
hkA = -hkB; hkA(1) = hkA(1)+1;
coef = [hkB(:) hkA(:)];

function coef = lowpass_cosine_filter_coef(Cf,M)
% Positive coeficients of cosine filter low-pass.
coef = Cf*[1 sin(pi*(1:M)*Cf)./(pi*(1:M)*Cf)];

function [window,Ff] = spectral_window(coef,N)
% Window of cosine filter in frequency space.
Ff = 0:2/N:1; window = zeros(length(Ff),1);
for i = 1:length(Ff)
 window(i) = coef(1) + 2*sum(coef(2:end).*cos((1:length(coef)-1)'*pi*Ff(i)));
end

function [y,Cx] = spectral_filtering(x,window)
% Filtering in frequency space is multiplication, (convolution in time
% space).
Nx  = length(x);
Cx  = fft(x(:)); Cx = Cx(1:floor(Nx/2)+1);
CxH = Cx.*window(:);
CxH(length(CxH)+1:Nx) = conj(CxH(Nx-length(CxH)+1:-1:2));
y = real(ifft(CxH));

% Carlos Adrián Vargas Aguilera. nubeobscura@hotmail.com

author: Sable

repo: mcbench-benchmarks


text:
%P2.9a
%T1[x(n)]=2^x(n)
n=[0:200];x=cos(0.2*pi*n)+0.5*cos(0.6*pi*n);
t1=2.^x
[x1,n1]=sigshift(t1,n,50); %x(n-k) -> k=50
[x2,n2]=sigfold(t1,n); %x(-n)
[x3,n3]=sigfold(x1,n1); %x(k-n)
[rxx,nrxx]=conv_m(t1,n,x2,n2);
[rxk,nrkx]=conv_m(x1,n1,x3,n3);
[y,m]=sigadd(rxx,nrxx,0.1*rxk,nrkx);
figure(1);
subplot(3,1,1); stem(m,y,'k');title('Autocorrelation ryy');xlabel('m');

%T2[x(n)]=3x(n)+4
n=[0:200];x=cos(0.2*pi*n)+0.5*cos(0.6*pi*n);
t2=3*x+4
[x1,n1]=sigshift(t2,n,50); %x(n-k) -> k=50
[x2,n2]=sigfold(t2,n); %x(-n)
[x3,n3]=sigfold(x1,n1); %x(k-n)
[rxx,nrxx]=conv_m(t2,n,x2,n2);
[rxk,nrkx]=conv_m(x1,n1,x3,n3);
[y,m]=sigadd(rxx,nrxx,0.1*rxk,nrkx);
subplot(3,1,2); stem(m,y,'k');title('Autocorrelation ryy');xlabel('m');

%T3[x(n)]=x(n)+2x(n-1)-x(n-2)
n=[0:200];x=cos(0.2*pi*n)+0.5*cos(0.6*pi*n);
[xa,nxa]=sigshift(x,n,1);
[xb,nxb]=sigshift(x,n,2);
[xa,nxa]=sigadd(2*xa,nxa,-1*xb,nxb);
[t3,n]=sigadd(x,n,xa,nxa)
[x1,n1]=sigshift(t3,n,50); %x(n-k) -> k=50
[x2,n2]=sigfold(t3,n); %x(-n)
[x3,n3]=sigfold(x1,n1); %x(k-n)
[rxx,nrxx]=conv_m(t3,n,x2,n2);
[rxk,nrkx]=conv_m(x1,n1,x3,n3);
[y,m]=sigadd(rxx,nrxx,0.1*rxk,nrkx);
subplot(3,1,3); stem(m,y,'k');title('Autocorrelation ryy');xlabel('m');

%P2.9b
%T1[x(n)]=2^x(n)
%x(n)=random sequence
nrand=[0:100];xrand=rand(size(nrand));
t1=2.^xrand
[x1,n1]=sigshift(t1,n,50); %x(n-k) k=50
[x2,n2]=sigfold(t1,n); %x(-n)
[x3,n3]=sigfold(x1,n1); %x(k-n)
[rxx,nrxx]=conv_m(t1,n,x2,n2);
[rxk,nrkx]=conv_m(x1,n1,x3,n3);
[y,m]=sigadd(rxx,nrxx,0.1*rxk,nrkx);
figure(2);
subplot(3,1,1); stem(m,y,'k');title('Autocorrelation ryy');xlabel('m');

%T2[x(n)]=3x(n)+4
%x(n)=random sequence
nrand=[0:100];xrand=rand(size(nrand));
t2=3*xrand+4
[x1,n1]=sigshift(t2,n,50); %x(n-k) -> k=50
[x2,n2]=sigfold(t2,n); %x(-n)
[x3,n3]=sigfold(x1,n1); %x(k-n)
[rxx,nrxx]=conv_m(t2,n,x2,n2);
[rxk,nrkx]=conv_m(x1,n1,x3,n3);
[y,m]=sigadd(rxx,nrxx,0.1*rxk,nrkx);
subplot(3,1,2); stem(m,y,'k');title('Autocorrelation ryy');xlabel('m');

%T3[x(n)]=x(n)+2x(n-1)-x(n-2)
%x(n)=Random Sequence
nrand=[0:100];xrand=rand(size(nrand));
[xa,nxa]=sigshift(xrand,nrand,1);
[xb,nxb]=sigshift(xrand,nrand,2);
[xa,nxa]=sigadd(2*xa,nxa,-1*xb,nxb);
[t3,n]=sigadd(xrand,nrand,xa,nxa);
[x1,n1]=sigshift(t3,n,50); %x(n-k) -> k=50
[x2,n2]=sigfold(t3,n); %x(-n)
[x3,n3]=sigfold(x1,n1); %x(k-n)
[rxx,nrxx]=conv_m(t3,n,x2,n2);
[rxk,nrkx]=conv_m(x1,n1,x3,n3);
[y,m]=sigadd(rxx,nrxx,0.1*rxk,nrkx);
subplot(3,1,3); stem(m,y,'k');title('Autocorrelation ryy');xlabel('m');

%P2.9c
%T1[x(n)]=2^x(n)
%x(n)=Gaussian random sequence
nrand=[0:100];xrand=randn(size(nrand))
t1=2.^xrand
[x1,n1]=sigshift(t1,n,50); %x(n-k) k=50
[x2,n2]=sigfold(t1,n); %x(-n)
[x3,n3]=sigfold(x1,n1); %x(k-n)
[rxx,nrxx]=conv_m(t1,n,x2,n2);
[rxk,nrkx]=conv_m(x1,n1,x3,n3);
[y,m]=sigadd(rxx,nrxx,0.1*rxk,nrkx);
figure(3);
subplot(3,1,1); stem(m,y,'k');title('Autocorrelation ryy');xlabel('m');

%T2[x(n)]=3x(n)+4
%x(n)=Gaussian random sequence
nrand=[0:100];xrand=randn(size(nrand))
t2=3*xrand+4
[x1,n1]=sigshift(t2,n,50); %x(n-k) -> k=50
[x2,n2]=sigfold(t2,n); %x(-n)
[x3,n3]=sigfold(x1,n1); %x(k-n)
[rxx,nrxx]=conv_m(t2,n,x2,n2);
[rxk,nrkx]=conv_m(x1,n1,x3,n3);
[y,m]=sigadd(rxx,nrxx,0.1*rxk,nrkx);
subplot(3,1,2); stem(m,y,'k');title('Autocorrelation ryy');xlabel('m');

%T3[x(n)]=x(n)+2x(n-1)-x(n-2)
%x(n)=Gaussian Random Sequence
nrand=[0:100];xrand=randn(size(nrand));
[xa,nxa]=sigshift(xrand,nrand,1);
[xb,nxb]=sigshift(xrand,nrand,2);
[xa,nxa]=sigadd(2*xa,nxa,-1*xb,nxb);
[t3,n]=sigadd(xrand,nrand,xa,nxa);
[x1,n1]=sigshift(t3,n,50); %x(n-k) -> k=50
[x2,n2]=sigfold(t3,n); %x(-n)
[x3,n3]=sigfold(x1,n1); %x(k-n)
[rxx,nrxx]=conv_m(t3,n,x2,n2);
[rxk,nrkx]=conv_m(x1,n1,x3,n3);
[y,m]=sigadd(rxx,nrxx,0.1*rxk,nrkx);
subplot(3,1,3); stem(m,y,'k');title('Autocorrelation ryy');xlabel('m');


text:
function a = ge( x, y )

%Disciplined convex programming information for GE (>=):
%   The right-hand side of a less-than constraint must be convex. The
%   left-hand side must be concave. Of course, real constant and affine
%   expressions are both convex and concave and can be used on either
%   side as well.
%
%Disciplined geometric programming information for GE (>=):
%   The right-hand side of a less-than constraint must be log-convex---
%   including positive constants, monomials, posynomials, generalized
%   posynomials, and products thereof. The left-hand side must be log-
%   concave---including positive constants, monomials, reciprocals of
%   log-convex expressions, and products thereof.
%
%Note that CVX does not distinguish between strict greater-than (>) and
%greater-than-or-equal (<=) constraints; they are treated identically.
%Feasible interior-point solvers tend to return points which satisfy
%strict inequality, but not all solvers do.

evalin( 'caller', 'cvx_verify' );
b = cvx_pushcnstr( x, y, '>=' );
if nargout, a = b; end

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: yu-jiang

repo: radpbook

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


text:
function [jsonStr] = jsonWrite(inputVar,saveStr,varargin)
	% Writes an input variable out as a JSON string to a file or sends back to parent function.
	% Biafra Ahanonu
	% started: 2020.06.09 [12:19:01]
	% inputs
		%
	% outputs
		%

	% changelog
		% 2021.08.08 [19:30:20] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package.
	% TODO
		%

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

	%========================
	% DESCRIPTION
	% options.exampleOption = '';
	% get options
	% options = getOptions(options,varargin);
	% display(options)
	% unpack options into current workspace
	% fn=fieldnames(options);
	% for i=1:length(fn)
	% 	eval([fn{i} '=options.' fn{i} ';']);
	% end
	%========================

	try
		if nargin==1
			saveStr = '';
		end

		% Encode the input variable as a JSON string
		jsonStr = jsonencode(inputVar);

		if ~isempty(saveStr)
			fprintf('Saving to %s.\n',saveStr);
			fID = fopen(saveStr,'w');
			fprintf(fID,'%s',jsonStr);
			fclose(fID);
		end

		if nargout==0
			disp(subfxnHumanReadable(jsonStr));
		end

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

function [jsonStr] = subfxnHumanReadable(jsonStr)
	% Make human readable
	jsonStr = strrep(jsonStr, ',', sprintf(',\r'));
	jsonStr = strrep(jsonStr, '[{', sprintf('[\r{\r'));
	jsonStr = strrep(jsonStr, '}]', sprintf('\r}\r]'));
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function value = index_to_level_closed ( dim_num, t, order, level_max )

%*****************************************************************************80
%
%% INDEX_TO_LEVEL_CLOSED determines the level of a point given its index.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 November 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    A Sparse Grid Stochastic Collocation Method for Partial Differential
%    Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2309-2345.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer T(DIM_NUM), the grid indices of a point in a 1D closed rule.
%    0 <= T(I) <= ORDER.
%
%    Input, integer ORDER, the order of the rule.
%
%    Input, integer LEVEL_MAX, the level with respect to which the
%    index applies.
%
%    Output, integer VALUE, the first level on which
%    the point associated with the given index will appear.
%
  value = 0;

  for dim = 1 : dim_num

    s = i4_modp ( t(dim), order );

    if ( s == 0 )

      level = 0;

    else

      level = level_max;

      while ( mod ( s, 2 ) == 0 )
        s = floor ( s / 2 );
        level = level - 1;
      end

    end

    if ( level == 0 )
      level = 1;
    elseif ( level == 1 )
      level = 0;
    end

    value = value + level;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%% RRT*FN for 3 DOF Mobile Robot
%
%

% Do not change anything in rrt.m,
% rrt_star.m and rrt_star_fn.m
%

%

%% Getting started
% Create *main_3dof_mobile_rrt_star_fn.m* file. You can create m-file with any other name, that is
% perfectly okay, this will not affect to the solution of a path/motion
% planning problem. All the sources could be found in examples/ directory
% of the distribution.
%

%% Step 1: Choosing the map
% Firstly, we should define what map to use. It is done by defining map
% structure the sample code is goes after.
%
%   map = struct('name', 'bench_june1.mat', 'start_point', [-14.5 -7.5], 'goal_point', [10 -0.65]);
%

%%
%
% * *name field* defines the file name in maps/ directory
% * *start_point* field defines where the initial point of the problem
% * *goal_point* field define where is the goal point on the map
%

%% Step 2: Setting maximum number of iterations
%
%   max_iter = 20e3;
%
% * *max_iter* variable defines how much iteration should be done to solve
% the path planning problem.

%% Step 3: Setting maximum nodes in the tree
%
%   max_nodes = 5e3;
%
% * *max_nodes* variable defines how many nodes should be stored in the
% tree.

%% Step 4: Do we have to benchmark
%
%   is_benchmark = false;
%
% * *is_benchmark* enables benchmarking. For more details please read the
% sources of rrt.m, rrt_star.m and rrt_star_fn.m

%% Step 5: Setting random seed
%
%   rand_seed = 40;
%
% * *rand_seed* variable is a random seed for random number generator, it
% is used for sampling nodes. It is usually used for benchmarking. We set
% the same map, however the random seed is different for the same map. You
% can use *now* if you don't care about random seed.
%
%   rand_seed = now;
%

%% Step 6: Choosing the class (model) we want
%
%   variant = 'FNSimple3D';
%
% * *variant* defines from what class we should instantiate the object. In
% other words it defines what model we choose for application of RRT.
%
% *FNSimple3D* is a name of a class which contains all necessary methods
% and fields in order to represent simple 3 DOF Mobile Robot model.

%% Step 7: RRT*FN
%
%   rrt_star_fn(map, max_iter, max_nodes, is_benchmark, rand_seed, variant);
%
% Line above runs RRT with given parameters. In addition, *rrt* function
% returns the class object with a certain solution.

%% Sources of *main_3dof_mobile_rrt_star_fn.m*
% Press
% <matlab:edit('examples/main_3dof_mobile_rrt_star_fn.m') here>
% to play with example code.
%
%   % 3 DOF mobile robot example.
%   % by Almaskhan Baimyshev
%   % 08/28/2013
%
%   map = struct('name', 'bench_june1.mat', 'start_point', [-14.5 -7.5], 'goal_point', [10 -0.65]);
%   max_iter = 20e3;
%   max_nodes = 3e3;
%   is_benchmark = false;
%   rand_seed = 40;
%   variant = 'FNSimple3D';
%   result = rrt_star_fn(map, max_iter, max_nodes, is_benchmark, rand_seed, variant);
%

author: olzhas


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

    properties
        tank LaunchVehicleTank

        lwrBnd(1,1) double = 0;
        uprBnd(1,1) double = 0;

        useTf(1,1) = false;
    end

    methods
        function obj = StageTankInitMassOptimVar(tank)
            obj.tank = tank;
            obj.tank.optVar = obj;

            obj.id = rand();
        end

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

            if(obj.useTf)
                x = obj.tank.initialMass;
            end
        end

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

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

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

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

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

        function updateObjWithVarValue(obj, x)
            obj.tank.initialMass = x;
        end

        function nameStrs = getStrNamesOfVars(obj, evtNum, varLocType)
            if(evtNum > 0)
                subStr = sprintf('Event %i',evtNum);
            else
                subStr = varLocType;
            end

            tankName = obj.tank.name;
            if(isempty(tankName))
                tankName = 'Untitled Tank';
            end

            nameStrs = {sprintf('%s Tank "%s" Initial Mass', subStr, tankName)};
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function[]=makefigs_slidetrans
%MAKEFIGS_SLIDETRANS  Makes a sample figure for SLIDETRANS.

figure
M=3000;
t=(0:M-1)';
N=500;
w=hermfun((-N:N)'./(N/4),0);
w=w./sqrt(w'*w);
fs=2*pi*(1:30)./1000;
clear x
x(1:M/2,1)=sin(2*pi*t(1:M/2)./70/3);
x(M/2:M,1)=sin(2*pi*t(M/2:M)./70);
y=slidetrans(x,w,fs,'zeros');
h=wavespecplot(t,x,1./fs,abs(y),1/2);
hlines(70*3/2/pi),hlines(70/2/pi)

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
clear;
%% load file
gcp;                            % start cluster
addpath(genpath('utilities'));
addpath(genpath('deconvolution'));

nam = 'demoMovie.tif';          % insert path to tiff stack here
sframe=1;						% user input: first frame to read (optional, default 1)
num2read=2000;					% user input: how many frames to read   (optional, default until the end)
Y = read_file(nam,sframe,num2read);

%Y = Y - min(Y(:));
if ~isa(Y,'single');    Y = single(Y);  end         % convert to single

[d1,d2,T] = size(Y);                                % dimensions of dataset
d = d1*d2;                                          % total number of pixels

%% Set parameters

K = 40;                                           % number of components to be found
tau = 5;                                          % std of gaussian kernel (half size of neuron)
p = 2;

options = CNMFSetParms(...
    'd1',d1,'d2',d2,...                         % dimensionality of the FOV
    'p',p,...                                   % order of AR dynamics
    'gSig',tau,...                              % half size of neuron
    'merge_thr',0.80,...                        % merging threshold
    'nb',2,...                                  % number of background components
    'min_SNR',3,...                             % minimum SNR threshold
    'space_thresh',0.5,...                      % space correlation threshold
    'cnn_thr',0.2...                            % threshold for CNN classifier
    );
%% Data pre-processing

[P,Y] = preprocess_data(Y,p);
%% fast initialization of spatial components using greedyROI and HALS

[Ain,Cin,bin,fin,center] = initialize_components(Y,K,tau,options,P);  % initialize

% display centers of found components
Cn =  correlation_image(Y); %reshape(P.sn,d1,d2);  %max(Y,[],3); %std(Y,[],3); % image statistic (only for display purposes)
figure;imagesc(Cn);
    axis equal; axis tight; hold all;
    scatter(center(:,2),center(:,1),'mo');
    title('Center of ROIs found from initialization algorithm');
    drawnow;

    %% manually refine components (optional)
refine_components = false;  % flag for manual refinement
if refine_components
    [Ain,Cin,center] = manually_refine_components(Y,Ain,Cin,center,Cn,tau,options);
end

%% update spatial components
Yr = reshape(Y,d,T);
[A,b,Cin] = update_spatial_components(Yr,Cin,fin,[Ain,bin],P,options);

%% update temporal components
P.p = 0;    % set AR temporarily to zero for speed
[C,f,P,S,YrA] = update_temporal_components(Yr,A,b,Cin,fin,P,options);

%% classify components

rval_space = classify_comp_corr(Y,A,C,b,f,options);
ind_corr = rval_space > options.space_thresh;           % components that pass the correlation test
                                        % this test will keep processes

%% further classification with cnn_classifier
try  % matlab 2017b or later is needed
    [ind_cnn,value] = cnn_classifier(A,[d1,d2],'cnn_model',options.cnn_thr);
catch
    ind_cnn = true(size(A,2),1);                        % components that pass the CNN classifier
end

%% event exceptionality

fitness = compute_event_exceptionality(C+YrA,options.N_samples_exc,options.robust_std);
ind_exc = (fitness < options.min_fitness);

%% select components

keep = (ind_corr | ind_cnn) & ind_exc;

%% display kept and discarded components
A_keep = A(:,keep);
C_keep = C(keep,:);
figure;
    subplot(121); montage(extract_patch(A(:,keep),[d1,d2],[30,30]),'DisplayRange',[0,0.15]);
        title('Kept Components');
    subplot(122); montage(extract_patch(A(:,~keep),[d1,d2],[30,30]),'DisplayRange',[0,0.15])
        title('Discarded Components');
%% merge found components
[Am,Cm,K_m,merged_ROIs,Pm,Sm] = merge_components(Yr,A_keep,b,C_keep,f,P,S,options);

%%
display_merging = 1; % flag for displaying merging example
if and(display_merging, ~isempty(merged_ROIs))
    i = 1; %randi(length(merged_ROIs));
    ln = length(merged_ROIs{i});
    figure;
        set(gcf,'Position',[300,300,(ln+2)*300,300]);
        for j = 1:ln
            subplot(1,ln+2,j); imagesc(reshape(A_keep(:,merged_ROIs{i}(j)),d1,d2));
                title(sprintf('Component %i',j),'fontsize',16,'fontweight','bold'); axis equal; axis tight;
        end
        subplot(1,ln+2,ln+1); imagesc(reshape(Am(:,K_m-length(merged_ROIs)+i),d1,d2));
                title('Merged Component','fontsize',16,'fontweight','bold');axis equal; axis tight;
        subplot(1,ln+2,ln+2);
            plot(1:T,(diag(max(C_keep(merged_ROIs{i},:),[],2))\C_keep(merged_ROIs{i},:))');
            hold all; plot(1:T,Cm(K_m-length(merged_ROIs)+i,:)/max(Cm(K_m-length(merged_ROIs)+i,:)),'--k')
            title('Temporal Components','fontsize',16,'fontweight','bold')
        drawnow;
end

%% refine estimates excluding rejected components

Pm.p = p;    % restore AR value
[A2,b2,C2] = update_spatial_components(Yr,Cm,f,[Am,b],Pm,options);
[C2,f2,P2,S2,YrA2] = update_temporal_components(Yr,A2,b2,C2,f,Pm,options);

%% do some plotting

[A_or,C_or,S_or,P_or] = order_ROIs(A2,C2,S2,P2); % order components
K_m = size(C_or,1);
[C_df,~] = extract_DF_F(Yr,A_or,C_or,P_or,options); % extract DF/F values (optional)

figure;
[Coor,json_file] = plot_contours(A_or,Cn,options,1); % contour plot of spatial footprints
%savejson('jmesh',json_file,'filename');        % optional save json file with component coordinates (requires matlab json library)

%% display components

plot_components_GUI(Yr,A_or,C_or,b2,f2,Cn,options);

%% make movie
if (0)
    make_patch_video(A_or,C_or,b2,f2,Yr,Coor,options)
end

author: flatironinstitute

repo: CaImAn-MATLAB


text:
% GET ROTATION MATRIX FROM QUATERNION
function [R] = R_q(q)
% This function defines the equivalent rotation matrix from the
% provided quaternion. (Associated block "Quaternion to DCM")
% INPUT:
% q    - The quaternion rotation
% OUTPUT:
% R_AB - The rotation matrix through A to B
% R_BA - The reverse rotation matrix from B to A

assert(numel(q) == 4,'Input quaternion has wrong dimensions');

R = zeros(3,3);
% Normalise the quaternion
q = unit(q);

% SAME AS THE ROBOTICS TOOL BOX....
R(1,1) = q(1)^2 + q(2)^2 - q(3)^2 - q(4)^2;
R(1,2) = 2*(q(2)*q(3) - q(1)*q(4));
R(1,3) = 2*(q(1)*q(3) + q(2)*q(4));
R(2,1) = 2*(q(2)*q(3) + q(1)*q(4));
R(2,2) = q(1)^2 - q(2)^2 + q(3)^2 - q(4)^2;
R(2,3) = 2*(q(3)*q(4) - q(1)*q(2));
R(3,1) = 2*(q(2)*q(4) - q(1)*q(3));
R(3,2) = 2*(q(1)*q(2) + q(3)*q(4));
R(3,3) = q(1)^2 - q(2)^2 - q(3)^2 + q(4)^2;
end

repo: OpenMAS


text:
function varargout = process_cohere1( varargin )
% PROCESS_COHERE1: Compute the coherence between one signal and all the others, in one file.
%
% USAGE:  OutputFiles = process_cohere1('Run', sProcess, sInputA)
%                       process_cohere1('Test', 1)
%                       process_cohere1('Test', 2)

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

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Coherence 1xN [Deprecated]';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = 'Connectivity';
    sProcess.Index       = 655;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/Connectivity';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'data',     'results',  'matrix'};
    sProcess.OutputTypes = {'timefreq', 'timefreq', 'timefreq'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;

    % === CONNECT INPUT
    sProcess = process_corr1n('DefineConnectOptions', sProcess, 0);
    % === REMOVE EVOKED REPONSE
    sProcess.options.removeevoked.Comment = 'Remove evoked response from each trial';
    sProcess.options.removeevoked.Type    = 'checkbox';
    sProcess.options.removeevoked.Value   = 0;
    sProcess.options.removeevoked.Group   = 'input';
    % === COHERENCE METHOD
    sProcess.options.cohmeasure.Comment = {...
        ['Magnitude-squared Coherence
' ...
        '|C|^2 = |Gxy|^2/(Gxx*Gyy)'], ...
        ['Imaginary Coherence (2019)
' ...
        'IC    = |imag(C)|'], ...
        ['Lagged Coherence (2019)
' ...
        'LC    = |imag(C)|/sqrt(1-real(C)^2)'], ...
        [' Imaginary Coherence (before 2019)
' ...
        ' IC    = imag(C)^2 / (1-real(C)^2) ']; ...
        'mscohere', 'icohere2019','lcohere2019', 'icohere'};
    sProcess.options.cohmeasure.Type    = 'radio_label';
    sProcess.options.cohmeasure.Value   = 'mscohere';
    % === Overlap
    sProcess.options.overlap.Comment = 'Overlap for PSD estimation:' ;
    sProcess.options.overlap.Type    = 'value';
    sProcess.options.overlap.Value   = {50, '%', []};
    % === MAX FREQUENCY RESOLUTION
    sProcess.options.maxfreqres.Comment = 'Maximum frequency resolution:';
    sProcess.options.maxfreqres.Type    = 'value';
    sProcess.options.maxfreqres.Value   = {2,'Hz',2};
    % === HIGHEST FREQUENCY OF INTEREST
    sProcess.options.maxfreq.Comment = 'Highest frequency of interest:';
    sProcess.options.maxfreq.Type    = 'value';
    sProcess.options.maxfreq.Value   = {60,'Hz',2};
    % === OUTPUT MODE
    sProcess.options.outputmode.Comment = {'Save individual results (one file per input file)', 'Concatenate input files before processing (one file)', 'Save average connectivity matrix (one file)'};
    sProcess.options.outputmode.Type    = 'radio';
    sProcess.options.outputmode.Value   = 1;
    sProcess.options.outputmode.Group   = 'output';
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    Comment = sProcess.Comment;
end

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputA) %#ok<DEFNU>
    % Input options
    OPTIONS = process_corr1n('GetConnectOptions', sProcess, sInputA);
    if isempty(OPTIONS)
        OutputFiles = {};
        return
    end

    % Metric options
    OPTIONS.Method = 'cohere';
    OPTIONS.RemoveEvoked  = sProcess.options.removeevoked.Value;
    OPTIONS.MaxFreqRes    = sProcess.options.maxfreqres.Value{1};
    OPTIONS.MaxFreq       = sProcess.options.maxfreq.Value{1};
    OPTIONS.CohOverlap    = 0.50;  % First pre-define the overlap
    OPTIONS.pThresh       = 0.05;
    OPTIONS.CohMeasure    = sProcess.options.cohmeasure.Value;

    % Change the overlap if it is specified
    if isfield(sProcess.options, 'overlap') && isfield(sProcess.options.overlap, 'Value') && ...
       iscell(sProcess.options.overlap.Value) && ~isempty(sProcess.options.overlap.Value) && ~isempty(sProcess.options.overlap.Value{1})
        OPTIONS.CohOverlap = sProcess.options.overlap.Value{1}/100 ;
    end

    % Compute metric
    OutputFiles = bst_connectivity({sInputA.FileName}, {sInputA.FileName}, OPTIONS);
end

%% ===== TEST FUNCTION =====
function Test(iTest) %#ok<DEFNU>
    % Start a new report
    bst_report('Start');
    % Select tests
    if (nargin < 1) || isempty(iTest)
        iTest = 1;
    end
    % Get test datasets
    switch iTest
        case 1,  sFile = process_simulate_matrix('Test');
        case 2,  sFile = process_simulate_ar('Test');
    end
    % Loop on frequency resolutions
    for freq = [1 2 3 5 10 20]
        % Coherence process
        sTmp = bst_process('CallProcess', 'process_cohere1', sFile, [], ...
            'timewindow',   [], ...    % All the time in input
            'src_rowname',  '1', ...   % Test signal #1 => other signals
            'cohmeasure',   1, ...     % 1=Magnitude-squared, 2=Imaginary
            'overlap',      3, ...     % 50%
            'maxfreqres',   freq, ...  % VARIES
            'maxfreq',      [], ...    % No maximum frequency
            'pThresh',      0.05, ...  % p-value thrshold
            'outputmode',   1);        % Save individual results (one file per input file)
        % Snapshot: spectrum
        bst_process('CallProcess', 'process_snapshot', sTmp, [], ...
            'target',       10, ...  % Frequency spectrum
            'modality',     1, 'orient', 1, 'time', 0, 'contact_time', [-40, 110], 'contact_nimage', 16, ...
            'Comment',      [sFile.Comment, ': ' sTmp.Comment]);
    end
    % Save and display report
    ReportFile = bst_report('Save', sFile);
    bst_report('Open', ReportFile);
end

author: brainstorm-tools


text:
function [pdag, G] = learn_struct_pdag_pc_mod(cond_indep, n, k, varargin)
% LEARN_STRUCT_PDAG_PC Learn a partially oriented DAG (pattern) using the PC algorithm
% P = learn_struct_pdag_pc(cond_indep, n, k, ...)
%
% n is the number of nodes.
% k is an optional upper bound on the fan-in (default: n)
% cond_indep is a boolean function that will be called as follows:
%   feval(cond_indep, x, y, S, ...)
% where x and y are nodes, and S is a set of nodes (positive integers),
% and ... are any optional parameters passed to this function.
%
% The output P is an adjacency matrix, in which
% P(i,j) = -1 if there is an i->j edge.
% P(i,j) = P(j,i) = 1 if there is an undirected edge i <-> j
%
% The PC algorithm does structure learning assuming all variables are observed.
% See Spirtes, Glymour and Scheines, "Causation, Prediction and Search", 1993, p117.
% This algorithm may take O(n^k) time if there are n variables and k is the max fan-in,
% but this is quicker than the Verma-Pearl IC algorithm, which is always O(n^n).

sep = cell(n,n);
ord = 0;
done = 0;
G = ones(n,n);
G=setdiag(G,0);
while ~done
 done = 1;
 [X,Y] = find(G);
 for i=1:length(X)
   x = X(i); y = Y(i);
   %nbrs = mysetdiff(myunion(neighbors(G, x), neighbors(G,y)), [x y]);
   nbrs = mysetdiff(neighbors(G, y), x);  % bug fix by Raanan Yehezkel <raanany@ee.bgu.ac.il> 6/27/04
   if length(nbrs) >= ord & G(x,y) ~= 0
     done = 0;
     %SS = subsets(nbrs, ord, ord); % all subsets of size ord
     SS = subsets1(nbrs, ord);
     for si=1:length(SS)
       S = SS{si};
       if feval(cond_indep, x, y, S, varargin{:})
         %if isempty(S)
         %  fprintf('%d indep of %d ', x, y);
         %else
         %  fprintf('%d indep of %d given ', x, y); fprintf('%d ', S);
         %end
         %fprintf('\n');

         % diagnostic
         %[CI, r] = cond_indep_fisher_z(x, y, S, varargin{:});
         %fprintf(': r = %6.4f\n', r);

         G(x,y) = 0;
         G(y,x) = 0;
         sep{x,y} = myunion(sep{x,y}, S);
         sep{y,x} = myunion(sep{y,x}, S);
         break; % no need to check any more subsets
       end
     end
   end
 end
 ord = ord + 1;
end

% Create the minimal pattern,
% i.e., the only directed edges are V structures.
pdag = G;
[X, Y] = find(G);
% We want to generate all unique triples x,y,z
% This code generates x,y,z and z,y,x.
for i=1:length(X)
 x = X(i);
 y = Y(i);
 Z = find(G(y,:));
 Z = mysetdiff(Z, x);
 for z=Z(:)'
   if G(x,z)==0 & ~ismember(y, sep{x,z}) & ~ismember(y, sep{z,x})
     %fprintf('%d -> %d <- %d\n', x, y, z);
     pdag(x,y) = -1; pdag(y,x) = 0;
     pdag(z,y) = -1; pdag(y,z) = 0;
   end
 end
end

% Convert the minimal pattern to a complete one,
% i.e., every directed edge in P is compelled
% (must be directed in all Markov equivalent models),
% and every undirected edge in P is reversible.
% We use the rules of Pearl (2000) p51 (derived in Meek (1995))

old_pdag = zeros(n);
iter = 0;
while ~isequal(pdag, old_pdag)
 iter = iter + 1;
 old_pdag = pdag;
 % rule 1
 [A,B] = find(pdag==-1); % a -> b
 for i=1:length(A)
   a = A(i); b = B(i);
   C = find(pdag(b,:)==1 & G(a,:)==0); % all nodes adj to b but not a
   if ~isempty(C)
     pdag(b,C) = -1; pdag(C,b) = 0;
     %fprintf('rule 1: a=%d->b=%d and b=%d-c=%d implies %d->%d\n', a, b, b, C, b, C);
   end
 end
 % rule 2
 [A,B] = find(pdag==1); % unoriented a-b edge
 for i=1:length(A)
   a = A(i); b = B(i);
   if any( (pdag(a,:)==-1) & (pdag(:,b)==-1)' );
     pdag(a,b) = -1; pdag(b,a) = 0;
     %fprintf('rule 2: %d -> %d\n', a, b);
   end
 end
 % rule 3
 [A,B] = find(pdag==1); % a-b
 for i=1:length(A)
   a = A(i); b = B(i);
   % Bug fix by Imme Ebert-Uphoff (ebert@tree.com), Jan 2007
   % C = find( (G(a,:)==1) & (pdag(:,b)==-1)' );
   C = find( (pdag(a,:)==1) & (pdag(:,b)==-1)' );
   % C contains nodes c s.t. a-c->ba
   G2 = setdiag(G(C, C), 1);
   if any(G2(:)==0) % there are 2 different non adjacent elements of C
     pdag(a,b) = -1; pdag(b,a) = 0;
     %fprintf('rule 3: %d -> %d\n', a, b);
   end
 end
end

%  % Test Rule 3 of PC algorithm
%
%  % Define PDAG
%  pdag = zeros(4);
%  pdag(2,1)=-1;
%  pdag(3,1)=-1;
%  pdag(2,4)=-1;
%  pdag(3,4)=-1;
%  pdag(1,4)=1;
%  pdag(4,1)=1;
%
%  fprintf('\nSample input PDAG:\n');
%  pdag
%
%  fprintf('Sample DAG generated from PDAG:\n');
%  dag = abs(pdag_to_dag(pdag))
%
%  fprintf('Output from current PC algorithm:\n');
%  pdag_PC = learn_struct_pdag_pc('dsep', 4, 3, dag)
%
%  % Problem can be fixed by changing Line 120 of learn_struct_pdag_pc.m
%  %    C = find( (G(a,:)==1) & (pdag(:,b)==-1)' );
%  % to
%  %    C = find( (pdag(a,:)==1) & (pdag(:,b)==-1)' );
%
%  fprintf('Correct version:\n');
%  pdag_PC_mod = learn_struct_pdag_pc_mod('dsep', 4, 3, dag)

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function I = sum3(f)
%SUM3   Triple integral of a BALLFUN over its domain.
%   I = SUM3(F) returns the double definite integral of a BALLFUN.
%
% See also SUM, SUM2.

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

% Empty check
if ( isempty(f) )
    I = [];
    return
end

% Second argument: integral over the sphere of radius -1
[m,n,p] = size(f);
m = m+2; p = p+2;
F = coeffs3(f,m,n,p);

% Extract the 0-th Fourier mode
F = reshape(F(:,floor(n/2)+1,:), m, p);

% Increase the discretization by 2 in the r and theta direction
F = [zeros(m,1),F,zeros(m,1);zeros(2,p+2)];
m = m+2;
p = p+2;

% Multiply f par r^2sin(theta) (= Jacobian)
Msin = trigspec.multmat(p, [0.5i;0;-0.5i] );
Mr2 = ultraS.multmat(m, [0.5; 0; 0.5], 0 );
F = Mr2*F*(Msin.');

% Coefficients of integration between 0 and 1 of the chebyshev polynomials
IntChebyshev = zeros(1,m);
for i = 0:m-1
    if mod(i,4)==0
        IntChebyshev(i+1) = -1/(i^2-1);
    elseif mod(i,4)==1
        IntChebyshev(i+1) = 1/(i+1);
    elseif mod(i,4)==2
        IntChebyshev(i+1) = -1/(i^2-1);
    else
        IntChebyshev(i+1) = -1/(i-1);
    end
end

% Coefficients of integration between 0 and pi of the theta Fourier
% function
Listp = (1:p).' - floor(p/2)-1;
IntTheta = -1i*((-1).^Listp-1)./Listp;
IntTheta(floor(p/2)+1) = pi;

% Integrate over lambda
IntTheta = 2*pi*IntTheta;

% Integrate over the sphere of radius -1
if nargin>1
    IntChebyshev = (-1).^(0:m-1).*IntChebyshev;
end

% Return the integral of f over the ballfun
I = IntChebyshev*F*IntTheta;

% Return real value if the function is real
if f.isReal
   I = real(I);
end
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
1     % problem
4     % grid parameter
1     % uniform/stretched grid
2     % domain size
1     % approximation

%% Data file for test problem P1


text:
classdef testHorizontalTensorRotatedVsVHP < ...
        testHorizontalTensorRotatedVsSequentialLaminate

    methods (Access = public)

        function obj = testHorizontalTensorRotatedVsVHP()
            obj.computeTest()
        end

        function hasPassed = hasPassed(obj)
            rotHor = obj.rotHorTensor.getValue();
            lTens  = obj.lamTensor.getValue();
            lb = norm(rotHor - lTens) > 1e-12;
            ub = norm(rotHor - lTens) < 1e-2;
            hasPassed = lb & ub;
        end

    end

    methods (Access = protected)

        function computeLaminateDirectly(obj)
            c0       = obj.C0;
            c1       = obj.C1;
            dir{1}   = obj.lamDir;
            m1       = obj.lamPar;
            frac     = obj.theta;
            lam      = VoigtHomogPlaneStressHomogenizer(c0,c1,dir,m1,frac);
            obj.lamTensor = lam.getPlaneStressHomogenizedTensor();
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function dfadjusttoolbar(dffig)
%DFADJUSTTOOLBAR Adjust contents of distribution fitting plot toolbar

%   $Revision: 1.1.6.5 $  $Date: 2004/01/24 09:35:14 $
%   Copyright 2003-2004 The MathWorks, Inc.

h0 = findall(dffig,'Type','uitoolbar');
h1 = findall(h0,'Parent',h0);
czoom = [];
for j=length(h1):-1:1
   mlabel = get(h1(j),xlate('TooltipString'));
   if ~isempty(findstr(mlabel,'Zoom')) || ~isempty(findstr(mlabel,'Pan'))
      czoom(end+1) = h1(j);
   elseif isempty(findstr(mlabel,'Print'))
      delete(h1(j));
      h1(j) = [];
   else
     c1 = h1(j);
   end
end

% Add more icons especially for distribution fitting
if exist('dficons.mat','file')==2
   icons = load('dficons.mat','icons');
   state = dfgetset('showlegend');
   if isempty(state), state = 'on'; end
   try
      % Try to get the default MATLAB legend icon
      legicon = load([matlabroot '/toolbox/matlab/icons/legend.mat']);
      cdata = legicon.cdata;
   catch
      cdata = icons.icons.legend; % in case of trouble, use older icon
   end
   c2 = uitoggletool(h0, 'CData',cdata,...
                    'State',state,...
                    'TooltipString', 'Legend On/Off',...
                    'Separator','on',...
                    'ClickedCallback','dfittool(''togglelegend'')',...
                    'Tag','showlegend');
   state = dfgetset('showgrid');
   if isempty(state), state = 'off'; end
   c3 = uitoggletool(h0, 'CData',icons.icons.grid,...
                    'State',state,...
                    'TooltipString', ('Grid On/Off'),...
                    'Separator','off',...
                    'ClickedCallback','dfittool(''togglegrid'')',...
                    'Tag','showgrid');
   c = get(h0,'Children');
   cnew = [c1 czoom c2 c3]';

   set(h0,'Children',cnew(end:-1:1));
end

author: zouchuhang

repo: LayoutNet


text:
function prtPlotUtilAlgorithmGui(connectivityMatrix, actionCell, algo)
% Internal function, for PRT use only, makes use of GraphViz
% xxx NEED HELP xxx

% Given the algorithms action cell extract the names of each of the blocks
% Also add in Input and Output blocks
algoStr = cellfun(@(c)c.nameAbbreviation,actionCell,'uniformoutput',false);
algoStr = cat(1,{'Input'}, algoStr(:), {'Output'});
actionCell = cat(1,{[]},actionCell(:),{[]});

% Call GraphViz to get a good initial layout
GraphLayoutInfo = prtPlotUtilGraphVizRun(connectivityMatrix');

% Scale the GraphViz layout to get something reasonable
nodePosMat = cat(1,GraphLayoutInfo.Nodes.pos);
nodePosMat = bsxfun(@minus,nodePosMat,min(nodePosMat));
nodePosMat = bsxfun(@rdivide,nodePosMat,max(nodePosMat));
nodePosMat(isnan(nodePosMat))=0;

% Decide how big the blocks can reasonably be given the number of blocks
% we have to put in the unit square
minDistBetweenBlocks = min(min(prtDistanceLNorm(nodePosMat,nodePosMat,2) + realmax*eye(size(connectivityMatrix))));
blockSize = minDistBetweenBlocks*2/3;

% Because nodePosMat gives us the left corner we have to modify the block
% size a little
nodePosMat = bsxfun(@rdivide,nodePosMat,max(nodePosMat)+[blockSize blockSize]);

strLengths = cellfun(@(s)length(s),algoStr);

textSizes = blockSize./(strLengths+1)/0.7;
textSizes = ones(size(textSizes))*min(textSizes);
originalTextSizes = textSizes;
% Add one character for the sides
% 0.8 is an approximate to the aspect ratio of fonts.

% Get the block drawing options
Options = localGetOptions();

Handles = localMakeFigure();

% Ready the outputs
BlankBlock = struct('handle',[],'textHandle',[], 'Object',[]);

Layout.nBlocks = size(connectivityMatrix,1);
Layout.Blocks = BlankBlock;
Layout.connectivity = connectivityMatrix;
Layout.edges = zeros(sum(connectivityMatrix(:)),1);

Gui.status = '';
Gui.previous_point = [];

for iBlockOuter = 1:Layout.nBlocks
    placeBlockFunction([], [], actionCell{iBlockOuter}, nodePosMat(iBlockOuter,:), iBlockOuter, algoStr{iBlockOuter}, textSizes(iBlockOuter));
end

for iEdge = 1:length(GraphLayoutInfo.Edges)
    startLoc = cat(2,nodePosMat(GraphLayoutInfo.Edges(iEdge).startIndex,1)+blockSize,nodePosMat(GraphLayoutInfo.Edges(iEdge).startIndex,2)+blockSize/2);
    stopLoc = cat(2,nodePosMat(GraphLayoutInfo.Edges(iEdge).stopIndex,1),nodePosMat(GraphLayoutInfo.Edges(iEdge).stopIndex,2)+blockSize/2);

    %Layout.edges(iEdge) = prtPlotUtilPlotArrow(cat(1,startLoc(1),stopLoc(1)),cat(1,startLoc(2),stopLoc(2)),[20 4]);
    %set(Layout.edges(iEdge),'facecolor',[0 0 0],'edgecolor',[0 0 0])

    Layout.edges(iEdge) = prtPlotUtilPlotArrow(cat(1,startLoc(1),stopLoc(1)),cat(1,startLoc(2),stopLoc(2)),[],'headWidth',0.01,'realHeadLength',0.01);

end

centerBlocks();

set(gcf,'NextPlot','new')

%% Begin Functions

    function patchWindowButtonMotion(hObject, eventData, blockIndex, offSet) %#ok
        cp = get(Handles.AxesPanel.mainAxes,'currentPoint');

        oldBlockPosition = get(Layout.Blocks(blockIndex).handle,'position');

        cp = cp - blockSize/2 + repmat(cat(2, offSet, 0),2,1);

        Layout.Blocks(blockIndex).polygonPosition = cp(1,1:2);

        set(Layout.Blocks(blockIndex).handle,'position',cat(2,cp(1,1:2),oldBlockPosition(3:4)));
        set(Layout.Blocks(blockIndex).textHandle,'Position',[cp(1,1)+Options.blockSize(2)/2 cp(1,2)+Options.blockSize(2)/2 0]);

        %Move Edges
        for jEdge = 1:length(GraphLayoutInfo.Edges)
            if blockIndex == GraphLayoutInfo.Edges(jEdge).startIndex || blockIndex == GraphLayoutInfo.Edges(jEdge).stopIndex

                startInd = GraphLayoutInfo.Edges(jEdge).startIndex;
                stopInd = GraphLayoutInfo.Edges(jEdge).stopIndex;

                cStart = get(Layout.Blocks(startInd).handle,'position');
                cStop = get(Layout.Blocks(stopInd).handle,'position');

                startLoc = cStart(1:2) + [blockSize blockSize/2];
                stopLoc = cStop(1:2) + [0 blockSize/2];

%                 try  %#ok<TRYNC>
%                     delete(Layout.edges(jEdge));
%                 end
%                 Layout.edges(jEdge) = prtPlotUtilPlotArrow(cat(1,startLoc(1),stopLoc(1)),cat(1,startLoc(2),stopLoc(2)));

                %Layout.edges(jEdge) = prtPlotUtilPlotArrow(cat(1,startLoc(1),stopLoc(1)),cat(1,startLoc(2),stopLoc(2)),[20 4]);
                %set(Layout.edges(jEdge),'facecolor',[0 0 0],'edgecolor',[0 0 0])

                Layout.edges(jEdge) = prtPlotUtilPlotArrow(cat(1,startLoc(1),stopLoc(1)),cat(1,startLoc(2),stopLoc(2)),Layout.edges(jEdge));
            end
        end
    end

    function patchButtonDownFunction(hObject, eventData, blockIndex) %#ok<INUSL>

        switch lower(get(Handles.handle, 'selectiontype'))
            case 'open'
                % Nothing for now
                BlockObject = Layout.Blocks(blockIndex).Object;
                if ~isempty(BlockObject) && BlockObject.isTrained && prtUtilIsMethodIncludeHidden(BlockObject,'plot') && BlockObject.dataSetSummary.nFeatures < 4
                    plot(BlockObject);
                elseif (blockIndex == length(Layout.Blocks)) && algo.isPlottableAsClassifier
                    plotAsClassifier(algo);
                elseif (blockIndex == length(Layout.Blocks)) && algo.isPlottableAsRegressor
                    plotAsRegressor(algo);
                end

            otherwise % case 'normal'

                % Get block clock offset
                cp = get(Handles.AxesPanel.mainAxes,'currentPoint');

                cPolyCorner = Layout.Blocks(blockIndex).polygonPosition;

                offSet = cPolyCorner + blockSize/2 - cp(1,1:2);

                set(Handles.handle,'WindowButtonMotionFcn',@(h,E)patchWindowButtonMotion(h,E, blockIndex, offSet));
                set(Handles.handle,'WindowButtonUpFcn',@(h,E)patchButtonUpFunction(h,E, blockIndex));
                set(Layout.Blocks(blockIndex).handle,'Selected','off')
        end
    end

    function patchButtonUpFunction(hObject, eventData, blockIndex)   %#ok<INUSL>
        set(Handles.handle,'WindowButtonMotionFcn',@mouseControlWindowButtonMotionFcn,...
                           'WindowButtonDownFcn'  , @mouseControlWindowButtonDownFcn,...
                           'WindowButtonUpFcn'    , @mouseControlWindowButtonUpFcn);
        Gui.status = '';

        set(Layout.Blocks(blockIndex).handle,'Selected','off')
    end

    function placeBlockFunction(hObject, eventData, BlockObject, position, iBlock, blockStr, textSize)  %#ok<INUSL>

        % Extract stuff from the options so that the GUI can use it
        if isempty(BlockObject)
            blockColor = Options.BlockColors.dataSet;
            textColor = Options.BlockTextColors.dataSet;
        elseif isa(BlockObject,'prtPreProc')
            blockColor = Options.BlockColors.preProcessor;
            textColor = Options.BlockTextColors.preProcessor;
        elseif isa(BlockObject,'prtFeatSel')
            blockColor = Options.BlockColors.featureSelector;
            textColor = Options.BlockTextColors.featureSelector;
        elseif isa(BlockObject,'prtClass')
            blockColor = Options.BlockColors.classifier;
            textColor = Options.BlockTextColors.classifier;
        elseif isa(BlockObject,'prtDecision')
            blockColor = Options.BlockColors.decision;
            textColor = Options.BlockTextColors.decision;
        elseif isa(BlockObject,'prtDecision')
            blockColor = Options.BlockColors.decision;
            textColor = Options.BlockTextColors.decision;
        elseif isa(BlockObject,'prtOutlierRemoval')
            blockColor = Options.BlockColors.outlierRemoval ;
            textColor = Options.BlockTextColors.outlierRemoval ;
        elseif isa(BlockObject,'prtCluster')
            blockColor = Options.BlockColors.cluster ;
            textColor = Options.BlockTextColors.cluster ;
        elseif isa(BlockObject,'prtRv')
            blockColor = Options.BlockColors.rv ;
            textColor = Options.BlockTextColors.rv ;
        else
            blockColor = [1 1 1];
            textColor = [0 0 0];
        end

        Layout.Blocks(iBlock).polygonPosition = position;
        Layout.Blocks(iBlock).handle = rectangle('Position',[position,blockSize,blockSize],'Curvature',[0.25, 0.25],'FaceColor',blockColor,'LineWidth',2);

        if (~isempty(BlockObject) && BlockObject.isTrained && prtUtilIsMethodIncludeHidden(BlockObject,'plot') && BlockObject.dataSetSummary.nFeatures < 4)
            %set(Layout.Blocks(iBlock).handle,'lineWidth',2,'EdgeColor',[1 1 0.1]);
            set(Layout.Blocks(iBlock).handle,'lineWidth',4);
        end

        if (iBlock == Layout.nBlocks) && (algo.isPlottableAsClassifier || algo.isPlottableAsRegressor)
            set(Layout.Blocks(iBlock).handle,'lineWidth',4);
        end

        Layout.Blocks(iBlock).textHandle = text(position(1)+blockSize/2,position(2)+blockSize/2,blockStr,'VerticalAlignment','Middle','HorizontalAlignment','Center','Color',textColor,'FontUnits','Normalized','FontSize',textSize,'FontWeight','Bold');

        set(Layout.Blocks(iBlock).handle,'ButtonDownFcn',@(h,E)patchButtonDownFunction(h,E,iBlock));
        set(Layout.Blocks(iBlock).textHandle,'ButtonDownFcn',@(h,E)patchButtonDownFunction(h,E,iBlock));

        Layout.Blocks(iBlock).Object = BlockObject;
    end

    function mouseControlWindowScrollWheelZoomFcn(varargin)
        scrollAmountFactor = 50;

        if ~ishandle(Handles.handle)
            return
        end

        % current axes?
        if get(Handles.handle, 'currentaxes') ~= Handles.AxesPanel.mainAxes
            return
        end

        % get the amount of scolls
        scrolls = varargin{2}.VerticalScrollCount;

        % Get the axes' x- and y-limits
        xlim = get(Handles.AxesPanel.mainAxes, 'xlim');
        ylim = get(Handles.AxesPanel.mainAxes, 'ylim');

        % get the current camera position, and save the [z]-value
        cam_pos_Z = get(Handles.AxesPanel.mainAxes, 'cameraposition');
        cam_pos_Z = cam_pos_Z(3);

        % get the current point
        old_position = get(Handles.AxesPanel.mainAxes, 'CurrentPoint');
        old_position(1,3) = cam_pos_Z;

        % calculate zoom factor
        zoomfactor = 1 - scrolls/scrollAmountFactor;

        % adjust camera position
        set(Handles.AxesPanel.mainAxes,...
            'cameratarget', [old_position(1, 1:2), 0],...
            'cameraposition', old_position(1, 1:3));

        % adjust the camera view angle (equal to zooming in)
        camzoom(zoomfactor);

        % zooming with the camera has the side-effect of
        % NOT adjusting the axes limits. We have to correct for this:
        x_lim1 = (old_position(1,1) - min(xlim))/zoomfactor;
        x_lim2 = (max(xlim) - old_position(1,1))/zoomfactor;
        xlim   = [old_position(1,1) - x_lim1, old_position(1,1) + x_lim2];
        y_lim1 = (old_position(1,2) - min(ylim))/zoomfactor;
        y_lim2 = (max(ylim) - old_position(1,2))/zoomfactor;
        ylim   = [old_position(1,2) - y_lim1, old_position(1,2) + y_lim2];
        set(Handles.AxesPanel.mainAxes, 'xlim', xlim)
        set(Handles.AxesPanel.mainAxes, 'ylim', ylim)

        % set new camera position
        new_position = get(Handles.AxesPanel.mainAxes, 'CurrentPoint');
        old_camera_target =  get(Handles.AxesPanel.mainAxes, 'CameraTarget');
        old_camera_target(3) = cam_pos_Z;
        new_camera_position = old_camera_target - ...
            (new_position(1,1:3) - old_camera_target(1,1:3));

        % adjust camera target and position
        set(Handles.AxesPanel.mainAxes, 'cameraposition', new_camera_position(1, 1:3),...
            'cameratarget', [new_camera_position(1, 1:2), 0]);

        % we also have to re-set the axes to stretch-to-fill mode
        set(Handles.AxesPanel.mainAxes, 'cameraviewanglemode', 'auto',...
            'camerapositionmode', 'auto',...
            'cameratargetmode', 'auto');

        % Resize the text since the block size has changed
        xRange = xlim(2)-xlim(1);
        yRange = ylim(2)-ylim(1);

        zoomFactorFromOriginal = max(xRange,yRange) ./ (Options.initialCanvasLimits(2)-Options.initialCanvasLimits(1));
        textSizes = 1./zoomFactorFromOriginal.*originalTextSizes;

        resizeText();
    end

    function mouseControlWindowButtonDownFcn(varargin)
        % If we are not a handle get out of here
        if ~ishandle(Handles.handle)
            return
        end

        % current axes?
        if get(Handles.handle, 'currentaxes') ~= Handles.AxesPanel.mainAxes
            return
        end

        % perform appropriate action
        switch lower(get(Handles.handle, 'selectiontype'))
            % start panning on left click
            case 'normal'
                cp = get(Handles.AxesPanel.mainAxes, 'CurrentPoint');

                deselectEverything();
                % Clicking on back canvas Allow dragging around
                Gui.status = 'down';
                Gui.previous_point = cp;

            case 'open' % double click (left or right)
                centerBlocks();
                Gui.status = '';
        end
    end

    function mouseControlWindowButtonUpFcn(varargin)
        % double check if these axes are indeed the current axes
        if get(Handles.handle, 'currentaxes') ~= Handles.AxesPanel.mainAxes
            return
        end

        %deselectEverything();
        % reset Gui.status
        Gui.status = '';
    end

    % The primary mouse button controler
    function mouseControlWindowButtonMotionFcn(varargin)
        if ~ishandle(Handles.handle)
            return
        end

        % double check if these axes are indeed the current axes
        if get(Handles.handle, 'currentaxes') ~= Handles.AxesPanel.mainAxes
            return
        end
        % return if there isn't a previous point
        if isempty(Gui.previous_point)
            return
        end
        % return if mouse hasn't been clicked
        if isempty(Gui.status)
            return
        end
        % get current location (in pixels)
        current_point = get(Handles.AxesPanel.mainAxes, 'CurrentPoint');

        switch lower(Gui.status)
            case 'down'
                % get current XY-limits
                xlim = get(Handles.AxesPanel.mainAxes, 'xlim');
                ylim = get(Handles.AxesPanel.mainAxes, 'ylim');

                % find change in position
                delta_points = current_point - Gui.previous_point;
                % adjust limits
                new_xlim = xlim - delta_points(1);
                new_ylim = ylim - delta_points(3);
                % set new limits
                set(Handles.AxesPanel.mainAxes, 'Xlim', new_xlim); set(Handles.AxesPanel.mainAxes, 'Ylim', new_ylim);
                % save new position
                Gui.previous_point = get(Handles.AxesPanel.mainAxes, 'CurrentPoint');
            otherwise % Window motion without anything going on
        end
    end

    function deselectEverything()
        for iBlock = 1:Layout.nBlocks
            set(Layout.Blocks(iBlock).handle,'Selected','off');
        end
    end

    function centerBlocks()
        if Layout.nBlocks == 0
            xLims = Options.initialCanvasLimits(1:2);
            yLims = Options.initialCanvasLimits(3:4);
            Options.blockFontSizeNormalized = Options.initialBlockFontSizeNormalized;

            xlim(xLims)
            ylim(yLims);

            return
        end

        xLims = [inf -inf];
        yLims = [inf -inf];
        for iBlock = 1:Layout.nBlocks
            maxPos = Layout.Blocks(iBlock).polygonPosition + blockSize;
            minPos = Layout.Blocks(iBlock).polygonPosition;
            xLims(1) = min(xLims(1),minPos(1));
            xLims(2) = max(xLims(2),maxPos(1));
            yLims(1) = min(yLims(1),minPos(2));
            yLims(2) = max(yLims(2),maxPos(2));
        end
        % Add some percentage to the edges
        edgePadPercent = 0.05;
        xLimsPadded = [-(xLims(2)-xLims(1)) (xLims(2)-xLims(1))]*edgePadPercent + xLims;
        yLimsPadded = [-(yLims(2)-yLims(1)) (yLims(2)-yLims(1))]*edgePadPercent + yLims;

        xRange = xLimsPadded(2)-xLimsPadded(1);
        yRange = yLimsPadded(2)-yLimsPadded(1);

        if xRange > yRange
            % We need to add some to the yRange to make equal aspect ratio
            xLims = xLimsPadded;
            yLims = yLimsPadded + (xRange-yRange)/2*[-1 1];
            yRange = xRange;
        else
            xLims = xLimsPadded + (yRange-xRange)/2*[-1 1];
            yLims = yLimsPadded;
            xRange = yRange;
        end

        zoomFactorFromOriginal = max(xRange,yRange) ./ (Options.initialCanvasLimits(2)-Options.initialCanvasLimits(1));
        textSizes = 1./zoomFactorFromOriginal.*originalTextSizes;
        resizeText();

        xlim(xLims)
        ylim(yLims);
    end

    function resizeText
        for iBlock = 1:Layout.nBlocks
            set(Layout.Blocks(iBlock).textHandle,'FontSize',textSizes(iBlock));
        end
    end

%     function blockInd = findBlockClick(cp)
%         blockInd = [];
%         for iBlock = 1:Layout.nBlocks
%             isInThisBlock = all(cp(1,1:2) < max(Layout.Blocks(iBlock).polygonNodes) & cp(1,1:2) > min(Layout.Blocks(iBlock).polygonNodes));
%             if isInThisBlock
%                 blockInd = iBlock;
%
%                 return
%             end
%         end
%     end

    function Options = localGetOptions()

        colors = prtPlotUtilClassColors(7);

        Options.initialCanvasLimits = [-0.1 1.1 -0.1 1.1];
        Options.blockSize = blockSize*[1 1];
        Options.initialBlockFontSizeNormalized = 0.03; %This is the default relative to the initialCanvas limits
        Options.blockFontSizeNormalized = Options.initialBlockFontSizeNormalized; %It can change as a function of the zoom level.

        Options.BlockColors.dataSet = [0.8 0.8 0.8];
        Options.BlockColors.preProcessor = colors(1,:);
        Options.BlockColors.featureSelector = colors(2,:);
        Options.BlockColors.classifier = colors(3,:);
        Options.BlockColors.decision = colors(4,:);
        Options.BlockColors.outlierRemoval = colors(5,:);
        Options.BlockColors.cluster = colors(6,:);
        Options.BlockColors.rv = colors(7,:);

        Options.BlockTextColors.dataSet = [0 0 0];
        Options.BlockTextColors.preProcessor = [1 1 1];
        Options.BlockTextColors.featureSelector = [1 1 1];
        Options.BlockTextColors.classifier = [1 1 1];
        Options.BlockTextColors.decision = [0 0 0];
        Options.BlockTextColors.rv = [0 0 0];
        Options.BlockTextColors.outlierRemoval = [1 1 1];
        Options.BlockTextColors.cluster = [1 1 1];

    end

    function Handles = localMakeFigure()

        ss = get(0,'screensize');

        windowSize = [754 600];

        % Center the window
        sizePads = round((ss(3:4)-windowSize));
        sizePads(1) = sizePads(1)/2; % We should use 2 right?
        sizePads(2) = sizePads(2)/2;
        pos = cat(2,sizePads,windowSize);

        % Create the figure an UIControls
        MainFigure.style = 'figure';
        MainFigure.units = 'pixels';
        MainFigure.position = pos;
        MainFigure.name = 'PRT Algorithm';
        MainFigure.NumberTitle = 'Off';
        MainFigure.Menu = 'none';
        MainFigure.toolbar = 'none';
        MainFigure.DockControls = 'off';

        MainFigure.Children.AxesPanel.style = 'panel';
        MainFigure.Children.AxesPanel.units = 'normalized';
        MainFigure.Children.AxesPanel.position = [0 0 1 1];
        MainFigure.Children.AxesPanel.Children.mainAxes.style = 'axes';
        MainFigure.Children.AxesPanel.Children.mainAxes.units = 'Normalized';
        MainFigure.Children.AxesPanel.Children.mainAxes.position = [0 0 1 1];
        MainFigure.Children.AxesPanel.Children.mainAxes.xtick = [];
        MainFigure.Children.AxesPanel.Children.mainAxes.ytick = [];
        MainFigure.Children.AxesPanel.Children.mainAxes.xlim = [0 1];
        MainFigure.Children.AxesPanel.Children.mainAxes.ylim = [0 1];
        MainFigure.Children.AxesPanel.Children.mainAxes.nextPlot = 'add';

        Handles = prtUtilSuicontrol(MainFigure);

        % define zooming with scrollwheel, and panning with mouseclicks
        set(Handles.handle, 'WindowScrollWheelFcn' , @mouseControlWindowScrollWheelZoomFcn,...
            'WindowButtonDownFcn'  , @mouseControlWindowButtonDownFcn,...
            'WindowButtonUpFcn'    , @mouseControlWindowButtonUpFcn,...
            'WindowButtonMotionFcn', @mouseControlWindowButtonMotionFcn);

    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function X = polar2xyz(P)
% POLAR2XYZ Transformation of polar coordinates to cartesian coordinates.
% ------------------------------------------------------------------------------
% DESCRIPTION/NOTES
% Transformation of polar coordinates [distance, horizontal angle, vertical
% angle] to cartesian coordinates [x, y, z].
% ------------------------------------------------------------------------------
% INPUT
% P
%   Matrix of polar coordinates distace, horizontal angle and vertical angle.
%   Each point in a row.
%   Horizontal angle is defined mathematically positive starting from x-axis.
%   Vertical   angle is defined as angle between point and zenith.
% ------------------------------------------------------------------------------
% OUTPUT
% X
%   Matrix of cartesian coordinates x, y and z. Each point in a row.
% ------------------------------------------------------------------------------
% EXAMPLES
% Transformation of randomly generated points into polar coordinates and back.
%   X1   = [rand(10,1)*10 rand(10,1)*5 rand(10,1)];
%   P    = xyz2polar(X1);
%   X2   = polar2xyz(P);
%   Diff = X1-X2;
% ------------------------------------------------------------------------------
% pg@geo.tuwien.ac.at
% ------------------------------------------------------------------------------

% Input parsing ----------------------------------------------------------------

p = inputParser;
p.addRequired('P', @(x) isreal(x) && size(x,2) == 3);
p.parse(P);
p = p.Results;
% Clear required inputs to avoid confusion
clear P

% Transformation ---------------------------------------------------------------

d  = p.P(:,1);
ha = p.P(:,2);
va = p.P(:,3);

x = d .* cosg(ha) .* sing(va);
y = d .* sing(ha) .* sing(va);
z = d .* cosg(va);

% Output -----------------------------------------------------------------------

X = [x y z];

end

author: pglira


text:
function [dat] = read_neuralynx_cds(filename, hdr, begsample, endsample, chanindx);

% READ_NEURALYNX_CDS reads selected samples and channels from a  combined Neuralynx dataset with seperate subdirectories for the LFP, MUA and spike channels
%
% Use as
%    hdr = read_neuralynx_cds(parentdir)
%    dat = read_neuralynx_cds(parentdir, hdr, begsample, endsample, chanindx)

% Copyright (C) 2006, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id: read_neuralynx_cds.m 945 2010-04-21 17:41:20Z roboos $

needhdr = (nargin==1);
needdat = (nargin>=2);

dirlist = dir(filename);
haslfp   = any(filetype_check_extension({dirlist.name}, 'lfp'));
hasmua   = any(filetype_check_extension({dirlist.name}, 'mua'));
hasspike = any(filetype_check_extension({dirlist.name}, 'spike'));
%hasttl  = any(filetype_check_extension({dirlist.name}, 'ttl'));   % seperate file with original Parallel_in
%hastsl  = any(filetype_check_extension({dirlist.name}, 'tsl'));   % seperate file with original TimeStampLow
%hastsh  = any(filetype_check_extension({dirlist.name}, 'tsh'));   % seperate file with original TimeStampHi

if needhdr
  % read the header from the LFP dataset
  if haslfp
    lfpdir = fullfile(filename, dirlist(find(filetype_check_extension({dirlist.name}, 'lfp'))).name);
    lfphdr = read_header(lfpdir);
    lfpfil = lfphdr.filename;
    for i=1:lfphdr.nChans
      lfplab{i} = sprintf('lfp_%s', lfphdr.label{i});
    end
  else
    lfphdr = [];
    lfpfil = {};
    lfplab = {};
  end

  % read the header from the MUA dataset
  if hasmua
    muadir = fullfile(filename, dirlist(find(filetype_check_extension({dirlist.name}, 'mua'))).name);
    muahdr = read_header(muadir);
    muafil = muahdr.filename;
    for i=1:muahdr.nChans
      mualab{i} = sprintf('mua_%s', muahdr.label{i});
    end
  else
    muahdr = [];
    muafil = {};
    mualab = {};
  end

  % read the header from the SPIKE dataset
  if hasspike
    spikedir = fullfile(filename, dirlist(find(filetype_check_extension({dirlist.name}, 'spike'))).name);
    spikehdr = read_header(spikedir);
    spikefil = spikehdr.filename;
    for i=1:spikehdr.nChans
      spikelab{i} = sprintf('spike_%s', spikehdr.label{i});
    end
  else
    spikehdr = [];
    spikefil = {};
    spikelab = {};
  end

  if haslfp
    % assume that the LFP header describes the MUA and the spikes as well
    hdr = lfphdr;
  else hasmua
    % assume that the MUA header describes the spikes as well
    hdr = muahdr;
  end

  % concatenate the lfp/mua/spike channels, they have already been renamed
  hdr.label    = cat(1, lfplab(:), mualab(:), spikelab(:));
  hdr.filename = cat(1, lfpfil(:), muafil(:), spikefil(:));
  hdr.nChans = length(hdr.label);
  % remember the original header details
  hdr.orig    = {};
  hdr.orig{1} = lfphdr;
  hdr.orig{2} = muahdr;
  hdr.orig{3} = spikehdr;

  % return the header
  dat = hdr;
else

  % use the original header details
  lfphdr   = hdr.orig{1};
  muahdr   = hdr.orig{2};
  spikehdr = hdr.orig{3};
  % the spike header does not contain these, but requires them from the LFP or MUA
  spikehdr.FirstTimeStamp     = hdr.FirstTimeStamp;
  spikehdr.TimeStampPerSample = hdr.TimeStampPerSample;

  if ~isempty(lfphdr)
    Nlfp = lfphdr.nChans;
  else
    Nlfp = 0;
  end

  if ~isempty(muahdr)
    Nmua = muahdr.nChans;
  else
    Nmua = 0;
  end

  if ~isempty(spikehdr)
    Nspike = spikehdr.nChans;
  else
    Nspike = 0;
  end

  % determine which files should be read from each respective dataset
  selchan = zeros(hdr.nChans,1);
  selchan(chanindx) = 1;
  sellfp   = find(selchan(1:Nlfp));
  selmua   = find(selchan((Nlfp+1):(Nlfp+Nmua)));
  selspike = find(selchan((Nlfp+Nmua+1):(Nlfp+Nmua+Nspike)));

  if haslfp
    % the header contains the relevant file names
    lfpdat = read_neuralynx_ds([], lfphdr, begsample, endsample, sellfp);
  else
    lfpdat = [];
  end

  if hasmua
    % the header contains the relevant file names
    muadat = read_neuralynx_ds([], muahdr, begsample, endsample, selmua);
  else
    muadat = [];
  end

  if hasspike
    % the header contains the relevant file names
    spikedat = read_neuralynx_ds([], spikehdr, begsample, endsample, selspike);
  else
    spikedat = [];
  end

  % combine all data into a single array
  dat = cat(1, lfpdat, muadat, spikedat);

end

author: PatternRecognition

repo: OpenBMI


text:
function Cba=quat2dcm(qba)

a=qba(1);
b=qba(2);
c=qba(3);
d=qba(4);

Cba=zeros(3);
Cba(1,1)=a^2+b^2-c^2-d^2;
Cba(2,1)=2*(b*c+a*d);
Cba(3,1)=2*(b*d-a*c);

Cba(1,2)=2*(b*c-a*d);
Cba(2,2)=a^2-b^2+c^2-d^2;
Cba(3,2)=2*(c*d+a*b);

Cba(1,3)=2*(b*d+a*c);
Cba(2,3)=2*(c*d-a*b);
Cba(3,3)=a^2-b^2-c^2+d^2;

repo: ekfmonoslam


text:
% <Purpose>
%    Computing the numerical rank of a downdated matrix.
%
% <Syntax>
%    [r,Basis,C] = NumericalRankDowndate(A,pth,C,RC)
%
% <Input Parameters>
%  1.   A -- the target matrix;
%  2. pth -- the index of row/column to be deleted;
%  3.   C -- cell array contains information required by updating/downdating;
%  4.  RC -- Set to 'row', then the pth row will be deleted.
%            Set to 'column', then the pth column will be deleted.
%
% <Output Parameters>
%  1. r     -- the numerical rank of the downdated matrix;
%  2. Basis --
%         For high rank cases...
%             a matrix whose columns form an orthonormal basis of
%             the numerical kernel;
%
%         For low rank cases...
%             a matrix whose columns form an orthonormal basis of
%             the numerical range;
%
%  3.     C -- Matlab cell array
%
%     For high rank cases...
%         C{1,1} = rank : the numerical rank of the downdated matrix;
%         C{2,1} = Basis : matrix whose columns form an orthonormal kernel basis;
%         C{3,1} = Q : the Q in the QR decomposition of the kernel stacked matrix;
%         C{4,1} = R : the R in the QR decomposition of the kernel stacked matrix;
%         C{5,1} = tau : scaling factor in the kernel stacked matrix;
%         C{6,1} = tol : the rank decision threshold;
%
%     For low rank cases...
%         C{1,1} = rank : the numerical rank of the downdated matrix;
%         C{2,1} = U : the U in the USV+E decomposition of the downdated matrix;
%         C{3,1} = V : the V in the USV+E decomposition of the downdated matrix;
%         C{4,1} = S : the S in the USV+E decomposition of the downdated matrix;
%         C{5,1} = tol : the rank decision threshold;
%
% <Reference>
%   [1] T.Y. Li and Z. Zeng, "A Rank-Revealing Method with Updating, Downdating
%       and Applications", SIAM J. Matrix Anal. and Appl., 26 (2005), pp. 918--946.
%
%   [2] T.L. Lee, T.Y. Li and Z. Zeng, "A Rank-Revealing Method with Updating,
%       Downdating and Applications, Part II", SIAM J. Matrix Anal. and Appl. (2009).

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
% funJn calculates and returns the In integral (eq (15) in [1])
% Jn=funJn(Tn,betaConst,n,numbMC)
% Jn is the J_{n,\beta} integral (scalar) value
% T_n is the n-based SINR threshold value (eq (17) in [1])
% betaConstant is path-loss exponent
% n is integer parameter
% betaConst, n, and numbMC are scalars. T_n can be a vector
% numbMc is number of sample (Sobol) points for  quasi-MC integration
%
% Author: H.P. Keeler, Inria Paris/ENS, 2013
%
% References
% [1] H.P. Keeler, B. Błaszczyszyn and M. Karray,
% 'SINR-based k-coverage probability in cellular networks with arbitrary
% shadowing', accepted at ISIT, 2013

function Jn=funJn(Tn,betaConst,n,numbMC)
% Calculates Jn with various integration methods
% n =2 and 3 uses quad and dblquad respectively
% n>3 uses quasi Monte Carlo based on Sobol points
% function is called by funProbCov; see Corollary 7 in [1]
if nargin==3
    numbMC=10^3; %set default number of qMC sample points
end
if n<=3
    numbMC=0; %monte carlo not used
end

Jn=zeros(size(Tn));
for k=1:length(Tn)
    %%% Use  quadrature methods for n=2 and n=3 cases
    if n==3
        fv=@(v1,v2)(1./((v1.*v2)+Tn(k))).*(1./((v1.*(1-v2))+Tn(k))).*(v1.*v2.*(1-v2).*(1-v1)).^(2/betaConst).*v1.^(2/betaConst+1);
        Jn(k)=dblquad(fv,0,1,0,1); %perform double qudarature
    elseif n==2

        fv=@(v1)(v1.*(1-v1)).^(2/betaConst)./(v1+Tn(k));
        Jn(k)=quad(fv,0,1); %perform single qudarature

    elseif n==1
        Jn=ones(size(Tn)); %return ones since J_1=1;
    else
        %%% Use QMC method
        numbMCn=(n-1)*numbMC; %scale number of points by dimension
        cubeVol=1; %hyper-cube volume
        eta_i=cell(1,n);
        %Use sobol points (can also use 'halton')
        q = qrandstream('halton',n-1,'Skip',1e3,'Leap',1e2);
        qRandAll=qrand(q,numbMCn);

        %create eta_i values
        eta_i{1}=prod((qRandAll(:,1:end)),2);
        for i=2:n-1
            eta_i{i}=(1-qRandAll(:,i-1)).*prod((qRandAll(:,i:end)),2);
        end
        eta_i{n}=(1-qRandAll(:,n-1));

        %create/sample nominator and denominator of integral kernel
        numProdv_i=ones(numbMCn,1);
        denomProdv_i=ones(numbMCn,1);
        for i=1:n-1
            viRand=qRandAll(:,i);
            numProdv_i=(viRand).^(i*(2/betaConst+1)-1).*(1-viRand).^(2/betaConst).*numProdv_i; %numerator
            denomProdv_i=(eta_i{i}+Tn(k)).*denomProdv_i; %denominator term
        end
        denomProdv_i=(eta_i{n}+Tn(k)).*denomProdv_i;

        %factor out one term, arbitrarily choose the j-th term
        j=n-1;
        denomProdv_i=(eta_i{j}+Tn(k))./denomProdv_i;

        kernelInt=numProdv_i.*denomProdv_i; %integral kernerl
        Jn(k)=mean(kernelInt)*cubeVol; % perform (q)MC step

    end

end

author: Sable


text:
classdef MimContextResultCombiner < CoreBaseClass
    % MimContextResultCombiner. Part of the internal framework of the TD MIM Toolkit.
    %
    %     This class is used to build up aggregate results combined from results for different contexts
    %
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties (Access = private)
        FirstRun = true
        CacheInfo
        ImageTemplates
        OutputImage
        OutputImageTemplate
        Result
        PluginHasBeenRun = false
    end

    methods
        function obj = MimContextResultCombiner(image_templates)
            obj.ImageTemplates = image_templates;
        end

        function AddParentResult(obj, parent_context, this_context_results, output_context, output_context_mapping, dataset_stack, reporting)
            obj.Result = MimReduceResultToContext(this_context_results.GetResult, output_context_mapping.Context, obj.ImageTemplates, dataset_stack, reporting);
            if ~isempty(this_context_results.GetOutputImage)
                obj.OutputImage = MimReduceResultToContext(this_context_results.GetOutputImage, output_context, obj.ImageTemplates, dataset_stack, reporting);
            else
                obj.OutputImage = [];
            end
            obj.CacheInfo = this_context_results.GetCacheInfo;
            obj.PluginHasBeenRun = obj.PluginHasBeenRun || this_context_results.GetPluginHasBeenRun;
        end

        function AddChildResult(obj, child_context, this_context_results, output_context, output_context_mapping, dataset_stack, reporting)
            if isempty(obj.OutputImageTemplate)
                obj.OutputImageTemplate = obj.ImageTemplates.GetTemplateImage(output_context, dataset_stack, reporting);
            end

            obj.Result = MimExpandResultToContext(obj.Result, this_context_results.GetResult, child_context, obj.OutputImageTemplate, obj.ImageTemplates, dataset_stack, reporting);
            if ~isempty(this_context_results.GetOutputImage)
                obj.OutputImage = MimExpandResultToContext(obj.OutputImage, this_context_results.GetOutputImage, child_context, obj.OutputImageTemplate, obj.ImageTemplates, dataset_stack, reporting);
            else
                obj.OutputImage = [];
            end
            if isempty(obj.CacheInfo)
                obj.CacheInfo = MimCompositeResult;
            end
            obj.CacheInfo.AddField(char(child_context), this_context_results.GetCacheInfo);
            obj.PluginHasBeenRun = obj.PluginHasBeenRun || this_context_results.GetPluginHasBeenRun;
        end

        function result = GetResult(obj)
            result = obj.Result;
        end

        function result = GetCacheInfo(obj)
            result = obj.CacheInfo;
        end

        function result = GetOutputImage(obj)
            result = obj.OutputImage;
        end

        function has_been_run = GetPluginHasBeenRun(obj)
            has_been_run = obj.PluginHasBeenRun;
        end

    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
%*********************************************************
%% svec: compute the vector svec(M),
%%
%%   x = svec(blk,M,isspx);
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%**********************************************************

  function x = svec(blk,M,isspx);

   if iscell(M)
      if (size(blk,1) ~= size(M,1))
         error('svec: number of rows in blk and M not equal');
      end
      if (nargin == 2)
         %%if (size(M,2) == 1)
         %%   isspx = zeros(size(blk,1),1);
         %%else
         %%   isspx = ones(size(blk,1),1);
         %%end
         isspx = ones(size(blk,1),1);
      else
         if (length(isspx) < size(blk,1))
            isspx = ones(size(blk,1),1);
         end
      end
      x = cell(size(blk,1),1);
      for p=1:size(blk,1)
         pblk = blk(p,:);
         n = sum(pblk{2});  m = size(M,2);
         if strcmp(pblk{1},'s')
            n2 = sum(pblk{2}.*(pblk{2}+1))/2;
            if (isspx(p));
               x{p} = sparse(n2,m);
            else
               x{p} = zeros(n2,m);
            end
            numblk = length(pblk{2});
            if (pblk{2} > 0)
               for k = 1:m
                  if (numblk > 1) & ~issparse(M{p,k});
                     x{p}(:,k) = mexsvec(pblk,sparse(M{p,k}),isspx(p));
                  else
                     x{p}(:,k) = mexsvec(pblk,M{p,k},isspx(p));
                  end
               end
            end
	 else
            if (isspx(p))
               x{p} = sparse(n,m);
            else
               x{p} = zeros(n,m);
            end
            for k = 1:m
                x{p}(:,k) = M{p,k};
            end
         end
      end
   else
      if strcmp(blk{1},'s')
         numblk = length(blk{2});
         if (numblk > 1) & ~issparse(M);
            x = mexsvec(blk,sparse(M),1);
         else
            x = mexsvec(blk,sparse(M));
         end
      else
         x = M;
      end
   end
%%**********************************************************

repo: REST


text:
function [h, hdata] = mlphess_weighted(net, x, t, eso_w, hdata)
%MLPHESS Evaluate the Hessian matrix for a multi-layer perceptron network.
%
%	Description
%	H = MLPHESS(NET, X, T) takes an MLP network data structure NET, a
%	matrix X of input values, and a matrix T of target values and returns
%	the full Hessian matrix H corresponding to the second derivatives of
%	the negative log posterior distribution, evaluated for the current
%	weight and bias values as defined by NET.
%
%	[H, HDATA] = MLPHESS(NET, X, T) returns both the Hessian matrix H and
%	the contribution HDATA arising from the data dependent term in the
%	Hessian.
%
%	H = MLPHESS(NET, X, T, HDATA) takes a network data structure NET, a
%	matrix X of input values, and a matrix T of  target values, together
%	with the contribution HDATA arising from the data dependent term in
%	the Hessian, and returns the full Hessian matrix H corresponding to
%	the second derivatives of the negative log posterior distribution.
%	This version saves computation time if HDATA has already been
%	evaluated for the current weight and bias values.
%
%	See also
%	MLP, HESSCHEK, MLPHDOTV, EVIDENCE
%

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

% Check arguments for consistency
errstring = consist(net, 'mlp', x, t);
if ~isempty(errstring);
  error(errstring);
end

if nargin == 4
  % Data term in Hessian needs to be computed
  hdata = datahess(net, x, t, eso_w);
end

[h, hdata] = hbayes(net, hdata);

% Sub-function to compute data part of Hessian
function hdata = datahess(net, x, t, eso_w)

hdata = zeros(net.nwts, net.nwts);

for v = eye(net.nwts);
  hdata(find(v),:) = mlphdotv_weighted(net, x, t, eso_w, v);
end

return


text:
function [dists, neighInds] = nndist(points)
%NNDIST Nearest-neighbor distances of each point in a set
%
%   DISTS = nndist(POINTS)
%   Returns the distance to the nearest neighbor of each point in an array
%   of points.
%   POINTS is an array of points, NP-by-ND.
%   DISTS is a NP-by-1 array containing the distances to the nearest
%   neighbor.
%
%   This functions first computes the Delaunay triangulation of the set of
%   points, then search for nearest distance in the set of each vertex
%   neighbors. This reduces the overall complexity, but difference was
%   noticed only for large sets (>10000 points)
%
%   Example
%     % Display Stienen diagram of a set of random points in unit square
%     pts = rand(100, 2);
%     [dists, inds] = nndist(pts);
%     figure; drawPoint(pts, 'k.');
%     hold on; drawCircle([pts dists/2], 'b');
%     axis equal; axis([-.1 1.1 -.1 1.1]);
%     % also display edges
%     drawEdge([pts pts(inds, :)], 'b');
%
%   See also
%     points2d, distancePoints, minDistancePoints, findPoint
%

% ------
% Author: David Legland
% e-mail: david.legland@nantes.inra.fr
% Created: 2011-12-01,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

% number of points
n = size(points, 1);

% allocate memory
dists = zeros(n, 1);
neighInds = zeros(n, 1);

% in case of few points, use direct computation
if n < 3
    inds = 1:n;
    for i = 1:n
        % compute minimal distance
        [dists(i), indN] = minDistancePoints(points(i,:), points(inds~=i, :));
        if indN >= i
            neighInds(i) = inds(indN) + 1;
        else
            neighInds(i) = inds(indN);
        end
    end
    return;
end

% use Delaunay Triangulation to facilitate computations
DT = delaunayTriangulation(points);

% compute distance to nearest neighbor of each point in the pattern
for i = 1:n
    % find indices of neighbor vertices in Delaunay Triangulation.
    % this set contains the nearest neighbor
    inds = unique(DT.ConnectivityList(sum(DT.ConnectivityList == i, 2) > 0, :));
    inds = inds(inds~=i);

    % compute minimal distance
    [dists(i), indN] = min(distancePoints(points(i,:), points(inds, :)));
    neighInds(i) = inds(indN);
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


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

beamlet = struct(...                      %This is a *2D* structure array with index:  {structNum, beamletNum}.
                    'structureName', '',...    %Structure name
                    'format',[], ...          %Influence storage format e.g., pre-scaled uint8.
                    'influence',[], ...       %Sparse influence/dose element storage, within the structure points: only non-zeros are stored
                    'beamNum',  [], ...       %Beam indentifier.
                    'fullLength',[], ...      %Number of dose elements in non-sparse structure; used to put dose back in.
                    'indexV',[], ...          %Index of nonzero influence/dose values into influence/dose vector
                    'maxInfluenceVal',[], ...  %The maximum influence matrix value for this PB.  Needed to scale dose back up after conversion from uint8 (or whatever) format.
                    'lowDosePoints', [], ...   %Compressed boolean vector indicating points that are scaled from 0...1/256 * maxDose.
                    'sampleRate', 1, ...
                    'strUID','' ...
                    );
return;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function test_issue1601

% WALLTIME 00:10:00
% MEM 2gb
% DEPENDENCY bids_sidecar ft_read_header ft_read_sens

%%

cd(dccnpath('/home/common/matlab/fieldtrip/data/test/issue1601'));

%%

hdr  = ft_read_header('sub-RESP0521/ses-1/ieeg/sub-RESP0521_ses-1_task-Sleep_run-060307_ieeg.vhdr');
dat  = ft_read_data  ('sub-RESP0521/ses-1/ieeg/sub-RESP0521_ses-1_task-Sleep_run-060307_ieeg.vhdr');
evt  = ft_read_event ('sub-RESP0521/ses-1/ieeg/sub-RESP0521_ses-1_task-Sleep_run-060307_ieeg.vhdr');
elec = ft_read_sens  ('sub-RESP0521/ses-1/ieeg/sub-RESP0521_ses-1_task-Sleep_run-060307_ieeg.vhdr');

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
classdef OFA < ALGORITHM
% <single> <real/integer> <large/none> <constrained/none>
% Optimal foraging algorithm

%------------------------------- Reference --------------------------------
% G. Zhu and W. Zhang, Optimal foraging algorithm for global optimization,
% Applied Soft Computing, 2017, 51: 294-313.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        function main(Algorithm,Problem)
            %% Generate random population
            Population = Problem.Initialization();
            [~,rank]   = sort(FitnessSingle(Population));
            Population = Population(rank);

            %% Optimization
            while Algorithm.NotTerminated(Population)
                PopDec     = Population.decs;
                OffDec     = PopDec + Problem.FE./Problem.maxFE.*(rand(size(PopDec))-rand(size(PopDec))).*(PopDec-PopDec([end,floor(unifrnd(ones(1,end-1),2:end))],:));
                Offspring  = Problem.Evaluation(OffDec);
                lambda     = rand(Problem.N,1);
                replace    = lambda.*FitnessSingle(Offspring)./(1+lambda*ceil(Problem.FE/Problem.N)) < FitnessSingle(Population)/ceil(Problem.FE/Problem.N);
                Population(replace) = Offspring(replace);
                [~,rank]   = sort(FitnessSingle(Population));
                Population = Population(rank);
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function handles = setCurrentSimGraphProperties(hObject, handles)
% setCurrentSimGraphProperties Updates the current similarity graph
% information

handles.currentSimGraphType      = handles.SimGraphType;
handles.currentSimGraphNeighbors = get(handles.edtSimGraphNeighbors, 'String');
handles.currentSimGraphEps       = get(handles.edtSimGraphEps, 'String');
handles.currentSimGraphSigma     = get(handles.edtSimGraphSigma, 'String');

guidata(hObject, handles);

author: beckel

repo: nilm-eval

save_path: github-repos/MATLAB/beckel-nilm-eval


text:
function hf_out = lhs_operation_gpu(hf, samplesf, reg_filter, sample_weights)

% This is the left-hand-side operation in Conjugate Gradient

% Get sizes
num_features = length(hf);
filter_sz = zeros(num_features,2);
for k = 1:num_features
    filter_sz(k,:) = [size(hf{k},1), size(hf{k},2)];
end
[~, k1] = max(filter_sz(:,1));  % Index for the feature block with the largest spatial size
block_inds = 1:num_features;
block_inds(k1) = [];
output_sz = [size(hf{k1},1), 2*size(hf{k1},2)-1];

% Compute the operation corresponding to the data term in the optimization
% (blockwise matrix multiplications)
%implements: A' diag(sample_weights) A f

% sum over all features and feature blocks
sh = sum(bsxfun(@times, samplesf{k1}, hf{k1}), 3);    % assumes the feature with the highest resolution is first
pad_sz = cell(1,1,num_features);
for k = block_inds
    pad_sz{k} = (output_sz - [size(hf{k},1), 2*size(hf{k},2)-1]) / 2;

    sh(1+pad_sz{k}(1):end-pad_sz{k}(1), 1+pad_sz{k}(2):end,1,:) = ...
        sh(1+pad_sz{k}(1):end-pad_sz{k}(1), 1+pad_sz{k}(2):end,1,:) + sum(bsxfun(@times, samplesf{k}, hf{k}), 3);
end

% weight all the samples and take conjugate
sh = conj(bsxfun(@times,sample_weights,sh));

% multiply with the transpose
hf_out = cell(1,1,num_features);
hf_out{k1} = conj(sum(bsxfun(@times, sh, samplesf{k1}), 4));
for k = block_inds
    hf_out{k} = conj(sum(bsxfun(@times, sh(1+pad_sz{k}(1):end-pad_sz{k}(1), 1+pad_sz{k}(2):end,1,:), samplesf{k}), 4));
end

% compute the operation corresponding to the regularization term (convolve
% each feature dimension with the DFT of w, and the tramsposed operation)
% add the regularization part
% hf_conv = cell(1,1,num_features);
for k = 1:num_features
    reg_pad = min(size(reg_filter{k},2)-1, size(hf{k},2)-1);

    % add part needed for convolution
    hf_conv = cat(2, hf{k}, conj(rot90(hf{k}(:, end-reg_pad:end-1, :), 2)));

    % do first convolution
    hf_conv = convn(hf_conv, reg_filter{k});

    % do final convolution and put toghether result
    hf_out{k} = hf_out{k} + convn(hf_conv(:,1:end-reg_pad,:), reg_filter{k}, 'valid');
end

end

author: martin-danelljan

repo: ECO