text:
function sw = BLanalyticalCoreyFoam(mug, muw)
%BLANALYTICAL Analytical solution of Buckley-Leverett equation
% Written by Ali A. Eftekhari
% See the license file
%muw = 10e-3;
%muo = 10e-3;
u = 1e-6/(pi()*0.038^2/4)/60;
phi = 0.2;
k = 2e-12; %[m^2]
sw_end = .5;
fmmob = 1.01889e+05;
fmdry = 1.20000e-01;
epdry = 5.0000e+03;
swc = 0.07;
sw_in = swc;
sgr = 0.0;
krg0 = 1;
ng = 2;
krw0 = 1;
nw = 2;
sws = @(sw)((sw>swc).*(sw-swc)/(1-sgr-swc));
kr = @(sw)(krg0*(1-sws(sw)).^ng);
% fm = @(sw)((sw>fmdry).*(1+fmmob*(0.5+atan(epdry.*(sw-fmdry))/pi()))+(sw<=fmdry));
fm = @(sw)(1+fmmob*(0.5+atan(epdry.*(sw-fmdry))/pi()));
krg = @(sw)(kr(sw)./fm(sw));
krw = @(sw)(krw0*sws(sw).^nw);
dkrwdsw = @(sw)(nw*krw0*(1/(1-sgr-swc))*sws(sw).^(nw-1));
dkrdsw = @(sw)((-krg0*ng*(1-sws(sw)).^(ng-1))/(-swc-sgr+1));
% fm = @(sw)((sw>fmdry).*(1+fmmob*(0.5+atan(epdry.*(sw-fmdry))/pi()))+(sw<=fmdry));
dfmdsw = @(sw)(((epdry*fmmob)./(pi*(epdry^2*(sw-fmdry).^2+1))));
% dfmdsw = @(sw)((epdry*fmmob)./(pi*(epdry^2*(sw-fmdry).^2+1)));
dkrgdsw = @(sw)((dkrdsw(sw).*fm(sw)-dfmdsw(sw).*kr(sw))./fm(sw).^2);
fw = @(sw)((krw(sw)/muw)./(krw(sw)/muw+krg(sw)/mug));
dfwdsw = @(sw)((dkrwdsw(sw)/muw.*(krw(sw)/muw+krg(sw)/mug)- ...
    (dkrwdsw(sw)/muw+dkrgdsw(sw)/mug).*krw(sw)/muw)./ ...
    (krg(sw)/mug+krw(sw)/muw).^2);
s = linspace(swc,1,1000);
figure(1)
subplot(2,2,1);
plot(s, krw(s), s, kr(s), s, krg(s));
xlabel('S_w'); ylabel('rel perms'); legend('liq', 'gas', 'foam');
F = @(sw)(dfwdsw(sw)-(fw(sw_end)-fw(sw))/(sw_end-sw));
sw_shock = fzero(F, [swc+eps,0.3]);
subplot(2,2,2);
plot(s, fw(s), [sw_shock sw_end], [fw(sw_shock) fw(sw_end)]);
xlabel('S_w'); ylabel('f_w');
% plot(s, fw(s), [sw_end sw_shock], [fw(sw_end) fw(sw_shock)]);
s1 = linspace(sw_in, sw_shock, 50);
xt_s1 = u/phi*dfwdsw(s1);
xt_s = u/phi*dfwdsw(s);
xt_shock = u/phi*dfwdsw(sw_shock);
subplot(2,2,3);
plot(xt_s, s, '--', ...
    [xt_s1 xt_shock xt_shock max(xt_s)], [s1 sw_shock sw_end sw_end])
xlabel('x/t [m/s]'); ylabel('S_w');
subplot(2,2,4);
plot([xt_s1 xt_shock xt_shock 10*xt_shock], [s1 sw_shock sw_end sw_end])
xlabel('x/t [m/s]'); ylabel('S_w');
sw_prf = [s1 sw_shock sw_end sw_end];
xt_prf = [xt_s1 xt_shock xt_shock max(xt_s)];
sw = @(xt)([interp1(xt_s1, s1, xt(xt<xt_shock)) sw_end*xt(xt>=xt_shock)]);
% Now I can use the data to calculate the pressure drop across my domain
L = 0.17; % [m]
t = eps:20:5*40*60; % [s] time
dp = zeros(1,length(t));
for i =1:length(t)
    x = 0:0.001:L;
    xt = x/t(i);
    dp(i) = trapz(x, u./(k*(krg(sw(xt))/mug+krw(sw(xt))/muw)));

end
figure(2);plot(t/60, dp/1e5, '.')
xlabel('time [min]'); ylabel('pressure drop [bar]');
end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function test_old_ft_freqanalysis

% MEM 1gb
% WALLTIME 00:10:00

function test_ft_freqanalysis
% DEPENDENCY_FT_FREQANALYSIS
% This script tests the ft_freqanalysis functions using simulated data
% A. Stolk and J.M. Schoffelen

% simulate one second of data, samplefreq = 1200 hz
t = (1:1200)/1200;
a = cos(2*pi*10*t);
b = sin(2*pi*25*t);
c = a + b;

% simulate preprocessed data
%cfg           = [];
%cfg.layout    = 'CTF275.lay';
%cfg.layout    = prepare_layout(cfg);
%data.label    =  cfg.layout.label(1:273,1);
% data.grad.pnt = zeros(595,3);
% data.grad.ori = zeros(595,3);
% data.grad.tra = zeros(302,595);
data.fsample  = 1200;
data.label    = {'chan01';'cos10';'sin25'};
%for j = 1:273
for j = 1
    %data.trial{1,1}(j,:) = c;
  data.trial{1,1} = [c;a;b];
  data.time{1,1}(j,:)  = t;
end

% ft_freqnalysis_mtmfft
cfg              = [];
cfg.output       = 'pow';
%cfg.channel      = 'MEG';
cfg.channel      = 'all';
cfg.method       = 'mtmfft';
cfg.taper        = 'hanning';
cfg.foilim       = [1 40];
cfg.keeptrials   = 'no';
cfg.keeptapers   = 'no';
mtmfft           = ft_freqanalysis(cfg, data);

% check whether the powerpeaks are at the given frequencies
if mtmfft.powspctrm(1,9) < mtmfft.powspctrm(1,10) > mtmfft.powspctrm(1,11) && ...
        mtmfft.powspctrm(1,24) < mtmfft.powspctrm(1,25) > mtmfft.powspctrm(1,26);
else
    error('test_ft_freqanalysis:notEqual', 'Incorrect output for ft_freqanalysis_mtmfft.');
end

% ft_freqnalysis_mtmconvol
cfg              = [];
cfg.output       = 'pow';
%cfg.channel      = 'MEG';
cfg.channel      = 'all';
cfg.method       = 'mtmconvol';
cfg.taper        = 'hanning';
cfg.keeptrials   = 'no';
cfg.keeptapers   = 'no';
cfg.foi          = 1:1:40;
cfg.t_ftimwin    = 0.5 * ones(1,length(cfg.foi)); % 500 ms
cfg.toi          = 0.3:0.05:0.75; % center 1/2 second
mtmconvol        = ft_freqanalysis(cfg, data);

% check whether the powerpeaks are at the given frequencies
if mtmconvol.powspctrm(1,9,1) < mtmconvol.powspctrm(1,10,1) > mtmconvol.powspctrm(1,11,1) && ...
        mtmconvol.powspctrm(1,24,1) < mtmconvol.powspctrm(1,25,1) > mtmconvol.powspctrm(1,26,1);
else
    error('test_ft_freqanalysis:notEqual', 'Incorrect output for ft_freqanalysis_convol.');
end

% test new implementation specest

% ft_freqnalysis_mtmfft
cfg              = [];
cfg.output       = 'fourier';
cfg.channel      = 'all';
cfg.method       = 'mtmfft';
cfg.taper        = 'hanning';
cfg.foilim       = [1 40];
cfg.keeptrials   = 'yes';
cfg.keeptapers   = 'yes';
mtmfft1          = ft_freqanalysis(cfg, data);
mtmfft2          = ft_freqanalysis(cfg, data, 1);

x = [mtmfft1.fourierspctrm(:,1,10) mtmfft2.fourierspctrm(:,1,10)];
y = [mtmfft1.fourierspctrm(:,1,25) mtmfft2.fourierspctrm(:,1,25)];

% x(2) should only have real component (angle = 0)
% y(2) should only have imag component (negative value, angle = -pi/2)
% abs(x(1))==abs(x(2))
% abs(y(1))==abs(y(2))

% ft_freqnalysis_mtmconvol
cfg              = [];
cfg.output       = 'fourier';
cfg.channel      = 'all';
cfg.method       = 'mtmconvol';
cfg.taper        = 'hanning';
cfg.keeptrials   = 'no';
cfg.keeptapers   = 'no';
cfg.foi          = 1:1:40;
cfg.t_ftimwin    = 0.5 * ones(1,length(cfg.foi)); % 500 ms
cfg.toi          = 0.25:1./1200:0.75; % center 1/2 second
mtmconvol1        = ft_freqanalysis(cfg, data);
mtmconvol2        = ft_freqanalysis(cfg, data, 1);

x = [squeeze(mtmconvol1.fourierspctrm(:,2,10,:)) ...
     squeeze(mtmconvol2.fourierspctrm(:,2,10,:))];
y = [squeeze(mtmconvol1.fourierspctrm(:,3,25,:)) ...
     squeeze(mtmconvol2.fourierspctrm(:,3,25,:))];

% observations:
%  Old implementation has 1 Nan at the beginning
%  New implementation has 1 Nan at the end

% expectations:
%  angle(x(1,2)) = -pi, (phase of cosine @10Hz @0.25 s: this is approximately true, but not exact
%  angle(y(1,2)) = 0 (phase of sine @25Hz @0.25 s = 6.25 cycle: this is approximately true, but not exact

% FIXME should we look into this?
% issues are probably related to even numbered t_ftimwins...

% ft_freqnalysis_mtmconvol
cfg              = [];
cfg.output       = 'fourier';
cfg.channel      = 'all';
cfg.method       = 'mtmconvol';
cfg.taper        = 'dpss';
cfg.keeptrials   = 'no';
cfg.keeptapers   = 'no';
cfg.foi          = 1:1:40;
cfg.t_ftimwin    = 0.5 * ones(1,length(cfg.foi)); % 500 ms
cfg.toi          = 0.25:1./1200:0.75; % center 1/2 second
cfg.tapsmofrq    = ones(1,numel(cfg.foi)).*4;
mtmconvol1        = ft_freqanalysis(cfg, data);
mtmconvol2        = ft_freqanalysis(cfg, data, 1);

x = [squeeze(mtmconvol1.fourierspctrm(:,2,10,:)); ...
     squeeze(mtmconvol2.fourierspctrm(:,2,10,:))];
y = [squeeze(mtmconvol1.fourierspctrm(:,3,25,:)); ...
     squeeze(mtmconvol2.fourierspctrm(:,3,25,:))];

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [ n_data, n, area ] = hypersphere_01_area_values ( n_data )

%*****************************************************************************80
%
%% HYPERSPHERE_01_AREA_VALUES returns some areas of the unit sphere in ND.
%
%  Discussion:
%
%    The formula for the surface area of the unit sphere in N dimensions is:
%
%      Sphere_Unit_Area ( N ) = 2 * PI^(N/2) / Gamma ( N / 2 )
%
%    Some values of the function include:
%
%       N   Area
%
%       2    2        * PI
%       3  ( 4 /    ) * PI
%       4  ( 2 /   1) * PI^2
%       5  ( 8 /   3) * PI^2
%       6  ( 1 /   1) * PI^3
%       7  (16 /  15) * PI^3
%       8  ( 1 /   3) * PI^4
%       9  (32 / 105) * PI^4
%      10  ( 1 /  12) * PI^5
%
%    For the unit sphere, Area(N) = N * Volume(N)
%
%    In Mathematica, the function can be evaluated by:
%
%      2 * Pi^(n/2) / Gamma[n/2]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.
%    On input, if N_DATA is 0, the first test data is returned, and
%    N_DATA is set to the index of the test data.  On each subsequent
%    call, N_DATA is incremented and that test data is returned.  When
%    there is no more test data, N_DATA is set to 0.
%
%    Output, integer N, the spatial dimension.
%
%    Output, real AREA, the area of the unit sphere
%    in that dimension.
%
  n_max = 20;

  area_vec = [ ...
     0.2000000000000000E+01, ...
     0.6283185307179586E+01, ...
     0.1256637061435917E+02, ...
     0.1973920880217872E+02, ...
     0.2631894506957162E+02, ...
     0.3100627668029982E+02, ...
     0.3307336179231981E+02, ...
     0.3246969701133415E+02, ...
     0.2968658012464836E+02, ...
     0.2550164039877345E+02, ...
     0.2072514267328890E+02, ...
     0.1602315322625507E+02, ...
     0.1183817381218268E+02, ...
     0.8389703410491089E+01, ...
     0.5721649212349567E+01, ...
     0.3765290085742291E+01, ...
     0.2396678817591364E+01, ...
     0.1478625959000308E+01, ...
     0.8858104195716824, ...
     0.5161378278002812 ];

  n_vec = [ ...
     1, ...
     2, ...
     3, ...
     4, ...
     5, ...
     6, ...
     7, ...
     8, ...
     9, ...
    10, ...
    11, ...
    12, ...
    13, ...
    14, ...
    15, ...
    16, ...
    17, ...
    18, ...
    19, ...
    20 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    n = 0;
    area = 0.0;
  else
    n = n_vec(n_data);
    area = area_vec(n_data);
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function eeg_write_brainstorm(fileprefix,data)

% eeg_write_brainstorm - Write EEG data into brainstorm format
%
% Useage: eeg_write_brainstorm(fileprefix, data)
%
% file  - the path + filename, eg "c:\subj_data.mat"
%         It is actually a fileprefix, as the _data.mat are
%         appended by this function (they are mandatory)
%
% data  - a matlab structure with brainstorm data fields.
%         There are 2 essential fields (others are initialised):
%
%         F       a matrix of voltage values (in Volts) with
%                 electrodes in rows and data points in columns.
%         Time    Time is a row vector with time points (in sec)
%                 for each column of the voltage data.
%
% This script does not convert the voltage data from uV to Volts
% or the Time from msec to sec.  Also, it assumes that the associated
% Channel struct for this EEG data has the last Channel.Type = 'EEG REF'
% so the ChannelFlag array output from this function contains a -1 at
% index = size(data.F,1) + 1.
%
% Comment:  See the brainstorm website for more info, at
%           http://neuroimage.usc.edu/.
%

% $Revision: 1.1 $ $Date: 2009-04-28 22:13:54 $

% Licence:  GNU GPL, no express or implied warranties
% Author:   07/2001, Darren.Weber_at_radiology.ucsf.edu
%           10/2002, Darren.Weber_at_radiology.ucsf.edu
%                    modified warnings and channelflag
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fprintf('EEG_WRITE_BRAINSTORM...\n'); tic;

if isfield(data,'F'),
    F = data.F;
else
    msg = strcat('data structure must contain an ''F'' field.\n',...
                 'F is a matrix of voltage values (in Volts)\n',...
                 'with electrodes in rows and data points in columns.\n');
    error(sprintf(msg));
end
if isfield(data,'Time'),
    Time = data.Time;
else
    msg = strcat('data structure must contain a ''Time'' field.\n',...
                 'Time is a row vector with time points (in sec)\n',...
                 'for each column of the voltage data.\n');
    error(sprintf(msg));
end
if isfield(data,'ChannelFlag'),
    ChannelFlag = data.ChannelFlag;
else
%     msg = strcat('data structure can contain a ''ChannelFlag'' field.\n',...
%                  'ChannelFlag is a row vector indicating the status\n',...
%                  'of the electrodes (-1 is dead, 0 is ignored, 1 is good).\n',...
%                  'Creating this vector = ones(1,size(data.F,1)).\n');
%     fprintf(msg);

    ChannelFlag = [ones(1,size(F,1)), -1];
    fprintf('...assume data Nchannels + 1 is ref, so ChannelFlag(Nchannels + 1) = -1\n');
end
if isfield(data,'NoiseCov'),
    NoiseCov = data.NoiseCov;
else
%     msg = strcat('data structure can contain a ''NoiseCov'' field.\n',...
%                  'This is a square matrix the size of the electrodes\n',...
%                  'Creating this vector = ones(size(data.F,1)).\n');
%     fprintf(msg);
    NoiseCov = [];
end
if isfield(data,'SourceCov'),
    SourceCov = data.SourceCov;
else
%     msg = strcat('data structure can contain a ''SourceCov'' field.\n',...
%                  'This is a square matrix the size of ImageGridLoc(?)\n',...
%                  'Creating this matrix = [].\n');
%     fprintf(msg);
    SourceCov = [];
end
if isfield(data,'Project'),
    Projector = data.Project;
else
%     msg = strcat('data structure can contain a ''Project'' field.\n',...
%                  'Not sure what it is, but it is created for the ascii\n',...
%                  'tutorial data.  Creating this matrix = [].\n');
%     fprintf(msg);
    Project = [];
end
if isfield(data,'Projector'),
    Projector = data.Projector;
else
%     msg = strcat('data structure can contain a ''Projector'' field.\n',...
%                  'This is a matrix of size length of electrodes by rank,\n',...
%                  'not necessarily orthogonal due to Good-Channel selections.\n',...
%                  'Orthogonalized as ''U'', then the data are to be projected away\n',...
%                  'from the projector as F'' = F - U*(U^t * F).\n',...
%                  'Creating this matrix = [].\n');
%     fprintf(msg);
    Projector = [];
end
if isfield(data,'Comment'),
    Comment = data.Comment;
else
%     msg = strcat('data structure can contain a ''Comment'' field.\n',...
%                  'This is a character string describing the data.\n',...
%                  'Creating this field = file.\n');
%     fprintf(msg);
    Comment = fileprefix;
end
if isfield(data,'Device'),
    Device = data.Device;
else
%     msg = strcat('data structure can contain a ''Device'' field.\n',...
%                  'This is a character string describing the device\n',...
%                  'used to acquire the data (eg, SYNAMPS)\n',...
%                  'Creating an empty field.\n');
%     fprintf(msg);
    Device = 'Neuromag_Planar_122'; % the only supported device!
end

[path,file,ext] = fileparts(fileprefix);
ext = '_data.mat';

file = fullfile(path,[file ext]);

save(file,'F','Time','ChannelFlag','NoiseCov','SourceCov',...
          'Project','Projector','Comment','Device');

fprintf('...wrote brainstorm data to:\n\t%s\n', file);

t = toc; fprintf('...done (%6.2f sec).\n\n',t);

return

author: PatternRecognition

repo: OpenBMI


text:
function plotActivityTrace( inferred, filtered, raw, plotControl)

% This code plots the activity trace for all the ROIs
% Author: Weijian Yang

% input parameters:
% inferred, filtered, raw (from the output of "signalExtraction.m")
% plotControl.frameRate: frame rate in fps
% plotControl.normalization: normalize the amplitude of the temporal trace for plotting purpose
% plotControl.sep: vertical (amplitude) offset between adjacent temporal traces
% plotControl.displayLabel: whether to display ROI ID for the temporal traces
% plotControl.plotInferred: whether the inferred signal is plotted
% plotControl.plotFiltered: whether the filtered signal is plotted
% plotControl.plotRaw: whether the raw signal is plotted
% plotControl.rollingView: whether the display is set such that only 20 temporal traces are displayed at a time, followed by any key to move onto the next 20 traces, and etc.

if isfield(plotControl,'frameRate') frameRate=plotControl.frameRate;
else frameRate=1; end

if isfield(plotControl,'normalization') normalization=plotControl.normalization;
else normalization=1; end

if isfield(plotControl,'sep') sep=plotControl.sep;
else sep=1.2; end

if isfield(plotControl,'displayLabel') displayLabel=plotControl.displayLabel;
else displayLabel=1; end

if isfield(plotControl,'plotInferred') plotInferred=plotControl.plotInferred;
else plotInferred=1; end

if isfield(plotControl,'plotFiltered') plotFiltered=plotControl.plotFiltered;
else plotFiltered=1; end

if isfield(plotControl,'plotRaw') plotRaw=plotControl.plotRaw;
else plotRaw=0; end

if isfield(plotControl,'rollingView') rollingView=plotControl.rollingView;
else rollingView=0; end

t=1/frameRate*(1:size(inferred,2));
K=size(inferred,1);

figure;
hold on;

if normalization==1
    for idx=1:K
        if plotRaw  plot(t,raw(idx,:)/max(raw(idx,:))-idx*sep,'color',[0.8 0.8 0.8],'linewidth',2); end
        if plotFiltered  plot(t,filtered(idx,:)/max(filtered(idx,:))-idx*sep,'color',[1 0.7 0.7],'linewidth',1.5); end
        if plotInferred  plot(t,inferred(idx,:)/max(inferred(idx,:))-idx*sep,'color',[0 0 1],'linewidth',1); end
    end
    ylabel('Normalized \DeltaF/F');
else
    for idx=1:K
        if plotRaw  plot(t,raw(idx,:)-idx*sep,'color',[0.8 0.8 0.8],'linewidth',2); end
        if plotFiltered  plot(t,filtered(idx,:)-idx*sep,'color',[1 0.7 0.7],'linewidth',1.5); end
        if plotInferred  plot(t,inferred(idx,:)-idx*sep,'color',[0 0 1],'linewidth',1); end
    end
    ylabel('\DeltaF/F');
end
xlabel('Time (s)');

if displayLabel
    for idx=1:K
        text(max(t)*1.02,-idx*sep, num2str(idx));
    end
end

h=gca;
if ~rollingView
    axis(h,[0 max(t) -sep*(K+1) sep]);
else
    idx=1;
    displayNum=20;
    while idx+displayNum<K
        axis(h,[0 max(t) -sep*(idx+displayNum) sep*(-idx+1)]);
        idx=idx+displayNum;
        input('Press ENTER to continue...');
    end
    axis(h,[0 max(t) -sep*(idx+displayNum) sep*(-idx+1)]);
    input('Press ENTER to continue...');
    axis(h,[0 max(t) -sep*(K+1) sep]);
end

author: flatironinstitute

repo: CaImAn-MATLAB


text: %EQ Mapping overload

author: marianux

repo: ecg-kit

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


text:
function d = dot_outer(S3G,q,varargin)
% return outer inner product of all nodes within a eps neighborhood
%
% Syntax
%   d = dot_outer(SO3G,nodes,'epsilon',radius)
%
% Input
%  SO3G   - @SO3Grid
%  nodes  - @quaternion
%  radius - double
%
% Output
%  d      - sparse matrix
%

if ~isa(S3G,'SO3Grid')
  d = dot_outer(q,S3G,varargin{:}).';
  return
end

epsilon = get_option(varargin,'epsilon',pi);

if check_option(varargin,{'full','all'})
  if check_option(varargin,'nospecimensymmetry')
    S3G.SS = specimenSymmetry;
  end
  d = dot_outer@orientation(S3G,q,varargin{:});

else

  d = sparse(length(S3G),length(q));

  % rotate q according to SO3Grid.center
  if ~isempty(S3G.center),q = inv(S3G.center) * q; end %#ok<MINV>

  % extract SO3Grid
  [ybeta,yalpha,ialphabeta,palpha] = getdata(S3G.alphabeta);

  ygamma = [S3G.gamma.points];
  sgamma = [S3G.gamma.min];
  pgamma = S3G.gamma(1).period;
  igamma = cumsum([0,GridLength(S3G.gamma)]);

  % correct for specimen symmetry
  if check_option(varargin,'nospecimensymmetry')
    qss = quaternion.id;
    palpha = 2*pi;
  else
    qss = quaternion(rotation_special(S3G.SS));
    palpha = max(palpha,pi);
  end

  % for finding the minimial beta angle
  qcs = quaternion(S3G.CS.rotation_special);

  [xalpha,xbeta,xgamma] = Euler( qss * quaternion(q) * qcs ,'ZYZ');

  ncs = length(qss)*length(qcs);
  cs = 0:length(q):ncs*length(q);

  for k=1:ncs

    ndx = cs(k)+1:cs(k+1);

    dist = SO3Grid_dist_region(yalpha,ybeta,ygamma, ...
      sgamma, int32(igamma), int32(ialphabeta), palpha, pgamma, ...
      xalpha(ndx), xbeta(ndx), xgamma(ndx), epsilon);

    if nnz(dist) > 0, d = max(d,dist); end

  end

end

author: mtex-toolbox

repo: mtex

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


text:
% function plot_srErrorOrigin
% %
% %
% %
% % Author: van Stiphout, Thomas, vanstiphout@sed.ethz.ch

fRes=10000;
vEdges=[-9.9:0.2:9.9];

% Plot-list
sPlotList=['mCompose.m05';...
           'mCompose.m06';...
           'mCompose.m07';...
           'mCompose.m08';...
%            'mCompose.m12';...
           'mCompose.m13'];

mColorLine=lines(size(sPlotList,1));

figure;
for i=1:size(sPlotList,1)
    if i==1
       subplot(2,1,1)
    else
        hold on;
    end
    eval(sprintf('mPlotit=%s.mResult1;',sPlotList(i,:)));
    vVal=reshape(squeeze(mPlotit(:,1,:)),...
    size(mPlotit,1)*size(mPlotit,3),1);
vRes=reshape(squeeze(mPlotit(:,2,:)),...
    size(mPlotit,1)*size(mPlotit,3),1);
% vSel=(vRes<fRes);
N=histc(vVal,vEdges);
plot(vEdges,N,'Color',mColorLine(i,:),'LineWidth',2);
clear X N mPlotit vVal vRes
end
legend(sPlotList,'location','NW');

% % probZ
% for i=1:size(sPlotList,1)
%     if i==1
%                 subplot(2,2,2)
%     else
%         hold on;
%     end
%     eval(sprintf('mPlotit=%s.mResult2;',sPlotList(i,:)));
%     vVal=reshape(squeeze(mPlotit(:,1,:)),...
%     size(mPlotit,1)*size(mPlotit,3),1);
% vRes=reshape(squeeze(mPlotit(:,2,:)),...
%     size(mPlotit,1)*size(mPlotit,3),1);
% vSel=(vRes<fRes);
% [N,X]=histogram(calc_ProbColorbar2Value(vVal(vSel)),20);
% plot(X,N,'Color',mColorLine(i,:),'LineWidth',2);
% clear X N;
% end

% % cdf plots Z
% for i=1:size(sPlotList,1)
%     if i==1
%                 subplot(2,2,3)
%     else
%         hold on;
%     end
%     eval(sprintf('mPlotit=%s.mResult1;',sPlotList(i,:)));
%     vVal=reshape(squeeze(mPlotit(:,1,:)),...
%     size(mPlotit,1)*size(mPlotit,3),1);
% vRes=reshape(squeeze(mPlotit(:,2,:)),...
%     size(mPlotit,1)*size(mPlotit,3),1);
% vSel=(vRes<fRes);
% [h,stats(i)] = cdfplot(vVal(vSel));
% set(h,...
%     'Color',mColorLine(i,:),'LineWidth',2);
% end
% sLegend=[];
% for i=1:size(sPlotList,1)
%     sTmp=sprintf('mean=%s,std=%s',...
%         num2str(stats(i).mean,'%6.4f'),num2str(stats(i).std));
%     sLegend=[sLegend;sTmp]
% end
% legend(sLegend,'location','NW');

% cdf prob Z
for i=1:size(sPlotList,1)
    if i==1
       subplot(2,1,2)
    else
        hold on;
    end
    eval(sprintf('mPlotit=%s.mResult2;',sPlotList(i,:)));
    vVal=reshape(squeeze(mPlotit(:,1,:)),...
    size(mPlotit,1)*size(mPlotit,3),1);
vRes=reshape(squeeze(mPlotit(:,2,:)),...
    size(mPlotit,1)*size(mPlotit,3),1);
% vSel=(vRes<fRes);
[h,stats(i)] = cdfplot(calc_ProbColorbar2Value(vVal));
set(h,...
    'Color',mColorLine(i,:),'LineWidth',2);
end

% sLegend=[];
% for i=1:size(sPlotList,1)
%     sTmp=sprintf('mean=%s,std=%smin=%s,max=%s',...
%         num2str(stats(i).mean,'%010.5f'),num2str(stats(i).std,'%010.5f'),...
%         num2str(stats(i).min,'%010.5f'), num2str(stats(i).max,'%010.5f'));
%     sLegend=[sLegend;sTmp];
% end
% legend(sLegend,'location','NW');

author: CelsoReyes


text:
function w_nulldiff = sphNullformer_diff(order, src_dirs)
%SPHNULLFORMER_DIFF Beamweights for a beamformer with nulls at specified directions,
%and an omnidirectional constraint
%
%   For a set of K directions, computes the beamforming weights that places
%   nulls at them while trying to keep the overall response as
%   omnidirectional as possible. Useful when trying to extract
%   reverberant/diffuse sound while suppressing strong directional sources.
%   From a statistical formulation, it is similar to the beamformer
%   proposed in
%
%       Thiergart, O. and Habets, E.A., 2014.
%       Extracting reverberant sound using a linearly constrained minimum variance spatial filter.
%       IEEE Signal Processing Letters, 21(5), pp.630-634.
%
%   but in the SHD instead of the microphone signal domain.
%
%   Inputs:
%       order:  order of SH signals
%       src_dirs:  Kx2 [azi elev] directions to place the nulls
%
%   Outputs:
%       w_nulldiff:  (order+1)^2xK matrix of beamweights, one set per
%           column for each beamforming direction
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% SPHNULLFORMER_DIFF.M - 5/10/2016
% Archontis Politis, archontis.politis@aalto.fi
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    nSH = (order+1)^2;
    nSrc = size(src_dirs,1);
    src_dirs2 = [src_dirs(:,1) pi/2-src_dirs(:,2)]; % convert from azi-elev to azi-incl

    % steering vectors in the SHD
    Y_nulldiff = getSH(order, src_dirs2, 'real');

    % compute weight vector
    c = zeros(nSrc+1,1);
    c(1) = 1; % constraint vector
    g = zeros(nSH,1);
    g(1) = 1;
    A = [g Y_nulldiff'];
    w_nulldiff = pinv(A')*c;
end

author: polarch

repo: Spherical-Array-Processing


text:
function [models] = featureSelectionTime_HN(X,Y,maxOrder,nBoot,Info,censoring,seed)
% -------------------------------------------------------------------------
% function [models] = featureSelection(X,Y,maxOrder,nBoot,Info,imbalance,batchNum)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function computes feature set selection according to the 0.632+
% bootstrap methodology for an input matrix of features and and input
% outcome vector, and for multiple model orders as defined by the user.
% See ref. [1] for more details. This function uses logistic regression
% utilities from DREES <http://www.cerr.info/drees>, and a fast
% implementation of AUC calculation by Enric Junqué de Fortuny that is
% available at: <http://www.mathworks.com/matlabcentral/fileexchange/41258-faster-roc-auc>
%
% NOTE: This function now maximizes 0.5*AUC + 0.5*(1-abs(SENSITIVITY-SPECIFICITY))
% -------------------------------------------------------------------------
% REFERENCE:
% [1] Vallieres, M. et al. (2015). A radiomics model from joint FDG-PET and
%     MRI texture features for the prediction of lung metastases in soft-tissue
%     sarcomas of the extremities. Physics in Medicine and Biology, 60(14),
%     5471-5496. doi:10.1088/0031-9155/60/14/5471
% -------------------------------------------------------------------------
% INPUTS:
% - X: Matrix of size [nInst X nFeat], specifying the numerical data of the
%      features of the input features, where 'nInst' refers to the number
%      of instances in X, and 'nFeat' to the number of features in X.
%      Each column is a different feature.
% - Y: Column vector of size [nInst X 1] specifying the outcome status
%      (1 or 0) for all instances.
% - maxOrder: Integer specifying the maximal model order to construct.
% - nBoot: Number of bootstrap samples to use.
% - Info: Cell of size [nFeat X 1] of strings specifying the name of each
%         feature in 'X'.
% - imbalance: String specifying the type of imbalance-adjustement strategy
%              employed. Either 'IABR' for imbalance-adjusted bootstrap
%              resampling (see ref.[1]), or 'IALR' for imbalance-adjusted
%              logistic regression.
% - batchNum: (optional input). If present, integer that specifies the
%             batch number for parallelization purposes.
% -------------------------------------------------------------------------
% OUTPUTS:
% - models: Structure specifying the resulting multivariable models for the
%           input feature set in 'data', for each model order. Example for
%           order 4:
%        --> models.Order4.Data: Matrix of size [nInst X order], specifying
%                                the selected features, in order of
%                                selection.
%        --> models.Order4.Name: Cell specifying the names of the selected
%                                features, in order of selection.
% -------------------------------------------------------------------------
% AUTHOR(S):
% - Martin Vallieres <mart.vallieres@gmail.com>
% - DREES development team <http://www.cerr.info/drees> (logistic regression)
% - Enric Junqué de Fortuny (fastAUC.cpp)
% -------------------------------------------------------------------------
% HISTORY:
% - Creation - May 2015
% - Revision I - July 2015: including imbalance-adjusted logistic regression
% - Revision II - July 2015: maximizing 0.5*AUC + 0.5*(1-abs(SENSITIVITY-SPECIFICITY))
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015  Martin Vallieres
% --> Copyright 2010, Joseph O. Deasy, on behalf of the DREES development team.
%
%    This package is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This package is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this package.  If not, see <http://www.gnu.org/licenses/>.
%
%    _______________________________________________________________
%
% --> Copyright (c) 2013, Enric Junqué de Fortuny
%     All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.
% -------------------------------------------------------------------------

% RANDOM NUMBER GENERATOR SEED
rng(seed);

% INITIALIZATION
top = 1-1/exp(1);
low = 1/exp(1);
nFeat = size(X,2);
models = struct; % Final model structure

% GETTING BOOTSTRAP SAMPLES FOR ALL EXPERIMENTS
nInst = numel(Y);
trainSets = ceil(nInst .* rand(nInst,nBoot));
testSets = findBootTestSet(trainSets);

% FORWARD FEATURE SELECTION (for all different starters)
modelMat = zeros(nFeat,maxOrder);
metricMat = zeros(nFeat,maxOrder);
for i = 1:nFeat
     indLeft = 1:nFeat;
     % Order 1
     indLeft(i) = [];
     modelMat(i,1) = i;
     Xtrain = X(:,i);
     [Xtrain] = normalizeZeroOne(Xtrain);

     % Cox regression code
     CIdata = applyCoxfit(Xtrain,Y,Xtrain,Y,censoring,censoring);

     CItemp = 0;
     for n = 1:nBoot
         Xtrain = X(trainSets(:,n),i);Xtest = X(testSets{n},i); Ytrain = Y(trainSets(:,n),1); Ytest = Y(testSets{n},1);
         [Xtrain,Xtest] = normalizeZeroOne(Xtrain,Xtest);

         % Cox regression code
         CIboot = applyCoxfit(Xtrain,Ytrain,Xtest,Ytest,censoring(trainSets(:,n)),censoring(testSets{n}));

         % FOR CI (concordance index)
         alpha = top/(1-low*(CIdata-CIboot)/(CIdata-0.5+eps));
         if alpha > 1
             alpha = 1;
         elseif alpha < top
             alpha = top;
         end
         if CIboot < 0.5
             CIboot = 0.5;
         end
         CItemp = CItemp + (1-alpha)*CIdata + alpha*CIboot;

     end
     CItemp = CItemp/nBoot;
     metricMat(i,1) = CItemp;

     % Going for orders 2 to maxOrder
     for j = 2:maxOrder
         maxMetric = 0;
         for k = 1:(nFeat-j+1)
             indexModel = [modelMat(i,1:(j-1)),indLeft(k)];
             Xtrain = X(:,indexModel);
             [Xtrain] = normalizeZeroOne(Xtrain);

             % Cox regression code
             CIdata = applyCoxfit(Xtrain,Y,Xtrain,Y,censoring,censoring);

             CItemp = 0;
             for n = 1:nBoot
                 Xtrain = X(trainSets(:,n),indexModel); Xtest = X(testSets{n},indexModel); Ytrain = Y(trainSets(:,n),1); Ytest = Y(testSets{n},1);
                 [Xtrain,Xtest] = normalizeZeroOne(Xtrain,Xtest);

                 % Cox regression code
                 CIboot = applyCoxfit(Xtrain,Ytrain,Xtest,Ytest,censoring(trainSets(:,n)),censoring(testSets{n}));

                 % FOR CI (concordance index)
                 alpha = top/(1-low*(CIdata-CIboot)/(CIdata-0.5+eps));
                 if alpha > 1
                     alpha = 1;
                 elseif alpha < top
                     alpha = top;
                 end
                 if CIboot < 0.5
                     CIboot = 0.5;
                 end
                 CItemp = CItemp + (1-alpha)*CIdata + alpha*CIboot;

             end
             CItemp = CItemp/nBoot;
             metricTemp = CItemp;
             if metricTemp >= maxMetric
                 maxMetric = metricTemp;
                 index = indLeft(k);
             end
         end
         modelMat(i,j) = index;
         metricMat(i,j) = maxMetric;
         indLeft(find(indLeft==index)) = [];
     end
end

% OBTAINING MAXIMUM RESULTS FOR EVERY MODEL ORDER (maximum from all different starters)
[~,indMax] = max(metricMat);
for i = 1:maxOrder
    nameOrder = ['Order',num2str(i)];
    models.(nameOrder).Data = X(:,modelMat(indMax(i),1:i));
    models.(nameOrder).Name = cell(i,1);
    for j = 1:i
        models.(nameOrder).Name{j} = Info{modelMat(indMax(i),j)};
    end
end

end

author: mvallieres

repo: radiomics


text:
function [ c_est, c1_err, c2_err ] = triangle_random_centralize ( n, v )

%*****************************************************************************80
%
%% TRIANGLE_RANDOM_CENTRALIZE estimates the centroid of a triangle.
%
%  Discussion:
%
%    We generate N points on the surface of a triangle in 2 dimensions.
%    We seek to estimate the centroid of the triangle.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 May 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of random points to generate.
%
%    Input, real V(2,3), the vertices of the triangle.
%
%    Output, real C_EST(M,1), the estimated centroid.
%
%    Output, real C1_ERR, the norm of the difference of C_EST and the
%    area centroid.
%
%    Output, real C2_ERR, the norm of the difference of C_EST and the
%    boundary centroid.
%
  x = triangle_surface_sample ( n, v );

  c = sum ( v, 2 ) / 3.0;
  c_est = sum ( x, 2 ) / n;
  c1_err = norm ( c - c_est );

  l1 = norm ( v(1:2,2) - v(1:2,1) );
  l2 = norm ( v(1:2,3) - v(1:2,2) );
  l3 = norm ( v(1:2,1) - v(1:2,3) );
  l_sum = l1 + l2 + l3;
  l1 = 0.5 * l1 / l_sum;
  l2 = 0.5 * l2 / l_sum;
  l3 = 0.5 * l3 / l_sum;

  c2(1:2,1) = ( l3 + l1 ) * v(1:2,1) ...
            + ( l1 + l2 ) * v(1:2,2) ...
            + ( l2 + l3 ) * v(1:2,3);

  c2_err = norm ( c2 - c_est );

  return
end


text:
function value = t_triple_product_integral ( i, j, k )

%*****************************************************************************80
%
%% T_TRIPLE_PRODUCT_INTEGRAL: integral (-1<=x<=1) T(i,x)*T(j,x)*T(k,x)/sqrt(1-x^2) dx
%
%  Discussion:
%
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 April 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    John Mason, David Handscomb,
%    Chebyshev Polynomials,
%    CRC Press, 2002,
%    ISBN: 0-8493-035509,
%    LC: QA404.5.M37.
%
%  Parameters:
%
%    Input, integer I, J, K, the polynomial indices.
%    0 <= I, J.
%
%    Output, real VALUE, the integral.
%
  if ( i < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'T_TRIPLE_PRODUCT_INTEGRAL - Fatal error!\n' );
    fprintf ( 1, '  0 <= I is required.\n' );
    error ( 'T_TRIPLE_PRODUCT_INTEGRAL - Fatal error!' );
  end

  if ( j < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'T_TRIPLE_PRODUCT_INTEGRAL - Fatal error!\n' );
    fprintf ( 1, '  0 <= J is required.\n' );
    error ( 'T_TRIPLE_PRODUCT_INTEGRAL - Fatal error!' );
  end

  if ( k < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'T_TRIPLE_PRODUCT_INTEGRAL - Fatal error!\n' );
    fprintf ( 1, '  0 <= K is required.\n' );
    error ( 'T_TRIPLE_PRODUCT_INTEGRAL - Fatal error!' );
  end

  value = 0.5 * ( ...
      t_double_product_integral (       i + j,   k ) + ...
    + t_double_product_integral ( abs ( i - j ), k ) );

  return
end


text:
function qrdemo (matrix_name)
%QRDEMO a simple demo of SuiteSparseQR
% This is the MATLAB equivalent of the qrdemo.cpp program.  It reads in either
% a single matrix or a set of matrices, held as files in Matrix Market format
% in the ../Matrix directory.
%
% Example
%   qrdemo
%   qrdemo ('../Matrix/young1c.mtx')
%
% See also spqr, spqr_solve, mldivide

%   Copyright 2008, Timothy A. Davis
%   http://www.cise.ufl.edu/research/sparse

if (nargin == 0)

    % try all matrices
    qrdemo ('../Matrix/a2.mtx') ;
    qrdemo ('../Matrix/r2.mtx') ;
    qrdemo ('../Matrix/a04.mtx') ;
    qrdemo ('../Matrix/a2.mtx') ;
    qrdemo ('../Matrix/west0067.mtx') ;
    qrdemo ('../Matrix/c2.mtx') ;
    qrdemo ('../Matrix/a0.mtx') ;
    qrdemo ('../Matrix/lfat5b.mtx') ;
    qrdemo ('../Matrix/bfwa62.mtx') ;
    qrdemo ('../Matrix/LFAT5.mtx') ;
    qrdemo ('../Matrix/b1_ss.mtx') ;
    qrdemo ('../Matrix/bcspwr01.mtx') ;
    qrdemo ('../Matrix/lpi_galenet.mtx') ;
    qrdemo ('../Matrix/lpi_itest6.mtx') ;
    qrdemo ('../Matrix/ash219.mtx') ;
    qrdemo ('../Matrix/a4.mtx') ;
    qrdemo ('../Matrix/s32.mtx') ;
    qrdemo ('../Matrix/c32.mtx') ;
    qrdemo ('../Matrix/lp_share1b.mtx') ;
    qrdemo ('../Matrix/a1.mtx') ;
    qrdemo ('../Matrix/GD06_theory.mtx') ;
    qrdemo ('../Matrix/GD01_b.mtx') ;
    qrdemo ('../Matrix/Tina_AskCal_perm.mtx') ;
    qrdemo ('../Matrix/Tina_AskCal.mtx') ;
    qrdemo ('../Matrix/GD98_a.mtx') ;
    qrdemo ('../Matrix/Ragusa16.mtx') ;
    qrdemo ('../Matrix/young1c.mtx') ;
    fprintf ('All tests passed\n') ;

else

    %---------------------------------------------------------------------------
    % turn off warnings
    %---------------------------------------------------------------------------

    s = warning ('query', 'all') ;
    warning ('off', 'MATLAB:singularMatrix') ;
    warning ('off', 'MATLAB:rankDeficientMatrix') ;

    % --------------------------------------------------------------------------
    % the MATLAB equivalent of qrdemo.cpp
    % --------------------------------------------------------------------------

    fprintf ('qrdemo %s\n', matrix_name) ;
    A = mread (matrix_name) ;
    [m n] = size (A) ;
    anorm = norm (A,1) ;
    fprintf ('Matrix %6d-by-%-6d nnz: %6d ', m, n, nnz(A)) ;
    B = ones (m,1) ;
    [X, info] = spqr_solve (A,B) ;
    rnk = info.rank_A_estimate ;
    rnorm = norm (A*X-B,1) ;
    xnorm = norm (X,1) ;
    if (m <= n && anorm > 0 && xnorm > 0)
        rnorm = rnorm / (anorm * xnorm) ;
    end
    fprintf ('residual: %8.1e rank: %6d', rnorm, rnk) ;

    %---------------------------------------------------------------------------
    % compare with MATLAB
    % --------------------------------------------------------------------------

    try
        if (m ~= n || rnk == min (m,n))
            X = A\B ;
        else
            % The matrix is square and rank deficient.  Rather than getting a
            % divide-by-zero, force MATLAB to use a sparse QR by appending on a
            % blank row, in an attempt to find a consistent solution.
            X = [A ; sparse(1,n)] \ [B ; 0] ;
            X = X (1:m,:) ;
        end
        rnorm_matlab = norm (A*X-B,1) ;
        xnorm = norm (X,1) ;
        if (m <= n && anorm > 0 && xnorm > 0)
            rnorm_matlab = rnorm_matlab / (anorm * xnorm) ;
        end
    catch                                                                   %#ok
        % the above will fail for complex matrices that are rectangular and/or
        % rank deficient
        rnorm_matlab = nan ;
    end

    % find the true rank; this is costly
    rank_exact = rank (full (A)) ;

    if (rnk ~= rank_exact)
        % This can happen for some matrices; rnk is just an estimate, not
        % guaranteed to be equal to the true rank.  This does not occur
        % for any matrix in the ../Matrix test set.
        fprintf (' ?\n ... rank estimate is not exact; true rank is %d\n', ...
            rank_exact) ;
        warning ('rank estimate is incorrect') ;                            %#ok
    elseif (m <= n && rnk == m)
        % A*x=b is a full-rank square or under-determined system which spqr
        % should solve with low relative residual.
        if (rnorm > 1e-12)
            fprintf (' ?\n ... spqr solution is poor\n') ;
            error ('qrdemo failed') ;
        else
            fprintf (' OK\n') ;
        end
    elseif (m > n && rnk == n)
        % A*x=b is a full-rank least-squares problem.  Both spqr and MATLAB
        % should find a similar residual.
        err = abs (rnorm - rnorm_matlab) / max (rnorm_matlab, 1) ;
        if (err > 1e-12)
            fprintf (' ?\n ... spqr solution is poor\n') ;
            error ('qrdemo failed') ;
        else
            fprintf (' OK\n') ;
        end
    else
        % A*x=b is not full-rank; any solution will do ...
        fprintf (' %8.1e OK\n', rnorm_matlab) ;
    end

    %---------------------------------------------------------------------------
    % restore warning status
    %---------------------------------------------------------------------------

    warning (s) ;

end


text:
function [ mHC ] = CircularExtension2D( mH, numRows, numCols )

kernelRadiusV = floor(size(mH, 1) / 2); %<! Vertical Radius
kernelRadiusH = floor(size(mH, 2) / 2); %<! Horizontal Radius

mHC = mH;
mHC(numRows, numCols) = 0;
mHC = circshift(mHC, [-kernelRadiusV, -kernelRadiusH]);

end

author: RoyiAvital

repo: StackExchangeCodes


text:
function str = scn_get_datetime(varargin)
% :Usage:
% ::
%
%     str = scn_get_datetime
%
% pass in 'ymd' to get the string in yyyy_mm_dd-HH_MM format, so that
% alphanumeric order will correspond to chronological order
%
% Returns a string with the date and time
% Useful for annotating data and output

if strmatch('ymd', varargin)
    datestring=datestr(date, 'yyyy-mm-dd');

else
    datestring=date;
end

    cval = fix(clock);

    str = sprintf('%s_%02d:%02d', datestring, cval(4),  cval(5));

    str(str == ':') = '_';
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function C = cs_transpose (A)                                               %#ok
%CS_TRANSPOSE transpose a sparse matrix.
%   C = cs_transpose(A), computes C = A'
%   C = cs_transpose(A,-1) computes C=A.'
%   C = cs_transpose(A,1) computes C=A'
%
%   Example:
%       Prob = UFget ('HB/ibm32') ; A = Prob.A ;
%       C = cs_transpose (A) ;
%       C-A'
%
%   See also TRANSPOSE, CTRANSPOSE.

%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

error ('cs_transpose mexFunction not found') ;


text:
function [meantheta,R,sigma,confangle,kappa]=anglemean(theta);
% calculates the mean of angles
%
%    [meantheta,anglestrength,sigma,confangle,kappa]=anglemean(theta);
%
% anglestrength: can be thought of as the inverse variance. [varies between 0 and one]
% sigma: circular standard deviation
% confangle: a 95% confidence angle (confidence of the mean value)
% kappa: an estimate of kappa in the Von Mises distribution
%
% check: http://www.cosy.sbg.ac.at/~reini/huber_dutra_freitas_igarss_01.pdf
%
% Aslak Grinsted 2002
theta=mod(theta(:),2*pi);
n=length(theta);
S=sum(sin(theta));
C=sum(cos(theta));
meantheta=atan2(S,C);

if nargout<2
    return
end
%if ((S>C)&(C>0))
%    meantheta=atan(S/C);
%elseif (C<0)
%    meantheta=atan(S/C)+pi;
%else
%    meantheta=atan(S/C)+2*pi;
%end
%meantheta=mod(meantheta,2*pi);
Rsum=sqrt(S^2+C^2);
R=Rsum/n;

if (R<.53)
    kappa=2*R+R^3+5*R^5/6;
elseif (R<.85)
    kappa=-0.4+1.39*R+0.43/(1-R);
else
    kappa=1/(R^3-4*R^2+3*R);
end

%
%
% circular standard deviation:
sigma=sqrt(-2*log(R));

if nargout<4
    return
end

%conflim=.95;
%a=(length(theta)-Rsum)/chi2pdf()

%this is true if the
chi2=3.841; % = chi2inv(.95,1)
if ((R<.9)&(R>sqrt(chi2/(2*n))))
    confangle=acos(sqrt(2*n*(2*Rsum^2-n*chi2)/(4*n-chi2))/Rsum);
elseif (R>.9)
    confangle=acos(sqrt(n^2-(n^2-Rsum^2)*exp(chi2/n))/Rsum);
else %R is really really small ...
    confangle=pi/2;
    warning('Confidence angle not well determined.')
    %this is not good, but not important because the confidence is so low anyway...
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function cdf = levy_cdf ( x, a, b )

%*****************************************************************************80
%
%% LEVY_CDF evaluates the Levy CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 April 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the argument of the CDF.
%    Normally, A <= X.
%
%    Input, real A, B, the parameters of the PDF.
%    0 < B.
%
%    Output, real CDF, the value of the PDF.
%
  if ( b <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  LEVY_PDF - Fatal error!\n' );
    fprintf ( 1, '  Input parameter B <= 0.0\n' );
    error ( 'LEVY_PDF - Fatal error!' );
  end

  if ( x <= a )
    cdf = 0.0;
  else
    cdf = 1.0 - error_f ( sqrt ( b / ( 2.0 * ( x - a ) ) ) );
  end

  return
end


text:
% Initialize
clear all; close all; clc;
png_export = true;
pdf_export = false;
% Create the directory for storing images
[status_code,message,message_id] = mkdir('bin');

% A Dirac DCT dictionary
dict = spx.dict.simple.dirac_dct_mtx(16);
spx.graphics.display.dictionary_atoms_as_images(dict, struct('title', 'Dirac DCT Dictionary'));
if png_export
export_fig bin\dirac_dct_mtx.png -r120 -nocrop;
end
if pdf_export
export_fig bin\dirac_dct_mtx.pdf;
end

% A Gaussian dictionary
dict = spx.dict.simple.gaussian_mtx(16, 64);
spx.graphics.display.dictionary_atoms_as_images(dict, struct('title', 'Gaussian Dictionary'));
if png_export
export_fig bin\gaussian_mtx.png -r120 -nocrop;
end
if pdf_export
export_fig bin\gaussian_mtx.pdf;
end

% A Rademacher dictionary
dict = spx.dict.simple.rademacher_mtx(16, 64);
spx.graphics.display.dictionary_atoms_as_images(dict, struct('title', 'Rademacher Dictionary'));
if png_export
export_fig bin\rademacher_mtx.png -r120 -nocrop;
end
if pdf_export
export_fig bin\rademacher_mtx.pdf;
end

author: indigits

repo: sparse-plex


text:
function [approach,action] = gui_configapproach(approach,dosave)
% bring up a modal configuration dialog for the given approach
% [Approach,Action] = gui_configapproach(Approach)
%
% In:
%   Approach : an approach; struct with fields 'paradigm' and 'parameters' (and optionally 'description' and 'name')
%              or cell array {paradigm, parameter1, parameter2, ...}
%
%   DoSave: Whether to bring up a save approach gui, after clicking okay (default: false)
%
% Out:
%   Result : a (re-)configured version of the Approach, or the unmodified input Approach (though possibly reformatted) if the user pressed 'Cancel'
%   Action : either 'OK' or 'Cancel'
%
%                           Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                           2010-10-25

if ~exist('approach','var') || isempty(approach)
    approach = gui_chooseapproach();
    if isempty(approach)
        return; end
end

if ~exist('dosave','var') || isempty(dosave)
    dosave = false; end

if iscell(approach)
    approach = struct('paradigm',approach{1}, 'parameters', {approach(2:end)}); end

calibrate_func = approach.paradigm;
if ischar(calibrate_func)
    instance = eval(calibrate_func); %#ok<NASGU>
    calibrate_func = eval('@instance.calibrate');
end

% bring up the GUI
result = arg_guidialog(calibrate_func,'params',approach.parameters,'title','BCILAB: Configure approach','Invoke',false);
if ~isempty(result)
    approach.parameters = {result};
    % get rid of hidden references...
    approach = utl_prune_handles(approach);
    if dosave
        gui_saveapproach(approach); end
    action = 'OK';
else
    action = 'Cancel';
end

repo: REST


text:
clear all; close all;
addpath(genpath('./visualize'));
tri = load('tri.mat');
aa = dir('./aflw-2000_shown/2DASL_vertices/*.mat');
tri = tri.tri;

for ii = 1:length(aa)
    img = imread(strcat('./aflw-2000_shown/oriImgs/', aa(ii).name(1:end-4), '.jpg'));
    vertex = load(strcat('./aflw-2000_shown/2DASL_vertices/', aa(ii).name(1:end-4)));
    figure
    imshow(img)
    im1 = imagesc(img);
    hold on

    vertex = vertex.vertex;
    vertex(3,:,:) = vertex(3,:,:) - min(vertex(3,:,:));
    pcshow(vertex')
    view(2)
    saveas(gca, strcat('./aflw-2000_shown/results/', 'dense_align_', aa(ii).name(1:end-4), '.jpg'))
    close all
    figure
    im1 = imagesc(img);
    hold on
    render_face_mesh_xgtu(vertex, tri);
    saveas(gca, strcat('./aflw-2000_shown/results/', 'recons_align_', aa(ii).name(1:end-4), '.jpg'))
    close all
end

print('done')

author: XgTu


text:
function [freq] = ft_datatype_freq(freq, varargin)

% FT_DATATYPE_FREQ describes the FieldTrip MATLAB structure for freq data
%
% The freq data structure represents frequency or time-frequency decomposed
% channel-level data. This data structure is usually generated with the
% FT_FREQANALYSIS function.
%
% An example of a freq data structure containing the powerspectrum for 306 channels
% and 120 frequencies is
%
%       dimord: 'chan_freq'          defines how the numeric data should be interpreted
%    powspctrm: [306x120 double]     the power spectum
%        label: {306x1 cell}         the channel labels
%         freq: [1x120 double]       the frequencies expressed in Hz
%          cfg: [1x1 struct]         the configuration used by the function that generated this data structure
%
% An example of a freq data structure containing the time-frequency resolved
% spectral estimates of power (i.e. TFR) for 306 channels, 120 frequencies
% and 60 timepoints is
%
%       dimord: 'chan_freq_time'     defines how the numeric data should be interpreted
%    powspctrm: [306x120x60 double]  the power spectum
%        label: {306x1 cell}         the channel labels
%         freq: [1x120 double]       the frequencies, expressed in Hz
%         time: [1x60 double]        the time, expressed in seconds
%          cfg: [1x1 struct]         the configuration used by the function that generated this data structure
%
% Required fields:
%   - freq, dimord, label or labelcmb
%
% Optional fields:
%   - powspctrm, fouriesspctrm, csdspctrm, cohspctrm, time, grad, elec, cumsumcnt, cumtapcnt, trialinfo
%
% Deprecated fields:
%   - <none>
%
% Obsoleted fields:
%   - <none>
%
% Revision history:
%
% (2011/latest) The description of the sensors has changed, see FT_DATATYPE_SENS
% for further information.
%
% (2008) The presence of labelcmb in case of crsspctrm became optional,
% from now on the crsspctrm can also be represented as Nchan * Nchan.
%
% (2006) The fourierspctrm field was added as alternative to powspctrm and
% crsspctrm. The fields foi and toi were renamed to freq and time.
%
% (2003v2) The fields sgn and sgncmb were renamed into label and labelcmb.
%
% (2003v1) The initial version was defined.
%
% See also FT_DATATYPE, FT_DATATYPE_COMP, FT_DATATYPE_DIP, FT_DATATYPE_FREQ,
% FT_DATATYPE_MVAR, FT_DATATYPE_RAW, FT_DATATYPE_SOURCE, FT_DATATYPE_SPIKE,
% FT_DATATYPE_TIMELOCK, FT_DATATYPE_VOLUME

% Copyright (C) 2011, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% 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$

% get the optional input arguments, which should be specified as key-value pairs
version = ft_getopt(varargin, 'version', 'latest');

if strcmp(version, 'latest')
  version = '2011';
end

if isempty(freq)
  return;
end

% do some sanity checks
assert(isfield(freq, 'freq') && (isfield(freq, 'label') || isfield(freq, 'labelcmb')), 'inconsistent freq data structure, some field is missing');
if isfield(freq, 'label')
  % it could also be that it has labelcmb instead of label
  assert(length(unique(freq.label))==length(freq.label), 'channel labels must be unique');
end

% ensure consistency between the dimord string and the axes that describe the data dimensions
freq = fixdimord(freq);

if ~isrow(freq.freq)
  freq.freq = freq.freq';
end
if isfield(freq, 'label') && ~iscolumn(freq.label)
  % this is not present if the dimord is chancmb_freq or chancmb_freq_time
  freq.label = freq.label';
end
if isfield(freq, 'time') && ~isrow(freq.time)
  freq.time = freq.time';
end
if ~isfield(freq, 'label') && ~isfield(freq, 'labelcmb')
  ft_warning('data structure is incorrect since it has no channel labels');
end

switch version
  case '2011'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % ensure that the sensor structures are up to date
    if isfield(freq, 'grad')
      freq.grad = ft_datatype_sens(freq.grad);
    end
    if isfield(freq, 'elec')
      freq.elec = ft_datatype_sens(freq.elec);
    end
    if isfield(freq, 'opto')
      freq.opto = ft_datatype_sens(freq.opto);
    end

    if isfield(freq, 'foi') && ~isfield(freq, 'freq')
      % this was still the case in early 2006
      freq.freq = freq.foi;
      freq = rmfield(freq, 'foi');
    end

    if isfield(freq, 'toi') && ~isfield(freq, 'time')
      % this was still the case in early 2006
      freq.time = freq.toi;
      freq = rmfield(freq, 'toi');
    end

    if isfield(freq, 'cumtapcnt') && isvector(freq.cumtapcnt)
      % ensure that it is a column vector
      freq.cumtapcnt = freq.cumtapcnt(:);
    end

    if isfield(freq, 'cumsumcnt') && isvector(freq.cumsumcnt)
      % ensure that it is a column vector
      freq.cumsumcnt = freq.cumsumcnt(:);
    end

    % ensure that the structure has all required fields
    % note that dimord is listed as required field, but it might also be xxxdimord, or dynamically determined with GETDIMORD
    for required={'freq'}
      assert(isfield(freq, required), 'required field "%s" is missing', required{:});
    end
    % either label or labelcmb should be present
    assert(any(ismember({'label', 'labelcmb'}, fieldnames(freq))), 'required field "label" or "labelcmb" is missing');

  case '2008'
    % there are no known conversions for backward or forward compatibility support

  case '2006'
    % there are no known conversions for backward or forward compatibility support

  case '2003v2'
    % there are no known conversions for backward or forward compatibility support

  case '2003v1'
    % there are no known conversions for backward or forward compatibility support

  otherwise
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    ft_error('unsupported version "%s" for freq datatype', version);
end


text:
% Fig. 6.25   Feedback Control of Dynamic Systems, 5e
%             Franklin, Powell, Emami
%

clear all;
close all;

num=1;
den=conv([1 0],[1 2 1]);
w=logspace(-2,3,100);
[m,p]=bode(num,den,w);
subplot(2,1,1)
loglog(w,m);
axis([.01 100 .0001 100])
hold on;
ii=[11 21 47 61];
loglog(w([ii]),m([ii]),'*');
loglog(1,.5,'*')
grid;
ylabel('Magnitude');
title('Fig. 6.25 Bode Plot for G=1/(s+1)^2 (a) magnitude');
hold off;
subplot(2,1,2)
semilogx(w,p);
axis([.01 100 -270 -90])
hold on;
w180=[.01 100];
p180=[-180 -180];
semilogx(w([11 21 47 61]),p([11 21 47 61]),'*');
semilogx(1,-180,'*')
semilogx(w180,p180)
grid;
xlabel('\omega (rad/sec)');
ylabel('Phase (deg)');
title('(b) phase ');
hold off;

author: Sable

repo: mcbench-benchmarks


text:
function [vars,eval] = km_akcca(pars,data)
% KM_AKCCA performs Alternating Kernel Canonical Correlation Analysis to
% blindly identify and equalize a single-input multiple-output Wiener.
%
% Input:	- pars: structure containing parameters of the alternating KCCA
%             algorithm
%			- data: structure containing the available data x, and the
%             unavailable data (source signal s, internal signals y, and
%             linear channels B)
% Output:	- vars: estimated variables
%           - eval: structure containing the resuls of the algorithm
%             evaluation
% USAGE: [vars,eval] = km_akcca(pars,data)
%
% Author: Steven Van Vaerenbergh (steven *at* gtas.dicom.unican.es), 2011.
%
% The algorithm in this file is based on the following publication:
% S. Van Vaerenbergh, J. Via and I. Santamaria, "Blind Identification of
% SIMO Wiener Systems based on Kernel Canonical Correlation Analysis",
% accepted for publication in IEEE Transactions on Signal Processing, 2013.
%
% This file is part of the Kernel Methods Toolbox for MATLAB.
% https://github.com/steven2358/kmbox

% AKCCA CORE
[vars, eval] = AKCCA_CORE_INIT(pars,data.x); % initialize
vars.it = 0;
eval.converged = false;
while ((vars.it < pars.it_max) && ~eval.converged)
    fprintf('.');
    vars.it = vars.it+1;
    vars = AKCCA_CORE_LINID(pars,vars); % CCA 1: estimate h and obtain z
    vars = AKCCA_CORE_NLINID(pars,vars); % CCA 2: estimate alpha and obtain y
    vars = AKCCA_CORE_EQ(pars,vars); % EQUALIZATION
    eval = AKCCA_CORE_EVAL(pars,vars,data,eval); % Error calculation and check convergence
end
fprintf('\n');

function [vars,eval] = AKCCA_CORE_INIT(pars,x)

% READ PARAMETERS AND VARIABLES
p = pars.p;
m = pars.m;
N = pars.data_N;
num_it = pars.it_max;
ktype = pars.kernel.type;
% ktype = pars.kernel{1};
kpar = pars.kernel.par;
diff_nlin = ~pars.identical_nonlin;	% boolean indicating different nonlinearities

% CONSTRUCT KERNEL MATRIX DECOMPOSITIONS
K = cell(p,1);
c = cell(p,1);
if diff_nlin,
    switch pars.decomp
        case 'KPCA'
            for i=1:p,
                if isa(kpar,'function_handle')	% determine kernel parameter
                    parfun = kpar;
                    c{i} = parfun(x{i});
                else c{i} = kpar(i);
                end
                Kfull = km_kernel_center(x{i},x{i},x{i},ktype,c{i});
                Kfull = 0.5*(Kfull+Kfull');	% avoid matlab rounding errors
                [V,D] = eig(Kfull);
                [D_sort,ind] = sort(real(diag(D)),'descend');
                if m<1,	% m represents the fraction of energy to be witheld
                    energy = sum(D_sort);
                    cumul_energy = cumsum(D_sort)/energy;
                    m = sum(cumul_energy < 1-m)+1;	% new m is number of required eigenvectors
                end
                K{i} = V(:,ind(1:m));
            end
        case 'ICD'

            for i=1:p,
                if isa(kpar,'function_handle')	% determine kernel parameter
                    parfun = kpar;
                    c{i} = parfun(x{i});
                else
                    c{i} = kpar(i);
                end
                if m<1,	% m represents the precision in ICD
                    precision = m;
                    mm = N;
                else
                    precision = 1E-10;
                    mm = m;
                end
                G = km_kernel_icd(x{i},ktype,c{i},mm,precision);
                G = G-repmat(mean(G),N,1);
                K{i} = G;
            end
        otherwise
            error 'wrong decomposition method';
    end

    y_est = x;

else
    Xf = [];
    for i=1:p,
        Xf = [Xf;x{i}]; %#ok<AGROW>
    end

    if isa(kpar,'function_handle')	% determine kernel parameter
        parfun = kpar;
        c = parfun(Xf);
    else c = kpar(i);
    end

    switch pars.decomp
        case 'KPCA'
            Kfull = km_kernel_center(Xf,Xf,Xf,ktype,c);
            Kfull = 0.5*(Kfull+Kfull');
            [V,D] = eig(Kfull);
            [D_sort,ind] = sort(real(diag(D)),'descend');

            if m<1,	% m represents the fraction of energy to be witheld
                energy = sum(D_sort);
                cumul_energy = cumsum(D_sort)/energy;
                m = sum(cumul_energy < 1-m)+1;	% new m is number of required eigenvectors
            else
                m=m*p;
            end
            for i=1:p,
                stind = (i-1)*N + 1;
                ndind = stind+N - 1;
                K{i} = V(stind:ndind,ind(1:m));
            end

        case 'ICD'

            if m<1,	% m represents the precision in ICD
                precision = m;
                mm = N;
            else
                precision = 1E-10;
                mm=m*p;
            end
            G = km_kernel_icd(Xf,ktype,c,mm,precision);
            G = G-repmat(mean(G),p*N,1);

            for i=1:p,
                stind = (i-1)*N + 1;
                ndind = stind+N - 1;
                K{i} = G(stind:ndind,:);
            end
        otherwise
            error 'wrong decomposition method';
    end

    y_est = x;
end

% WRITE EVAL AND VARIABLES
eval.MSE_y = zeros(1,num_it);
eval.MSE_z = zeros(2,num_it);
eval.MSE_h = zeros(1,num_it);
eval.result = zeros(1,num_it);

vars.y_est = y_est;		% initialize without taking into account the nonlinearity
vars.K = K;		% reduced kernel matrices
vars.x = x;		% original output data
vars.c = c;     % kernel parameter

function vars = AKCCA_CORE_LINID(pars,vars)
% estimate linear filters

% COPY PARAMETERS
p = pars.p;
L = pars.model_L;
N = pars.data_N;
K = vars.K;

y_est = vars.y_est;

% PROGRAM
n = N-L+1;

% get time-embedded data matrices
Ka = cell(p,1);
for i=1:p,
	Ka{i} = zeros(n,L);
% 	for n_ind = 1:n,    % too slow, better fill differently
% 		Ka{i}(n_ind,:) = y_est{i}(n_ind+L-1:-1:n_ind)';
% 	end
    for l_ind = 1:L,
        Ka{i}(:,l_ind) = y_est{i}(L-l_ind+1:L-l_ind+n)';
    end
end

h = SIMO_id_CCA(Ka);

zh_est = cell(p);
W = cell(p);
% rescale solution to fulfill restriction sumsqr = 1
sumnrg_h = 0;
for i=1:p,
	for j=1:p,
		if(i~=j)
			sumnrg_h = sumnrg_h + sum((Ka{i}*h{j}).^2);
		end
	end
end
for i=1:p,
	h{i} = h{i}/sqrt(sumnrg_h);
end

for i=1:p,
	for j=1:p,
		if (i~=j)
            mi = size(K{i},2);
            W{i,j} = zeros(n,mi);
			zh_est{i,j} = Ka{i}*h{j};	% system output
%             for n_ind=1:n,  % too slow, better fill differently
%                 kimat = K{i}(n_ind+L-1:-1:n_ind,:);
%                 W{i,j}(n_ind,:) = h{j}'*kimat;	% auxiliary variable
%             end
            for l_ind = 1:L,    % superfast
                W{i,j} = W{i,j} + h{j}(l_ind)*K{i}(L-l_ind+1:L-l_ind+n,:);
            end
		end
	end
end

% COPY VARIABLES
vars.h = h;		% filter estimation
vars.zh_est = zh_est;	% new output estimate
vars.W = W;		% auxiliary variable (filtered kernel matrices)

function vars = AKCCA_CORE_NLINID(pars,vars)
% estimate inverse nonlinearities

% READ PARAMETERS AND VARIABLES
p = pars.p;
diff_nlin = ~pars.identical_nonlin;	% boolean indicating different nonlinearities

K = vars.K;
W = vars.W;
reg = pars.reg;

% PROGRAM: 1 iteration
za_est = cell(p);
y_est = cell(p,1);
if diff_nlin
	% standard case: no identical nonlinearity
	a = SIMO_id_CCA_dual(W,reg);

	% rescale solution to fulfill restriction sumsqr = 1
	sumnrg_a = 0;
	for i=1:p,
		for j=1:p,
			if(i~=j)
				sumnrg_a = sumnrg_a + sum((W{i,j}*a{i}).^2);
			end
		end
	end

	for i=1:p,
		a{i} = a{i}/sqrt(sumnrg_a);
		y_est{i} = K{i}*a{i};	% get y estimates
		y_est{i} = y_est{i}/sqrt(norm(y_est{i}))*sqrt(norm(vars.x{i}));	% normalize
		for j=1:p,
			if (i~=j)
			za_est{i,j} = W{i,j}*a{i};	% system output
			end
		end
	end
else
	% identical nonlinearity
	a = SIMO_id_CCA_dual_identical(W,reg);
	for i=1:p,
		y_est{i} = K{i}*a;	% get y estimates
		for j=1:p,
			if (i~=j)
			za_est{i,j} = W{i,j}*a;	% system output
			end
		end
	end
end

% WRITE VARIABLES
vars.alpha = a;
vars.y_est = y_est;
vars.za_est = za_est;

function vars = AKCCA_CORE_EQ(pars,vars)
% equalize linear channels (zero-forcing)

% COPY PARAMETERS
p = pars.p;
L = pars.model_L;
N = pars.data_N;
k = pars.zf_k;

y_est = vars.y_est;
h = vars.h;

% PROGRAM
% use k observations to generate filter matrix and estimate equalizers
Hc = cell(p,1); Hr = cell(p,1);
TH = zeros(p*k,k+L-1);
for i=1:p,
	Hc{i} = [h{i}(1);zeros(k-1,1)];
	Hr{i} = [h{i}' zeros(1,k-1)];
	TH(i:p:end,:) = toeplitz(Hc{i},Hr{i});
end

% zero-forced equalizing
ZF = pinv(TH);

% S_est = zeros(N-k+1,k);
% Yi = zeros(p*k,1);
% for i=1:N-k+1,  % too slow, better fill differently
% 	for j=1:p,
% 		yi = y_est{j}(i+k-1:-1:i);
% 		Yi(j:p:end) = yi;
% 	end
% 	s_esti = ZF*Yi;
% 	S_est(i,:) = s_esti(1:k).';
% end
YY = zeros(N-k+1,p*k);
for j = 1:p,
    for k_ind=1:k,
        col = (k_ind-1)*p + j;
        YY(:,col) = y_est{j}(k-k_ind+1:N-k_ind+1);
    end
end
S_est = YY*ZF(1:k,:)';

cnorms = zeros(k,1);
for i=1:k,
	cnorms(i) = norm(ZF(:,i));
end
[mm,zfi] = min(cnorms); %#ok<ASGLU>
s_est = S_est(:,zfi);

inds = k-zfi+1:N-zfi+1;

% COPY VARIABLES
vars.s_est = s_est;
vars.ind_s = inds;
vars.ZF = ZF;

function eval = AKCCA_CORE_EVAL(pars,vars,data,eval)
% evaluate solutions: calculate errors

% COPY PARAMETERS
p = pars.p;
vb = pars.verbose;

y_est = vars.y_est;
s_est = vars.s_est;
h = vars.h;
it = vars.it;
inds = vars.ind_s;

s = data.s; % true source signal
B = data.B; % true channels
y = data.y; % true internal channels

% PROGRAM
MSE_y = 0; MSE_h = 0; MSE_za = 0; MSE_zh = 0;
for i=1:p,
	% check if system internals are known
	if (iscell(y))
		% MSE Y
		MSE_y = MSE_y + norm_compare(y{i},y_est{i})/p;
		% MSE H
		numz = length(h{i})-length(B{i});
		if (numz>0)
			% channel overestimation case
			B{i} = [B{i};zeros(numz,1)];
		end
		MSE_h = MSE_h + norm_compare(B{i},h{i})/p;
	end

	% MSE Z
	if isfield(vars,'za_est')
		za_est = vars.za_est;
	end
	zh_est = vars.zh_est;

	for j=i+1:p,
		if (i~=j)
			MSE_zh = MSE_zh + sum((zh_est{i,j}-zh_est{j,i}).^2);
			if isfield(vars,'za_est')
				MSE_za = MSE_za + sum((za_est{i,j}-za_est{j,i}).^2);
			end
		end
	end
end

signal_test = s(inds);
signal_est = s_est;

switch lower(pars.data_type)
	case {'gaussian'}
		% result = MSE between s and s_est
		result = norm_compare(signal_test,signal_est);
	case 'bits'
		% result = BER
		[MSEsb,sci] = scale_compare(signal_test,signal_est); %#ok<ASGLU>
		result = sum(sign(signal_test)~=sign(signal_est/sci))/length(signal_test);
        s_est = sign(signal_est/sci);
	otherwise
		disp('Unknown signal type.')
end

% OUTPUT
if vb, fprintf(1,'MSE z_a: %f\n',MSE_za); end;
if vb, fprintf(1,'MSE z_h: %f\n',MSE_zh); end;
if vb, fprintf(1,'MSE y: %f\n',MSE_y); end;
if vb, fprintf(1,'Result: %f\n',result); end;

% COPY VARIABLES
eval.MSE_h(it) = MSE_h;
eval.MSE_y(it) = MSE_y;
eval.MSE_z(1,it) = MSE_zh;
eval.MSE_z(2,it) = MSE_za;
eval.result(it) = result;
eval.s_est = s_est;

% CHECK STOP
if (it>2)
	% stop when both sub-iterations obtain consecutive equal costs
	change1 = abs(eval.MSE_z(2*it)-eval.MSE_z(2*it-2));
	change2 = abs(eval.MSE_z(2*it-1)-eval.MSE_z(2*it-3));
	if ((change1 < pars.it_stop) && (change2 < pars.it_stop))
		eval.converged = true;
		num = pars.it_max - it;
		eval.MSE_h(it+1:pars.it_max) = repmat(MSE_h,num,1);
		eval.MSE_y(it+1:pars.it_max) = repmat(MSE_y,num,1);
		eval.MSE_z(2*it+1:2*pars.it_max) = repmat(MSE_zh,2*num,1);
		eval.result(it+1:pars.it_max) = repmat(result,num,1);
	end
end

eval.finalresult = result;

function [h,alphas,betas] = SIMO_id_CCA(X)
% SIMO_ID_CCA estimates various channels of a SIMO system using CCA.
% Finds the optimal h_i such that X_i h_j = X_j H_i s.t. sum (i = j)
% ||X_i h_j||^2 = 1.
% Input:	- X: cell containing p sets of length N and dimension L (N x L)
% Output:	- h: identified channels (cell)
%			- alphas: sorted eigenvectors
%			- betas: sorted eigenvalues
%
% Steven Van Vaerenbergh 2008

p = length(X);
L = size(X{1},2);

Rf = zeros(p*L); Df = zeros(p*L);
for i=1:p,
	D = zeros(L);
	for j=1:p,
		if i~=j
			sti = (i-1)*L+1;
			ndi = i*L;
			stj = (j-1)*L+1;
			ndj = j*L;
			Rf(sti:ndi,stj:ndj) = X{j}'*X{i};
			D = D + X{j}'*X{j};
		end
	end
	Df(sti:ndi,sti:ndi) = D;
end

[alphas,betas] = eig(Rf,Df);
[betas,ind] = sort(real(diag(betas)),'ascend');
alpha = alphas(:,ind(end));
% alpha = alpha/norm(alpha);

h = cell(p,1);
for i=1:p,
	sti = (i-1)*L+1;
	ndi = i*L;
	h{i} = alpha(sti:ndi);
end

function [h,alphas,betas] = SIMO_id_CCA_dual_identical(X,reg)
% SIMO_ID_CCA estimates various channels of a SIMO system using CCA.
% Finds the optimal h_i such that X_ij h = X_ji h
% Input:
%	X: cell containing pxp data sets of length N and dimension m (N by m)
%   reg: regularization
%
% outputs:
%	h: identified channels
%	alphas: sorted eigenvectors
%	betas: sorted eigenvalues
%
% Reference: Via, Santamaria [...]
%
% Steven Van Vaerenbergh 2008

if nargin<2
    reg = 1E-8;
end

p = size(X,1);
m = size(X{1,2},2);

Rf = zeros(m); Df = zeros(m);
for i=1:p,
	for j=1:p,
		if i~=j
			Rf = Rf + X{i,j}'*X{j,i};
			Df = Df + X{i,j}'*X{i,j};
		end
	end
end
Df = Df + reg*eye(size(Df));

% solve eigenvalue problem
[alphas,betas] = eig(Rf,Df);
[betas_a,ind] = sort(real(diag(betas))); %#ok<ASGLU>
h = alphas(:,ind(end));

function [h,alphas,betas] = SIMO_id_CCA_dual(X,reg)
% SIMO_ID_CCA estimates various channels of a SIMO system using CCA.
% Finds the optimal h_i such that X_ij h_i = X_ji h_j s.t. sum (i~=j)
% ||X_ij h_i||^2 = 1.
% Input:	- X: cell containing pxp sets of length N and dimension mi (N x mi).
%           - reg: regularization
% Outputs:	- h: identified channels (cell)
%			- alphas: sorted eigenvectors
%			- betas: sorted eigenvalues
%
% Steven Van Vaerenbergh 2008

if nargin<2
    reg = 1E-8;
end

p = size(X,1);
m = zeros(p,1);
for i=1:p,
	j = mod(i,p)+1;
	m(i) = size(X{i,j},2);
end
mcum = [0;cumsum(m)];

Rf = zeros(mcum(end)); Df = zeros(mcum(end));
for i=1:p,
	D = zeros(m(i));
	for j=1:p,
		if i~=j
			sti = mcum(i)+1;
			ndi = sti+m(i)-1;
			stj = mcum(j)+1;
			ndj = stj+m(j)-1;
			Rf(sti:ndi,stj:ndj) = X{i,j}'*X{j,i};
			D = D + X{i,j}'*X{i,j};
		end
	end
	Df(sti:ndi,sti:ndi) = D + reg*eye(size(D));
end

[alphas,betas] = eig(Rf,Df);
[betas,ind] = sort(real(diag(betas)),'ascend');
alpha = alphas(:,ind(end));
% alpha = alpha/norm(alpha);

h = cell(p,1);
for i=1:p
	sti = mcum(i)+1;
	ndi = mcum(i+1);
	h{i} = alpha(sti:ndi);
end

function MSE = norm_compare(data_1,data_2)
% NORM_COMPARE calculates the squared error between 2 normalized signals

data_1_norm = data_1/norm(data_1);
data_2_norm = data_2/norm(data_2);

% take into account sign ambiguity
err_vecA = data_1_norm - data_2_norm;
err_vecB = data_1_norm + data_2_norm;

MSEA = err_vecA'*err_vecA;
MSEB = err_vecB'*err_vecB;

MSE = min(MSEA,MSEB);

function [MSE,sc] = scale_compare(data_1,data_2)

data_1_norm = data_1/norm(data_1);
sc = data_1_norm'*data_2/(data_2'*data_2);

err_vec = data_1_norm - sc*data_2;

MSE = err_vec'*err_vec;
sc = sc*norm(data_1);

author: epfl-lasa


text:
%% RATE OF CONVERGENCE OF MIXED FINITE ELEMENT METHOD (RT0-P0) FOR POISSON EQUATION
%
% This example is to show the rate of convergence of RT0-P0 mixed finite
% element approximation of the Poisson equation on the unit square with the
% following boundary conditions:
%
% - pure Dirichlet boundary condition.
% - Pure Neumann boundary condition.
% - mxied boundary condition.
%
% The basis, data structure and numerical test is summarized in <a
% href="matlab:ifem PoissonRT0femrate">PoissonRT0femrate</a>.
%
% The optimal rates of convergence for u and sigma are observed, namely,
% 1st order for L2 norm of u, L2 norm of sigma and H(div) norm of sigma.
% The 2nd order convergent rates between two discrete functions ||uI-uh||
% and ||sigmaI-sigmah|| are known as superconvergence.
%
% Triangular preconditioned GMRES (default solver) and Uzawa preconditioned
% CG converges uniformly in all cases. Traingular preconditioner is two
% times faster than PCG although GMRES is used.
% See also PoissonBDM1mfemrate, Poissonfemrate, Poissonafemrate
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

%% Setting
[node,elem] = squaremesh([0,1,0,1],0.25);
mesh = struct('node',node,'elem',elem);
option.L0 = 1;
option.maxIt = 4;
option.printlevel = 1;
option.elemType = 'RT0';
option.solver = 'tri';
pde = sincosNeumanndata;

%% Pure Neumann boundary condition.
mesh.bdFlag = setboundary(node,elem,'Neumann');
mfemPoisson(mesh,pde,option);

%% Pure Dirichlet boundary condition.
mesh.bdFlag = setboundary(node,elem,'Dirichlet');
mfemPoisson(mesh,pde,option);

%% Mix Dirichlet and Neumann boundary condition.
option.solver = 'uzawapcg';
mesh.bdFlag = setboundary(node,elem,'Dirichlet','~(x==0)','Neumann','x==0');
mfemPoisson(mesh,pde,option);

repo: ifem


text:
function clearDirectory(rootdir)

if exist(fullfile(rootdir,'movie.avi'),'file')

  rmdir(rootdir,'s');

else

  dd = dir(rootdir);
  for ndx =1:numel(dd)
    if strcmp(dd(ndx).name(1),'.'), continue; end
    curd = fullfile(rootdir,dd(ndx).name);
    if isdir(curd)
       clearDirectory(curd);
    end

  end

end

author: kristinbranson

repo: JAABA


text:
%--- help for generic/is_stable_system ---
%
%  Checks the stability of a linear markov switching system.
%
%  ::
%
%     flag = is_stable_system(obj)
%     flag = is_stable_system(obj,varargin)
%
%  Args:
%
%     obj (dsge | rise | svar | rfvar): model object
%
%     varargin (name,value): pairwise valid options for RISE. The most
%       relevant in this case are
%
%        - **stability_criterion** [numeric\|{1.000001}]: stability criterion.
%          All eigenvalues must be smaller than this criterion for the system to
%          be MSS
%        - **stability_algorithm** ['cfm'\|{'hmg'}]: CFM stands for
%          Costa-Fragoso-Marques while HMG stands for Hassibi-Murray-Gupta.
%
%  Returns:
%     :
%
%     - **flag** [false\|true]: result of the investigation on whether the
%       system is stable or not.
%
%  Note:
%
%     RISE implements two algorithms from the engineering literature to check
%     for the stability. They are
%
%        - Costa-Fragoso-Marques :cite:`costa2006discrete`
%        - Hassibi-Murray-Gupta :cite:`gupta2003control`
%
%     Refer to the references to see the specific algorithms. However, for most
%     applications, one can just use the default options.
%
%  References:
%
%     - :cite:`costa2006discrete`
%     - :cite:`gupta2003control`
%
%

author: jmaih


text:
function out = spm_run_realignunwarp(job)
% SPM job execution function
% takes a harvested job data structure and call SPM functions to perform
% computations on the data.
% Input:
% job    - harvested job data structure (see matlabbatch help)
% Output:
% out    - computation results, usually a struct variable.
%__________________________________________________________________________
% Copyright (C) 2005-2013 Wellcome Trust Centre for Neuroimaging

% Darren R. Gitelman
% $Id: spm_run_realignunwarp.m 6554 2015-09-11 17:21:23Z guillaume $

%-Assemble flags
%--------------------------------------------------------------------------
flags.quality = job.eoptions.quality;
flags.fwhm    = job.eoptions.fwhm;
flags.sep     = job.eoptions.sep;
flags.rtm     = job.eoptions.rtm;
flags.PW      = char(job.eoptions.weight);
flags.interp  = job.eoptions.einterp;
flags.wrap    = job.eoptions.ewrap;

uweflags.order     = job.uweoptions.basfcn;
uweflags.regorder  = job.uweoptions.regorder;
uweflags.lambda    = job.uweoptions.lambda;
uweflags.jm        = job.uweoptions.jm;
uweflags.fot       = job.uweoptions.fot;

if ~isempty(job.uweoptions.sot)
    cnt = 1;
    for i=1:size(job.uweoptions.sot,2)
        for j=i:size(job.uweoptions.sot,2)
            sotmat(cnt,1) = job.uweoptions.sot(i);
            sotmat(cnt,2) = job.uweoptions.sot(j);
            cnt = cnt+1;
        end
    end
else
    sotmat = [];
end
uweflags.sot       = sotmat;
uweflags.fwhm      = job.uweoptions.uwfwhm;
uweflags.rem       = job.uweoptions.rem;
uweflags.noi       = job.uweoptions.noi;
uweflags.exp_round = job.uweoptions.expround;

uwrflags.interp    = job.uwroptions.rinterp;
uwrflags.wrap      = job.uwroptions.wrap;
uwrflags.mask      = job.uwroptions.mask;
uwrflags.which     = job.uwroptions.uwwhich(1);
uwrflags.mean      = job.uwroptions.uwwhich(2);
uwrflags.prefix    = job.uwroptions.prefix;

if uweflags.jm == 1
    uwrflags.udc = 2;
else
    uwrflags.udc = 1;
end

%-Assemble files
%--------------------------------------------------------------------------
P   = cell(size(job.data));
sfP = cell(size(job.data));
for i = 1:numel(job.data)
    P{i} = char(job.data(i).scans{:});
    if ~isempty(job.data(i).pmscan)
        sfP{i} = job.data(i).pmscan{1};
    else
        sfP{i} = [];
    end
end

%-Realign
%--------------------------------------------------------------------------
spm_realign(P,flags);

%-Unwarp Estimate
%--------------------------------------------------------------------------
for i = 1:numel(P)
    uweflags.sfP = sfP{i};
    P1 = deblank(P{i}(1,:));
    if isempty(spm_file(P1,'number'))
        P1 = spm_file(P1,'number',1);
    end
    VP1 = spm_vol(P1);
    uweflags.M = VP1.mat;
    ds = spm_uw_estimate(P{i},uweflags);
    sess(i).ds = ds;
    dsfile = spm_file(P{i}(1,:), 'suffix','_uw', 'ext','.mat');
    save(dsfile,'ds', spm_get_defaults('mat.format'));
end

%-Unwarp Write - Sessions should be within subjects
%--------------------------------------------------------------------------
spm_uw_apply(cat(2,sess.ds),uwrflags);

%-Dependencies
%--------------------------------------------------------------------------
for i=1:numel(P)
    out.sess(i).rpfile{1} = spm_file(P{i}(1,:), 'prefix','rp_', 'ext','.txt');
    out.sess(i).dsfile{1} = spm_file(P{i}(1,:), 'suffix','_uw', 'ext','.mat');
end

switch job.uwroptions.uwwhich(1)
    case 0
        out.sess.uwrfiles  = {};
    case 2
        for i=1:numel(P)
            out.sess(i).uwrfiles = spm_file(cellstr(P{i}), ...
                'prefix',job.uwroptions.prefix);
        end
end
if job.uwroptions.uwwhich(2)
    out.meanuwr{1} = spm_file(P{1}(1,:), ...
        'prefix',['mean' job.uwroptions.prefix], ...
        'number','');
end

author: spm


text:
function lm = energy_awgn_normapx(en0, epsil)
% Compute normapx on log m^*(E/N_0, \epsilon)

loge = log2(exp(1));
lm = en0 * loge + loge * sqrt(2*en0) * norminv(epsil) + log2(en0)/2;
%lm = en0 * loge + loge * sqrt(2*en0) * norminv(epsil) + log2(2*en0)/2;

author: yp-mit

repo: spectre

save_path: github-repos/MATLAB/yp-mit-spectre


text:
function value = r8_gamic ( a, x )

%*****************************************************************************80
%
%% R8_GAMIC evaluates the complementary incomplete gamma function.
%
%  Discussion:
%
%    GAMIC = integral ( x <= t < oo ) exp(-t) * t^(a-1) dt
%
%    GAMIC is evaluated for arbitrary real values of A and non-negative
%    values X (even though GAMIC is defined for X < 0.0), except that
%    for X = 0 and A <= 0.0, GAMIC is undefined.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 October 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%    Walter Gautschi,
%    A Computational Procedure for Incomplete Gamma Functions,
%    ACM Transactions on Mathematical Software,
%    Volume 5, Number 4, December 1979, pages 466-481.
%
%  Parameters:
%
%    Input, real A, the parameter.
%
%    Input, real X, the evaluation point.
%
%    Output, real VALUE, the value of the incomplete gamma function.
%
  persistent alneps
  persistent bot
  persistent eps
  persistent sqeps

  if ( isempty ( eps ) )
    eps = 0.5 * r8_mach ( 3 );
    sqeps = sqrt ( r8_mach ( 4 ) );
    alneps = - log ( r8_mach ( 3 ) );
    bot = log ( r8_mach ( 1 ) );
  end

  if ( x < 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8_GAMIC - Fatal error!\n' );
    fprintf ( 1, '  X < 0.\n' );
    error ( 'R8_GAMIC - Fatal error!' )
  end

  if ( x == 0.0 )

    if ( a <= 0.0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'R8_GAMIC - Fatal error!\n' );
      fprintf ( 1, '  X = 0 and A <= 0.\n' );
      error ( 'R8_GAMIC - Fatal error!' )
    end

    value = exp ( r8_lngam ( a + 1.0 ) - log ( a ) );

    return
  end

  alx = log ( x );
  if ( a < 0.0 )
    sga = - 1.0;
  else
    sga = + 1.0;
  end

  ainta = r8_aint ( a + 0.5 * sga );
  aeps = a - ainta;

  izero = 0;

  if ( x < 1.0 )

    if ( a <= 0.5 && abs ( aeps ) <= 0.001 )

      if ( - ainta <= 1.0 )
        e = 2.0;
      else
        e = 2.0 * ( - ainta + 2.0 ) / ( ainta * ainta - 1.0 );
      end

      e = e - alx * x^( - 0.001 );

      if ( e * abs ( aeps ) <= eps )
        value = r8_gmic ( a, x, alx );
        return
      end

    end

    [ algap1, sgngam ] = r8_lgams ( a + 1.0 );
    gstar = r8_gmit ( a, x, algap1, sgngam, alx );

    if ( gstar == 0.0 )
      izero = 1;
    else
      alngs = log ( abs ( gstar ) );
      sgngs = r8_sign ( gstar );
    end

  else

    if ( a < x )
      value = exp ( r8_lgic ( a, x, alx ) );
      return
    end

    sgngam = 1.0;
    algap1 = r8_lngam ( a + 1.0 );
    sgngs = 1.0;
    alngs = r8_lgit ( a, x, algap1 );

  end

  h = 1.0;

  if ( izero ~= 1 )

    t = a * alx + alngs;

    if ( alneps < t )
      sgng = - sgngs * sga * sgngam;
      t = t + algap1 - log ( abs ( a ) );
      value = sgng * exp ( t );
      return
    end

    if ( - alneps < t )
      h = 1.0 - sgngs * exp ( t );
    end

  end

  sgng = r8_sign ( h ) * sga * sgngam;
  t = log ( abs ( h ) ) + algap1 - log ( abs ( a ) );
  value = sgng * exp ( t );

  return
end

author: johannesgerer


text:
function [ CatchAline ] = RunElastic2D_SteelBlock(  )
% SoundSim - 2D Elastic Simulation - For Academic Use Only
% -------------------------------------------------------------------------------------
% Limited Feature Beta Release of a SoundSim (www.SoundSim.com)
% Please Send all Questions, Suggestions, and Comments to kevinrudd@SoundSim.com
% No Error Checking - May become unstable if incorrect parameters are used
% -------------------------------------------------------------------------------------
% Steel Block Example -
% A 0.8 MHz Impulse on a top of a 7cm by 7cm steel block.
% The longitudinal, sheer, and Raleigh waves are visible
% The mode conversions are also visible as these waves
% reflect from the edges of the block

%                                                                                Aprox. Material Values From
%                                                                               Stress Waves in Solids (H. Kolsky) Dover 1963
% Material Properties                                                                   Aluminum  Steel  Copper  Glass
p.cl = 5940;                         % Longitudinal Wave Speed (m/s)                     6320      5940   4560    5800
p.cs = 3220;                         % Sheer Wave Speed (m/s)                            3100      3220   2250    3350
p.den = 7800;                        % Material density (kg/m^3)                         2700      7800   8900    2500
p.platethickness = 0.07;             % block thickness (meters)  (in x2-direction)
p.platelength    = 0.07;             % block length (meters)     (in x1-direction)

% Driving Transducer Properties
% -- These parameters adust the transducer position, size, frequency, and drive function length
% -- of the driving transducer.  The transducer is on the top of the plate.
p.tpos = 0.025;                      % transducer position  (meters) - from left side of the block
p.tthickness = 0.001;                % transducer diameter  (meters) - this is 2D
p.tfreq = 800000;                    % frequency of transducer (Hertz)
p.tpulselength = (.5)*(1/(p.tfreq)); % transducer pulse length (seconds)

% Catch Transducer Properties
% -- These parameters adjust the catch transducer location and size.  Again, this transducer
% -- can only be placed on the top of the block
p.catchtpos = 0.05;                  % catch transducer position  (meters) - from left side of the plate
p.catchtthickness = 0.001;           % catch transducer diameter  (meters) - this is 2D

% Other Params
p.abc = 0;                           % thickness of absorbing boundary layer (in simulation units ds) Only on left and right
p.plotevery = 3;                     % update the plot every <plotevery> timesteps
p.SimulationTime =.000013;           % siulation runtime (seconds)

p.plotmode = 'abs_velocity';         % plot mode ('abs_velocity' = absolute velocity);
                                     %           ('x1_velocity'  = velocity in the x1 position (left and right))
                                     %           ('x2_velocity'  = velocity in the x2 position (up and down))

CatchAline=SoundSim_ElasticEngine2D( p );  % Run the Simulation


text:
function SHOWELEMENTS(ShowElements,ShowNodes,coordinates,X,Y,nel,nodes,nnode)
%--------------------------------------------------------------------------
%   Purpose:
%           To display and undisplay the Element numbers
%--------------------------------------------------------------------------
ShowElements = get(ShowElements,'Value') ;
ShowNodes = get(ShowNodes,'Value') ;
% Diaply only Element Numbers
if ShowElements==1 && ShowNodes==0
    cla(gcf) ;
    patch(X,Y,'w') ;
    for i = 1:nel
        EX = coordinates(nodes(i,:),1) ;EY = coordinates(nodes(i,:),2) ;
        pos = [sum(EX)/4,sum(EY)/4] ;
        text(pos(1),pos(2),int2str(i),'fontsize',8, 'fontweight','bold','color','b');
    end
end
% Display only Node Numbers
if ShowElements==0 && ShowNodes==1
    cla(gcf) ;
    patch(X,Y,'w')
    for i = 1:nnode
        text(coordinates(i,1),coordinates(i,2),int2str(i),'fontsize',8,....
            'fontweight','bold','Color','r');
    end
end
% Display both Element Numbers and Node Numbers
if ShowElements==1 && ShowNodes==1
    cla(gcf) ;
    patch(X,Y,'w') ;
    for i = 1:nel
        EX = coordinates(nodes(i,:),1) ;EY = coordinates(nodes(i,:),2) ;
        pos = [sum(EX)/4,sum(EY)/4] ;
        text(pos(1),pos(2),int2str(i),'fontsize',8, 'fontweight','bold','color','b');
    end
    for i = 1:nnode
        text(coordinates(i,1),coordinates(i,2),int2str(i),'fontsize',8,....
            'fontweight','bold','Color','r');
    end
end
% Display None
if ShowElements==0 && ShowNodes==0
    cla(gcf) ;
    patch(X,Y,'w')
end
return


text:
%
%M-file for compiling and linking the FGI source files
%
%Date: 01-06-2000
%Copyright (c) 2000 R.Heil, C. L. Wauters
%SYNTAX buildall;
%
%
%DESCRIPTION File for easily compiling and linking the FGI source files.
%
%INPUT The user is asked to input the source directory and output directory
%
%EXAMPLE:
%buildall;
%

INC = '-Iy:\software\dt\fgsdk32\include';
INC52='-Iy:\software\dt\fgsdk32\include\dt3152';
INC57='-Iy:\software\dt\fgsdk32\include\dt3157';
INCCO='-Iy:\software\dt\fgsdk32\include\colorsdk';
LIB = 'y:\software\dt\fgsdk32\lib\olimg32.lib y:\software\dt\fgsdk32\lib\olfg32.lib';

disp(' ');
%SOURCE = input('Enter the source directory (example: g:\\sourcedir ) -> ','s');
SOURCE = '..\Src';
%OUTDIR = input('Enter the output directory (example: g:\\matlab ) -> ','s');
OUTDIR = '..\Bin';
SOURCE = strcat(SOURCE,'\\');

disp(' ');
disp('Building FGI source files');
disp(' ');

disp('Compiling setfg');
com = strcat(SOURCE,'setfg.c');
mex(com,'-O',INC,INC52,INC57,INCCO,LIB,'-outdir',OUTDIR)

disp('Compiling getfg');
com = strcat(SOURCE,'getfg.c');
mex(com,'-O',INC,INC52,INC57,INCCO,LIB,'-outdir',OUTDIR)

disp('Compiling loadfg');
com = strcat(SOURCE,'loadfg.c');
mex(com,'-O',INC,INC52,INC57,INCCO,LIB,'-outdir',OUTDIR)

disp('Compiling resetfg');
com = strcat(SOURCE,'resetfg.c');
mex(com,'-O',INC,INC52,INC57,INCCO,LIB,'-outdir',OUTDIR)

disp('Compiling grabfg');
com = strcat(SOURCE,'grabfg.c');
mex(com,'-O',INC,INC52,INC57,INCCO,LIB,'-outdir',OUTDIR)

disp('Compiling isopenfg');
com = strcat(SOURCE,'isopenfg.c');
mex(com,'-O',INC,INC52,INC57,INCCO,'-outdir',OUTDIR)

disp('Compiling gethandlefg');
com = strcat(SOURCE,'gethandlefg.c');
mex(com,'-O',INC,INC52,INC57,INCCO,'-outdir',OUTDIR)

disp('Compiling openfg');
com = strcat(SOURCE,'openfg.c');
mex(com,'-O',INC,INC52,INC57,INCCO,LIB,'-outdir',OUTDIR)

disp('Compiling closefg');
com = strcat(SOURCE,'closefg.c');
mex(com,'-O',INC,INC52,INC57,INCCO,LIB,'-outdir',OUTDIR)

disp('Compiling mainfg');
com = strcat(SOURCE,'mainfg.c');
mex(com,'-O',INC,INC52,INC57,INCCO,'-outdir',OUTDIR)


text:
function test_suite=test_align
% tests for cosmo_align
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function test_align_basics
    x=cosmo_synthetic_dataset('type','meeg','size','big');
    y=cosmo_slice(x,x.fa.chan<=6,2);
    orig_ds=cosmo_dim_transpose(y,'time');

    [nsamples,nfeatures]=size(orig_ds.samples);

    ds1=cosmo_slice(orig_ds,randperm(nsamples));
    ds2=cosmo_slice(orig_ds,randperm(nsamples));

    [mp,pm]=cosmo_align({ds1.sa.targets,ds1.sa.chunks,ds1.sa.time},...
                            {ds2.sa.targets,ds2.sa.chunks,ds2.sa.time});
    assertEqual(cosmo_slice(ds1,mp),ds2);
    assertEqual(cosmo_slice(ds2,pm),ds1);

    [mp,pm]=cosmo_align([2 3 4],[4 2 3]);
    assertEqual(mp,[3 1 2]);
    assertEqual(pm,[2 3 1]);

    [mp,pm]=cosmo_align({{'b','c','c'},[3 2 3]},...
                            {{'c','b','c'},[3 3 2]});
    assertEqual(mp,[3 1 2]);
    assertEqual(pm,[2 3 1]);

    % test structs
    p=struct();
    p.x={'b','c','c'};
    p.y=[3 2 3];
    q=struct();
    q.y=[3 3 2];
    q.x={'c','b','c'};
    [mp,pm]=cosmo_align(p,q);
    assertEqual(mp,[3 1 2]);
    assertEqual(pm,[2 3 1]);

    % test NaN
    [mp,pm]=cosmo_align([2 NaN 4],[4 2 NaN]);
    assertEqual(mp,[3 1 2]);
    assertEqual(pm,[2 3 1]);

    % test exceptions
    ds_small=cosmo_slice(ds1,1:nsamples-1);
    aet=@(varargin)assertExceptionThrown(@()...
                    cosmo_align(varargin{:}),'');
    aet(ds1,ds2);
    aet([2 3 4 3],[4 2 3 4]);
    aet([2 3 4 3],[4 2 3 4]);
    aet([2 3 4 3],[4 2 3]);
    aet([2 3 4],{[2,3,4],[2,3,4]});

    aet([2 3 NaN NaN],[2 3 NaN Inf]);
    aet([2 NaN 3],[4 2 NaN]);

    aet(struct,struct('a',1));
    aet(struct('a',1),1);

    aet({'b','c','d','d'},{'d','b','c','d'});

function test_align_multiple_nans
    n_rows_half=20+ceil(rand()*20);
    x=ceil(rand(n_rows_half*2,2)*sqrt(n_rows_half));
    rp=randperm(2*n_rows_half);

    rp1=rp(1:n_rows_half);
    rp2=rp(n_rows_half+(1:n_rows_half));

    % half of the rows become NaN
    x(rp1,1)=NaN;
    x(rp1,2)=1:n_rows_half;
    x(rp2,1)=1:n_rows_half;
    x(rp2,2)=1:n_rows_half;

    rp2=randperm(n_rows_half*2);
    y=x(rp2,:);

    assertExceptionThrown(@()cosmo_align(...
                            {x(:,1),x(:,2)},{y(:,1),y(:,2)}),'');

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function varargout=notBoxPlot(y,x,jitter,style)
% notBoxPlot - Doesn't plot box plots!
%
% function notBoxPlot(y,x,jitter,style)
%
%
% Purpose
% An alternative to a box plot, where the focus is on showing raw
% data. Plots columns of y as different groups located at points
% along the x axis defined by the optional vector x. Points are
% layed over a 1.96 SEM (95% confidence interval) in red and a 1 SD
% in blue. The user has the option of plotting the SEM and SD as a
% line rather than area. Raw data are jittered along x for clarity. This
% function is suited to displaying data which are normally distributed.
% Since, for instance, the SEM is meaningless if the data are bimodally
% distributed.
%
%
% Inputs
% y - each column of y is one variable/group. If x is missing or empty
%     then each column is plotted in a different x position.
%
% x - optional, x axis points at which y columns should be
%     plotted. This allows more than one set of y values to appear
%     at one x location. Such instances are coloured differently.
% Note that if x and y are both vectors of the same length this function
% behaves like boxplot (see Example 5).
%
% jitter - how much to jitter the data for visualization
%          (optional). The width of the boxes are automatically
%          scaled to the jitter magnitude.
%
% style - a string defining plot style of the data.
%        'patch' [default] - plots SEM and SD as a box using patch
%                objects.
%        'line' - create a plot where the SD and SEM are
%                constructed from lines.
%        'sdline' - a hybrid of the above, in which only the SD is
%                replaced with a line.
%
%
% Outputs
% H - structure of handles for plot objects.
%
%
% Example 1 - simple example
% clf
% subplot(2,1,1)
% notBoxPlot(randn(20,5));
% subplot(2,1,2)
% h=notBoxPlot(randn(10,40));
% d=[h.data];
% set(d(1:4:end),'markerfacecolor',[0.4,1,0.4],'color',[0,0.4,0])
%
% Example 2 - overlaying with areas
% clf
% x=[1,2,3,4,5,5];
% y=randn(20,length(x));
% y(:,end)=y(:,end)+3;
% y(:,end-1)=y(:,end-1)-1;
% notBoxPlot(y,x);
%
% Example 3 - lines
% clf
% H=notBoxPlot(randn(20,5),[],[],'line');
% set([H.data],'markersize',10)
%
% Example 4 - mix lines and areas [note that the way this function
% sets the x axis limits can cause problems when combining plots
% this way]
%
% clf
% h=notBoxPlot(randn(10,1)+4,5,[],'line');
% set(h.data,'color','m')
% h=notBoxPlot(randn(50,10));
% set(h(5).data,'color','m')
%
% Example 5 - x and y are vectors
% clf
% x=[1,1,1,3,2,1,3,3,3,2,2,3,3];
% y=[7,8,6,1,5,7,2,1,3,4,5,2,4];
% notBoxPlot(y,x);
%
% Note: an alternative to the style used in Example 5 is to call
% notBoxPlot from a loop in an external function. In this case, the
% user will have to take care of the x-ticks and axis limits.
%
% Example 6 - replacing the SD with bars
% clf
% y=randn(50,1);
% clf
% notBoxPlot(y,1,[],'sdline')
% notBoxPlot(y,2)
% xlim([0,3])
%
%
% Rob Campbell - January 2010
%
% also see: boxplot

% Check input arguments
error(nargchk(0,4,nargin))
if nargin==0
    help(mfilename)
    return
end

if isvector(y), y=y(:); end

if nargin<2 || isempty(x)
    x=1:size(y,2);
end

if nargin<3 || isempty(jitter)
    jitter=0.3; %larger value means greater amplitude jitter
end

if nargin<4
  style='patch'; %Can also be 'line' or 'sdline'
end
style=lower(style);

if jitter==0 && strcmp(style,'patch')
    warning('A zero value for jitter means no patch object visible')
end

if isvector(y) & isvector(x) & length(x)>1
    x=x(:);

    if length(x)~=length(y)
        error('length(x) should equal length(y)')
    end

    u=unique(x);
    for ii=1:length(u)
        f=find(x==u(ii));
        h(ii)=notBoxPlot(y(f),u(ii),jitter,style);
    end

    %Make plot look pretty
    if length(u)>1
        xlim([min(u)-1,max(u)+1])
        set(gca,'XTick',u)
    end

    if nargout==1
        varargout{1}=h;
    end

    return

end

if length(x) ~= size(y,2)
    error('length of x doesn''t match the number of columns in y')
end

%We're going to render points with the same x value in different
%colors so we loop through all unique x values and do the plotting
%with nested functions. No clf in order to give the user more
%flexibility in combining plot elements.
hold on
[uX,a,b]=unique(x);

h=[];
for ii=1:length(uX)
    f=find(b==ii);
    h=[h,myPlotter(x(f),y(:,f))];
end

hold off

%Tidy up plot: make it look pretty
if length(x)>1
    set(gca,'XTick',unique(x))
    xlim([min(x)-1,max(x)+1])
end

if nargout==1
    varargout{1}=h;
end

%Nested functions follow

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function h=myPlotter(X,Y)

SEM=SEM_calc(Y); %Supplied external function
SD=nanstd(Y);  %Requires the stats toolbox
mu=nanmean(Y); %Requires the stats toolbox

%The plot colors to use for multiple sets of points on the same x
%location
cols=hsv(length(X)+1)*0.5;
cols(1,:)=0;
jitScale=jitter*0.55; %To scale the patch by the width of the jitter

for k=1:length(X)
    thisY=Y(:,k);
    thisY=thisY(~isnan(thisY));
    thisX=repmat(X(k),1,length(thisY));

    if strcmp(style,'patch')
      h(k).sdPtch=patchMaker(SD(k),[0.6,0.6,1]);
    end

    if strcmp(style,'patch') || strcmp(style,'sdline')
      h(k).semPtch=patchMaker(SEM(k),[1,0.6,0.6]);
      h(k).mu=plot([X(k)-jitScale,X(k)+jitScale],[mu(k),mu(k)],'-r',...
           'linewidth',2);
    end

    %Plot jittered raw data
    C=cols(k,:);
    J=(rand(size(thisX))-0.5)*jitter;

    h(k).data=plot(thisX+J, thisY, 'o', 'color', C,...
                   'markerfacecolor', C+(1-C)*0.65);
end

if strcmp(style,'line') | strcmp(style,'sdline')
  for k=1:length(X)
    %Plot SD
    h(k).sd=plot([X(k),X(k)],[mu(k)-SD(k),mu(k)+SD(k)],...
                 '-','color',[0.2,0.2,1],'linewidth',2);
    set(h(k).sd,'ZData',[1,1]*-1)
  end
end

if strcmp(style,'line')
    for k=1:length(X)
        %Plot mean and SEM
        h(k).mu=plot(X(k),mu(k),'o','color','r',...
            'markerfacecolor','r',...
            'markersize',10);

        h(k).sem=plot([X(k),X(k)],[mu(k)-SEM(k),mu(k)+SEM(k)],'-r',...
            'linewidth',2);
        h(k).xAxisLocation=x(k);
    end
end

function ptch=patchMaker(thisInterval,color)
    l=mu(k)-thisInterval;
    u=mu(k)+thisInterval;
    ptch=patch([X(k)-jitScale, X(k)+jitScale, X(k)+jitScale, X(k)-jitScale],...
           [l,l,u,u], 0);
    set(ptch,'edgecolor','none','facecolor',color)
end %function patchMaker

end %function myPlotter

end %function notBoxPlot