text:
classdef dmce_reserve_zone_mpc2 < mp.dmc_element % & mp.dmce_reserves
%MP.DMCE_RESERVE_ZONE_MPC2  Data model converter for reserve zones for MATPOWER case v2.

%   MATPOWER
%   Copyright (c) 2022, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%     properties
%     end     %% properties

    methods
        function name = name(obj)
            name = 'reserve_zone';
        end

        function df = data_field(obj)
            df = 'reserves';
        end

        function s = data_subs(obj)
            s = struct('type', {'.', '.'}, 'subs', {obj.data_field(), 'zones'});
        end

        function vmap = table_var_map(obj, dme, mpc)
            vmap = table_var_map@mp.dmc_element(obj, dme, mpc);

            import_req_fcn   = @(a, b, c, d)import_req(a, b, c, d);
            import_zones_fcn = @(a, b, c, d)import_zones(a, b, c, d);

            %% mapping for each name, default is {'col', []}
            vmap.uid        = {'IDs'};      %% consecutive IDs, starting at 1
            vmap.name       = {'cell', ''}; %% empty char
            vmap.status     = {'num', 1};   %% ones
            vmap.source_uid = {'cell', ''}; %% empty char
            vmap.req        = {'fcn', import_req_fcn};
            vmap.zones      = {'fcn', import_zones_fcn};
            vmap.prc        = {'num', 0};   %% zeros
        end

        function val = import_req(obj, mpc, spec, vn)
            val = mpc.reserves.req;
        end

        function val = import_zones(obj, mpc, spec, vn)
            val = mpc.reserves.zones;
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function contour_plot = plot_contour(date, pixel_dir, eddy_indexes, backgroundData, eddy_type)
%PLOT_CONTOUR plot contours for current eddies in the track viewer

temp = load([pixel_dir 'lat.mat']);
lat = temp.lat;
temp = load([pixel_dir 'lon.mat']);
lon = temp.lon;

contour_mask = false(length(lat), length(lon));

if strcmp(eddy_type, 'ant')
    pixel_file = [pixel_dir 'anticyclonic/pixels_' num2str(date)];
else
    pixel_file = [pixel_dir 'cyclonic/pixels_' num2str(date)];
end
temp = load(pixel_file);
pixel_list = temp.data;

for i = 1:length(eddy_indexes)
    contour_mask(pixel_list{eddy_indexes(i)}) = true;
end

contour_mask = bwperim(contour_mask);

contour = nan(length(lat), length(lon));

if strcmp(eddy_type, 'ant')
    if ~isempty(backgroundData)
        contour(contour_mask) = max(backgroundData(:)) + 0.2 * range(backgroundData(:));
    else
        contour(contour_mask) = 1;
    end
elseif strcmp(eddy_type, 'cyc')
    if ~isempty(backgroundData)
        contour(contour_mask) = min(backgroundData(:)) - 0.2 * range(backgroundData(:));
    else
        contour(contour_mask) = 0;
    end
else
    error('Not recognized eddy type');
end

contour_plot = pcolorm(lat, lon, contour);

end

author: jfaghm

repo: OceanEddies

save_path: github-repos/MATLAB/jfaghm-OceanEddies


text:
function linpack_c_test11 ( )

%*****************************************************************************80
%
%% TEST11 tests CHICO.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 June 2009
%
%  Author:
%
%    John Burkardt
%
  n = 3;
  lda = n;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST11\n' );
  fprintf ( 1, '  For a single precision complex (C)\n' );
  fprintf ( 1, '  Hermitian matrix (HI):\n' );
  fprintf ( 1, '  CHICO factors the matrix and estimates\n' );
  fprintf ( 1, '  the reciprocal condition number.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The matrix order is N = %d\n', n );
%
%  Set the values of the matrix A.
%
  seed = 123456789;

  for i = 1 : n
    [ a(i,i), seed ] = r4_uniform_01 ( seed );
    for j = i+1 : n
      [ a(i,j), seed ] = c4_uniform_01 ( seed );
      a(j,i) = conj ( a(i,j) );
    end
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The matrix A:\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    for j = 1 : n
      fprintf ( 1, '  (%8f  %8f)', real ( a(i,j) ), imag ( a(i,j) ) );
    end
    fprintf ( 1, '\n' );
  end
%
%  Factor the matrix A.
%
  [ a, ipvt, rcond ] = chico ( a, lda, n );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Estimated reciprocal condition RCOND = %f\n', rcond );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function y = RASTA(x,parameter)

y = rastafilt(x',parameter)';

author: singaxiong

repo: SignalGraph


text:
function glopts = loc2glob_3d ( cospitch, cosroll, cosyaw, sinpitch, sinroll, ...
  sinyaw, globas, locpts )

%*****************************************************************************80
%
%% LOC2GLOB_3D converts from a local to global coordinate system in 3D.
%
%  Discussion:
%
%    A global coordinate system is given.
%
%    A local coordinate system has been translated to the point with
%    global coordinates GLOBAS, and rotated through a yaw, a pitch, and
%    a roll.
%
%    A point has local coordinates LOCPTS, and it is desired to know
%    the point's global coordinates GLOPTS.
%
%    The transformation may be written as
%
%      GLOB = GLOBAS + N_YAW * N_PITCH * N_ROLL * LOC
%
%    where
%
%               (  cos(Yaw)   -sin(Yaw)        0      )
%    N_YAW    = (  sin(Yaw)    cos(Yaw)        0      )
%               (      0           0           1      )
%
%               (  cos(Pitch)      0       sin(Pitch) )
%    N_PITCH =  (      0           1           0      )
%               ( -sin(Pitch)      0       cos(Pitch) )
%
%               (      1           0           0      )
%    N_ROLL =   (      0       cos(Roll)  -sin(Roll)  )
%               (      0       sin(Roll)   cos(Roll)  )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real COSPITCH, COSROLL, COSYAW, the cosines of the
%    pitch, roll and yaw angles.
%
%    Input, real SINPITCH, SINROLL, SINYAW, the sines of the pitch,
%    roll and yaw angles.
%
%    Input, real GLOBAS(3,1), the global coordinates of the base
%    vector.
%
%    Input, real LOCPTS(3,1), the local coordinates of the point.
%
%    Output, real GLOPTS(3,1), the global coordinates of the point.
%
  glopts(1,1) = globas(1,1) + (  cosyaw * cospitch ) * locpts(1,1) ...
    + (  cosyaw * sinpitch * sinroll - sinyaw * cosroll ) * locpts(2,1) ...
    + (  cosyaw * sinpitch * cosroll + sinyaw * sinroll ) * locpts(3,1);

  glopts(2,1) = globas(2,1) + (  sinyaw * cospitch ) * locpts(1,1) ...
    + (  sinyaw * sinpitch * sinroll + cosyaw * cosroll ) * locpts(2,1) ...
    + (  sinyaw * sinpitch * cosroll - cosyaw * sinroll ) * locpts(3,1);

  glopts(3,1) = globas(3,1) + ( -sinpitch ) * locpts(1,1) ...
    + (  cospitch * sinroll ) * locpts(2,1) ...
    + (  cospitch * cosroll ) * locpts(3,1);

  return
end


text:
%% OFDM Syetem
% BER performance
% PAPR evaluation using Amplitude Clipping
%
clear all;close all;clc;
cdata_o=zeros(1,72);
%% Define Parameters
%--------------------------------------------------------------------------

N=input(' Number of Subcarriers N = ');
k=input(' Up samping factor k = ' );
n=input(' Number of all samples n = ');
M=input(' Mapping Order M = ');
P=input(' Phase Offset P = ');
%--------------------------------------------------------------------------
%Modems Creation
% ht=modem.qammod('M',M,'SymbolOrder','gray','PhaseOffset',P); % Tx - Modem
% hr=modem.qamdemod('M',M,'SymbolOrder','gray','PhaseOffset',P); % Rx - Modem
ht=modem.pskmod('M',M,'SymbolOrder','gray','PhaseOffset',P); % Tx - Modem
hr=modem.pskdemod('M',M,'SymbolOrder','gray','PhaseOffset',P); % Rx - Modem
%--------------------------------------------------------------------------
d=randi([0 M-1],1,n); % Data Generation
%--------------------------------------------------------------------------
% Generating the symbols before the mapping
[numberOfSymbols numberOfZeros symbols]=v2syms(d,N/k-16);

%--------------------------------------------------------------------------
% Mapping the data using the Tx-modem
mod_data=modulate(ht,symbols);
%--------------------------------------------------------------------------
%Zero Padding
for ct=1:size(mod_data,2)
    [z(ct) z_data(:,ct)]=zpp(mod_data(:,ct),N/k);
%
end
%Up sampling by the factor of k
up_data=rectpulse(z_data,k);
%--------------------------------------------------------------------------
% Basic Information displaying
disp('---------------------------------------------------------- ')
disp(['number of OFDM-Symbols = ' num2str(numberOfSymbols)])
disp(['each OFDM-Symbol consists of N = ' num2str(size(up_data,1))])
disp(['number of added zeros = ' num2str(numberOfZeros)])
%% --------------------------------------------------------------------------
% IFFT modulation
ifft_data=ifft(up_data,N);

%--------------------------------------------------------------------------
% find the PAPR
% par_o=zeros(size(ifft_data,2));
for a=1:size(ifft_data,2)
par_o(a)=papr2(ifft_data(:,a));

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

%% CCDF
iterations=length(par_o);
for ii=1:iterations
 % PAPR Original
        t=1;
    for zk=0:0.25:18-0.25
        if par_o(ii)>10^(zk/10)
            cdata_o(t)=cdata_o(t)+1;
        end
        t=t+1;
    end

end
cdata_o=cdata_o/iterations;
t=0:0.25:18-0.25;
%% Original Case
figure;
semilogy(t,cdata_o,'-ro','LineWidth',1,'MarkerEdgeColor','r','MarkerFaceColor','r','MarkerSize',5),grid;
hold on;
% xlim([4.5 13.5]);
xlabel('\gamma in dB');
ylabel('CCDF(\xi) = Pr(\xi >\gamma )');
legend('Original');
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
%------------------------------Adding AWGN---------------------------------
%--------------------------------------------------------------------------
% Adding AWGN
endsnr=40;
start=0;
inc=2;
lo=0;
for L=start:inc:endsnr
    sn=awgn(ifft_data,L,'measured');

    %----------------------------------------------------------------------
%     % FFT demodulation
    sfft=fft(sn,N);

    %----------------------------------------------------------------------
    % Down sample by the factor of k
    sd0=intdump(sfft,k);

    %----------------------------------------------------------------------
    % removing the zero padding
    for q=1:size(sd0,2)
        sd(:,q)=zppr(sd0(:,q),z(q));

    end

    % Demapping
    sdemap=demodulate(hr,sd);

    %----------------------------------------------------------------------
    % convert the received data from symbols to vector
    sv=syms2v(sdemap,numberOfZeros);

    % BER calculation
    lo=lo+1;
    [n0 r(lo)]=symerr(d,sv);

end
   snrv=start:inc:endsnr;
   figure;
semilogy(snrv,r,'-o'),grid;hold on;

legend('Original');
xlabel('SNR in dB');
ylabel('BER');
hold off;

% Power Spectral Density

fsMHz = 2;
[Pxx,W] = pwelch(ifft_data(:,10),[],[],2048,40);

figure;
plot([-1024:1023]*fsMHz/2048,10*log10(fftshift(Pxx)),'g');grid;hold on;

legend('Original');
xlabel('frequency, MHz')
ylabel('power spectral density')
hold off;
% Power Spectral Density after the SSPA
ifft_data_sspa=ifft_data./((1+ifft_data.^4).^(1/4));
figure;
[Pxx_sspa,W] = pwelch(ifft_data_sspa(:,10),[],[],2048,100);
plot([-1024:1023]*fsMHz/2048,10*log10(fftshift(Pxx_sspa)),'g');grid;hold on;

legend('Original');
xlabel('frequency, MHz')
ylabel('power spectral density')
hold off;

author: Sable

repo: mcbench-benchmarks


text:
function [rho,rhou,Ener] = Euler1D(rho, rhou, Ener, FinalTime)

% function [rho, rhou, Ener] = Euler1D(rho, rhou, Ener, FinalTime)
% Purpose  : Integrate 1D Euler equations until FinalTime starting with
%            initial conditions [rho, rhou, Ener]

Globals1D;

% Parameters
gamma = 1.4; CFL = 1.0; time = 0;

% Prepare for adaptive time stepping
mindx = min(x(2,:)-x(1,:));

% Limit initial solution
rho =SlopeLimitN(rho); rhou=SlopeLimitN(rhou); Ener=SlopeLimitN(Ener);

% outer time step loop
while(time<FinalTime)

  Temp = (Ener - 0.5*(rhou).^2./rho)./rho;
  cvel = sqrt(gamma*(gamma-1)*Temp);
  dt = CFL*min(min(mindx./(abs(rhou./rho)+cvel)));

  if(time+dt>FinalTime)
    dt = FinalTime-time;
  end

  % 3rd order SSP Runge-Kutta

  % SSP RK Stage 1.
  [rhsrho,rhsrhou,rhsEner]  = EulerRHS1D(rho, rhou, Ener);
  rho1  = rho  + dt*rhsrho;
  rhou1 = rhou + dt*rhsrhou;
  Ener1 = Ener + dt*rhsEner;

  % Limit fields
  rho1  = SlopeLimitN(rho1); rhou1 = SlopeLimitN(rhou1); Ener1 = SlopeLimitN(Ener1);

  % SSP RK Stage 2.
  [rhsrho,rhsrhou,rhsEner]  = EulerRHS1D(rho1, rhou1, Ener1);
  rho2   = (3*rho  + rho1  + dt*rhsrho )/4;
  rhou2  = (3*rhou + rhou1 + dt*rhsrhou)/4;
  Ener2  = (3*Ener + Ener1 + dt*rhsEner)/4;

  % Limit fields
  rho2  = SlopeLimitN(rho2); rhou2 = SlopeLimitN(rhou2); Ener2 = SlopeLimitN(Ener2);

  % SSP RK Stage 3.
  [rhsrho,rhsrhou,rhsEner]  = EulerRHS1D(rho2, rhou2, Ener2);
  rho  = (rho  + 2*rho2  + 2*dt*rhsrho )/3;
  rhou = (rhou + 2*rhou2 + 2*dt*rhsrhou)/3;
  Ener = (Ener + 2*Ener2 + 2*dt*rhsEner)/3;

  % Limit solution
  rho =SlopeLimitN(rho); rhou=SlopeLimitN(rhou); Ener=SlopeLimitN(Ener);

  % Increment time and adapt timestep
  time = time+dt;
end
return

repo: Aero-matlab


text:
function [cmds, errs, msgs] = gitSetBranch(repo, branch, preview, quiet)
cmds = {};

if ~exist('repo','var') || isempty(repo)
    repo = [pwd, '/'];
end
if ~exist('branch','var') || isempty(branch)
    branch = 'development';
end
if ~exist('preview','var')
    preview = false;
end
if ~exist('quiet','var') || isempty(quiet)
    quiet = 1;
end

currdir = pwd;

repoFull = filesepStandard_startup(repo,'full');
ii = 1;

cmds{ii,1} = sprintf('cd %s', repoFull); ii = ii+1;
cmds{ii,1} = sprintf('git checkout %s', branch); ii = ii+1;
cmds{ii,1} = sprintf('git pull'); ii = ii+1;

[errs, msgs] = exeShellCmds(cmds, preview, quiet);

% If requested branch was not checkout it was because it doesn't exist
% so we create it
branchnew = gitGetBranch(repo);
if ~strcmp(branchnew, branch)
    cmds{ii,1} = sprintf('git checkout -b %s', branch);
    [errs, msgs] = exeShellCmds(cmds, preview, quiet);
end

cd(currdir);

author: BUNPC


text:
classdef Tx < adi.AD9371.Base & adi.common.Tx
    % adi.AD9371.Tx Transmit data from the AD9371 transceiver
    %   The adi.AD9371.Tx System object is a signal sink that can tranmsit
    %   complex data from the AD9371.
    %
    %   tx = adi.AD9371.Tx;
    %   tx = adi.AD9371.Tx('uri','192.168.2.1');
    %
    %   AD9371 Datasheet

    properties
        %AttenuationChannel0 Attenuation Channel 0
        %   Attentuation specified as a scalar from -89.75 to 0 dB with a
        %   resolution of 0.25 dB.
        AttenuationChannel0 = -30;
        %AttenuationChannel1 Attenuation Channel 1
        %   Attentuation specified as a scalar from -89.75 to 0 dB with a
        %   resolution of 0.25 dB.
        AttenuationChannel1 = -30;
    end

    properties (Hidden, Nontunable, Access = protected)
        isOutput = true;
    end

    properties(Nontunable, Hidden, Constant)
        Type = 'Tx';
        channel_names = {'voltage0','voltage1','voltage2','voltage3'};
    end

    properties (Nontunable, Hidden)
        devName = 'axi-ad9371-tx-hpc';
    end

    methods
        %% Constructor
        function obj = Tx(varargin)
            coder.allowpcode('plain');
            obj = obj@adi.AD9371.Base(varargin{:});
        end
        % Check Attentuation
        function set.AttenuationChannel0(obj, value)
            validateattributes( value, { 'double','single' }, ...
                { 'real', 'scalar', 'finite', 'nonnan', 'nonempty', '>=', -89.75,'<=', 0}, ...
                '', 'Attenuation');
            assert(mod(value,1/4)==0, 'Attentuation must be a multiple of 0.25');
            obj.AttenuationChannel0 = value;
            if obj.ConnectedToDevice
                id = 'voltage0';
                obj.setAttributeLongLong(id,'hardwaregain',value,true);
            end
        end
        % Check Attentuation
        function set.AttenuationChannel1(obj, value)
            validateattributes( value, { 'double','single' }, ...
                { 'real', 'scalar', 'finite', 'nonnan', 'nonempty', '>=', -89.75,'<=', 0}, ...
                '', 'Attenuation');
            assert(mod(value,1/4)==0, 'Attentuation must be a multiple of 0.25');
            obj.AttenuationChannel1 = value;
            if obj.ConnectedToDevice
                id = 'voltage1';
                obj.setAttributeLongLong(id,'hardwaregain',value,true);
            end
        end

    end

    %% API Functions
    methods (Hidden, Access = protected)

        function setupInit(obj)
            % Write all attributes to device once connected through set
            % methods
            % Do writes directly to hardware without using set methods.
            % This is required sine Simulink support doesn't support
            % modification to nontunable variables at SetupImpl

            if obj.EnableCustomProfile
                writeProfileFile(obj);
            end

            id = sprintf('altvoltage%d',strcmp(obj.Type,'Tx'));
            obj.setAttributeLongLong(id,'TX_LO_frequency',obj.CenterFrequency ,true);
            obj.setAttributeLongLong('voltage0','hardwaregain',obj.AttenuationChannel0,true);
            obj.setAttributeLongLong('voltage1','hardwaregain',obj.AttenuationChannel1,true);

            obj.ToggleDDS(strcmp(obj.DataSource,'DDS'));
            if strcmp(obj.DataSource,'DDS')
                obj.DDSUpdate();
            end
        end

    end

end

author: analogdevicesinc


text:
function [ y, j ] = j_borrow_republican ( y, j )

%*****************************************************************************80
%
%% J_BORROW_REPUBLICAN borrows year-days from years in a Republican date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 June 2000
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input/output, integer Y, J, a YJ date.
%
  while ( j <= 0 )

    y = y - 1;

    days = year_length_republican ( y );

    j = j + days;

  end

  return
end


text:
function delete_count = file_sequence_delete ( file_name )

%*****************************************************************************80
%
%% FILE_SEQUENCE_DELETE deletes a file sequence.
%
%  Discussion:
%
%    We suppose the user has a set of files whose names differ only
%    in some numeric tag that is sequentially increasing, as, perhaps,
%    "file001.txt", "file002.txt" through "file137.txt", say.
%
%    The user specifies FILE_NAME as the name of the first file in the
%    sequence.  This function deletes that file, generates the next
%    name in the sequence, and, if a file with that name exists, it
%    deletes it as well.  The process continues until a file name is
%    reached for which there is no existing file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 October 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string FILE_NAME, the name of the first file in the sequence.
%
%    Output, integer DELETE_COUNT, the number of files deleted.
%
  delete_count = 0;

  while ( exist ( file_name, 'file' );

    delete ( file_name );

    delete_count = delete_count + 1;

    file_name = filename_inc ( file_name );

  end

  return
end


text:
function disk_grid_test01 ( )

%*****************************************************************************80
%
%% DISK_GRID_TEST01 tests DISK_GRID.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 November 2011
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST01:\n' );
  fprintf ( 1, '  DISK_GRID can define a grid of points\n' );
  fprintf ( 1, '  with N+1 points on a horizontal or vertical radius,\n' );
  fprintf ( 1, '  based on any disk.\n' );

  n = 20;
  r = 2.0;
  c = [ 1.0, 5.0 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, '  We use N = %d\n', n );
  fprintf ( 1, '  Radius R = %f\n', r );
  fprintf ( 1, '  Center C = (%f,%f)\n', c(1), c(2) );

  ng = disk_grid_count ( n, r, c );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of grid points will be %d\n', ng );

  cg = disk_grid ( n, r, c, ng );

  r82vec_print_part ( ng, cg, 20, '  Part of the grid point array:' );

  filename = 'disk_grid_test01.xy';

  r8mat_write ( filename, 2, ng, cg );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Data written to the file "%s".\n', filename );
%
% Plot the points.
%
  i = 0 : 50;
  t = 2.0 * pi * i / 50.0;
  cx = c(1);
  cx = cx + r .* cos ( t );
  cy = c(2);
  cy = cy + r .* sin ( t );

  figure
  clf
  hold on
  plot ( cx, cy, 'r-', 'Linewidth', 3 )
  plot ( cg(1,1:ng), cg(2,1:ng), 'b.' )
  hold off
  axis equal
  xlabel ( '<---X--->', 'Fontsize', 16 )
  ylabel ( '<---Y--->', 'Fontsize', 16 )
  title ( sprintf ( 'Circle grid with %d points', ng ), 'Fontsize', 16 );
  filename = 'disk_grid_test01.png';
  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Graphics data saved in file "%s"\n', filename );

  return
end


text:
function I = undistort(Idistorted, params)
fx = params.fx;
fy = params.fy;
cx = params.cx;
cy = params.cy;
k1 = params.k1;
k2 = params.k2;
p1 = params.p1;
p2 = params.p2;

Idistorted = single(Idistorted);

K = [fx 0 cx; 0 fy cy; 0 0 1];

I = zeros(size(Idistorted));
[i j] = find(~isnan(I));

% Xp = the xyz vals of points on the z plane
Xp = inv(K)*[j i ones(length(i),1)]';

% Now we calculate how those points distort i.e forward map them through the distortion
r2 = Xp(1,:).^2+Xp(2,:).^2;
x = Xp(1,:);
y = Xp(2,:);

x = x.*(1+k1*r2 + k2*r2.^2) + 2*p1.*x.*y + p2*(r2 + 2*x.^2);
y = y.*(1+k1*r2 + k2*r2.^2) + 2*p2.*x.*y + p1*(r2 + 2*y.^2);

% u and v are now the distorted cooridnates
u = reshape(fx*x + cx,size(I));
v = reshape(fy*y + cy,size(I));

% Now we perform a backward mapping in order to undistort the warped image coordinates
I = interp2(Idistorted, u, v);

end

repo: MSCKF


text:
function [texture structure] = structure_texture_decomposition_rof(im, theta, nIters, alp)
%
% Decompose the input IMAGE into structure and texture parts using the
% Rudin-Osher-Fatemi method. The final output is a linear combination
% of the decomposed texture and the structure parts.
%
% According to Wedel etal "An Improved Algorithm for TV-L1 optical flow"
%       equations (8)-(10)
%
% Test code
% [im1, im2, tu, tv] = read_image_flow_tune_para(4,0);
% [t1 s1] = structure_texture_decomposition_rof(im1);
% [t2 s2] = structure_texture_decomposition_rof(im2);
% indx = ~isnan(tu) | ~isnan(tv);
% uv = cat(3, tu, tv);
% uv(isnan(uv)) = 0;
% figure; imshow(-abs(partial_deriv(cat(3,im1, im2), uv)).*indx, []); title('original');
% figure; imshow(-abs(partial_deriv(cat(3,s1, s2), uv)).*indx, []);  title('structure');
% figure; imshow(-abs(partial_deriv(cat(3,t1, t2), uv)).*indx, []);  title('texture');
% tmp = partial_deriv(cat(3,t1, t2, uv);
% figure; imshow(t1, []); figure; imshow(s1, []);

%   Author: Deqing Sun, Department of Computer Science, Brown University
%   Contact: dqsun@cs.brown.edu
%   $Date: 2009 $
%
% Copyright 2009-2010, Brown University, Providence, RI. USA
%
%                          All Rights Reserved
%
% All commercial use of this software, whether direct or indirect, is
% strictly prohibited including, without limitation, incorporation into in
% a commercial product, use in a commercial service, or production of other
% artifacts for commercial purposes.
%
% Permission to use, copy, modify, and distribute this software and its
% documentation for research purposes is hereby granted without fee,
% provided that the above copyright notice appears in all copies and that
% both that copyright notice and this permission notice appear in
% supporting documentation, and that the name of the author and Brown
% University not be used in advertising or publicity pertaining to
% distribution of the software without specific, written prior permission.
%
% For commercial uses contact the Technology Venture Office of Brown University
%
% THE AUTHOR AND BROWN UNIVERSITY DISCLAIM ALL WARRANTIES WITH REGARD TO
% THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
% FITNESS FOR ANY PARTICULAR PURPOSE.  IN NO EVENT SHALL THE AUTHOR OR
% BROWN UNIVERSITY BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL
% DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
% PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
% ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
% THIS SOFTWARE.

if nargin == 1
    theta   = 1/8;
    nIters  = 100;
    alp     = 0.95;  % alp = 0.75 results in 4:1
end;

% Rescale the input image to [-1 1]
IM   = scale_image(im, -1,1);

% Backup orginal images
im   = IM;

% stepsize
delta = 1.0/(4.0*theta);

for iIm = 1:size(im,3)

    % Initialize dual variable p to be 0
    p = zeros([size(im,1) size(im,2) 2]);

    % Gradient descend
    I = squeeze(IM(:,:,iIm));

    for iter = 1:nIters

        % Compute divergence        eqn(8)
        div_p = imfilter(p(:,:,1), [-1 1 0], 'corr', 0)+ ...
                imfilter(p(:,:,2), [-1 1 0]', 'corr', 0);

        I_x = imfilter(I+theta*div_p, [-1 1], 'replicate');

        I_y = imfilter(I+theta*div_p, [-1 1]', 'replicate');

        % Update dual variable      eqn(9)
        p(:,:,1) = p(:,:,1) + delta*(I_x);
        p(:,:,2) = p(:,:,2) + delta*(I_y);

        % Reproject to |p| <= 1     eqn(10)
        reprojection = max(1.0, sqrt(p(:,:,1).^2 + p(:,:,2).^2));
        p(:,:,1) = p(:,:,1)./reprojection;
        p(:,:,2) = p(:,:,2)./reprojection;

    end

    % compute divergence
    div_p = imfilter(p(:,:,1), [-1 1 0], 'corr', 0)+ ...
            imfilter(p(:,:,2), [-1 1 0]', 'corr', 0);

    % compute structure component
    IM(:,:,iIm) = I + theta*div_p;

end;

texture   = squeeze(scale_image(im - alp*IM, 0, 255));

if nargout == 2
    structure = squeeze(scale_image(IM, 0, 255)); %(u-min(u(:)))/(max(u(:))-min(u(:))) - 1;
end;

author: kristinbranson

repo: JAABA


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Test ALM in the paper ''Fast Alternating Linearization Methods for
%       Minimizing the Sum of Two Convex Functions'', Donald Goldfarb,
%       Shiqian Ma and Katya Scheinberg, Tech. Report, Columbia University,
%       2009 - 2010.
%
% Author: Shiqian Ma
% Date  : Apr. 20, 2010
% IEOR, Columbia University, Copyright (2010)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

randn('state',0); rand('state',0);
% addpath('C:\Mywork\Optimization\work\code\ADM\Sparse_LowRank_Mat\Matdata');

% get data
tic;
dataformat = 'surveillance-video-Hall';
opts = getdata(dataformat);
time_getdata = toc;
fprintf('%f seconds to get data ! \n', time_getdata);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Call ALM to solve the problem
tic; out_ALM = ALM_SADAL_smoothed(opts.D,opts); time_ALM = toc;
fprintf('*******************************************************************\n');
%%%%%%%%% print stats %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% fprintf('*******************************************************************\n');
% fprintf('ALM  : iter: %d, relX: %3.2e, relY: %3.2e, StopCrit: %3.2e, time: %f\n', ...
%     out_ALM.iter, out_ALM.relX, out_ALM.relY, out_ALM.StopCrit, time_ALM);

% plot the figures

j = 10;
subplot(3,3,1); imshow(reshape(opts.D(:,j),144,176),[]);
subplot(3,3,2); imshow(reshape(out_ALM.X(:,j),144,176),[]);
subplot(3,3,3); imshow(reshape(out_ALM.Y(:,j),144,176),[]);

j = 80;
subplot(3,3,4); imshow(reshape(opts.D(:,j),144,176),[]);
subplot(3,3,5); imshow(reshape(out_ALM.X(:,j),144,176),[]);
subplot(3,3,6); imshow(reshape(out_ALM.Y(:,j),144,176),[]);

j = 150;
subplot(3,3,7); imshow(reshape(opts.D(:,j),144,176),[]);
subplot(3,3,8); imshow(reshape(out_ALM.X(:,j),144,176),[]);
subplot(3,3,9); imshow(reshape(out_ALM.Y(:,j),144,176),[]);

author: andrewssobral

repo: lrslibrary


text:
function [r shift]=rs_dl(slot_num,sym_num,port_num,n_id_cell,n_rb_dl,cp_type)

% [r shift]=rs_dl(slot_num,sym_num,port_num,n_id_cell,n_rb_dl,cp_type)
%
% This function supplies the reference symbols to be used on OFDM symbol
% 'sym_num' of slot 'slot_num' for antenna port 'port_num'.
%
% The reference symbols should be located on subcarriers (1+shift):6:n_rb_dl*12
%
% If a particular OFDM symbol of a particular port does not contain any
% reference symbols, this function will return the empty matrix in 'r'
% and NaN in 'shift'.
%
% 'cp_type' can be either 'normal' or 'extended'.

% Copyright 2012 Evrytania LLC (http://www.evrytania.com)
%
% Written by James Peroulas <james@evrytania.com>
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

% Constants
n_rb_maxdl=110;

error(nargchk(6,6,nargin));
error(chk_param(slot_num,'slot_num','scalar','real','integer','>=',0,'<=',19));
error(chk_param(sym_num,'sym_num','scalar','real','integer','>=',0,'<=',6));
error(chk_param(port_num,'port_num','scalar','real','integer','>=',0,'<=',3));
error(chk_param(n_id_cell,'n_id_cell','scalar','real','integer','>=',0,'<=',503));
error(chk_param(n_rb_dl,'n_rb_dl','scalar','real','integer','>=',0,'<=',110));
error(chk_param(cp_type,'cp_type','string'));

% Second level processing
if (strcmpi(cp_type,'normal'))
  n_symb_dl=7;
  n_cp=1;
elseif (strcmpi(cp_type,'extended'))
  n_symb_dl=6;
  n_cp=0;
else
  error('Unrecognized cp_type specified');
end

% Quick return in some cases
r=[];
shift=NaN;
if ((port_num==0)||(port_num==1))
  if ((sym_num~=0)&&(sym_num~=n_symb_dl-3))
    return
  end
end
if ((port_num==2)||(port_num==3))
  if (sym_num~=1)
    return
  end
end

% Create the source sequence.
c_init=2^10*(7*(slot_num+1)+sym_num+1)*(2*n_id_cell+1)+2*n_id_cell+n_cp;
c=lte_pn(c_init,4*n_rb_maxdl);

% Create symbols to be used if n_rb_maxdl RB's were used on the DL.
r_l_ns=(1/sqrt(2))*complex(1-2*c(1:2:end),1-2*c(2:2:end));

r=r_l_ns(1+n_rb_maxdl-n_rb_dl:2*n_rb_dl+n_rb_maxdl-n_rb_dl);

if ((port_num==0)&&(sym_num==0))
  v=0;
elseif ((port_num==0)&&(sym_num~=0))
  v=3;
elseif ((port_num==1)&&(sym_num==0))
  v=3;
elseif ((port_num==1)&&(sym_num~=0))
  v=0;
elseif (port_num==2)
  v=3*mod(slot_num,2);
elseif (port_num==3)
  v=3+3*mod(slot_num,2);
end

v_shift=mod(n_id_cell,6);

shift=mod(v+v_shift,6);

author: JiaoXianjun

repo: rtl-sdr-LTE


text:
% EEG_POINT2LAT - convert latency in data points to latency in ms relative
%                   to the time locking. Used in EEGLAB.
% Usage:
%       >> [newlat ] = eeg_point2lat( lat_array, [], srate);
%       >> [newlat ] = eeg_point2lat( lat_array, epoch_array,...
%                                 srate, timelimits, timeunit);
% Inputs:
%   lat_array   - latency array in data points assuming concatenated
%                 data epochs (see EEGLAB event structure)
%   epoch_array - epoch number corresponding to each latency value
%   srate       - data sampling rate in Hz
%   timelimits  - [min max] timelimits in 'timeunit' units (see below)
%   timeunit    - time unit in second. Default is 1 = seconds.
%
% Outputs:
%   newlat      - converted latency values (in 'timeunit' units) for each epoch
%
% Example:
%   tmpevent = EEG.event;
%   eeg_point2lat( [ tmpevent.latency ], [], EEG.srate, [EEG.xmin EEG.xmax]);
%   % returns the latency of all events in second for a continuous
%   % dataset EEG
%
%   eeg_point2lat( [ tmpevent.latency ], [ tmpevent.epoch ],
%                 EEG.srate, [EEG.xmin EEG.xmax]*1000, 1E-3);
%   % returns the latency of all events in millisecond for a dataset
%   % containing data epochs.
%
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2 Mai 2002
%
% See also: EEG_LAT2POINT, EEGLAB, POP_EDITIEVENTVALS, POP_LOADDAT

% Copyright (C) 2 Mai 2002 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. 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 HOLDER 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.

function newlat = eeg_point2lat( lat_array, epoch_array, srate, timewin, timeunit);

if nargin <3
    help eeg_point2lat;
    return;
end
if isempty( epoch_array )
    epoch_array = ones( size(lat_array) );
end
if nargin <4
    timewin = 0;
end
if nargin <5
	timeunit = 1;
end

if length(lat_array) ~= length(epoch_array)
	if length(epoch_array)~= 1
		disp('eeg_point2lat: latency and epoch arrays must have the same length'); return;
	else
		epoch_array = ones(1,length(lat_array))*epoch_array;
	end
end
if length(timewin) ~= 2
    disp('eeg_point2lat: timelimits array must have length 2'); return;
end
if iscell(epoch_array)
	epoch_array = [ epoch_array{:} ];
end
if iscell(lat_array)
	lat_array = [ lat_array{:} ];
end

timewin = timewin*timeunit;

if length(timewin) == 2
    pnts = (timewin(2)-timewin(1))*srate+1;
else
    pnts = 0;
end
newlat  = ((lat_array - (epoch_array-1)*pnts-1)/srate+timewin(1))/timeunit;
newlat = round(newlat*1E9)/1E9;

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function L = cmpndNoiseLogLikelihood(noise, mu, varsigma, y)

% CMPNDNOISELOGLIKELIHOOD Log likelihood of the data under the CMPND noise model.
% FORMAT
% DESC returns the log likelihood of a data set under the  compound noise model.
% ARG noise : the noise structure for which the log likelihood is required.
% ARG mu : input mean locations for the log likelihood.
% ARG varSigma : input variance locations for the log likelihood.
% ARG y : target locations for the log likelihood.
%
% SEEALSO : cmpndNoiseParamInit, cmpndNoiseLikelihood, noiseLogLikelihood
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005

% NOISE

L = 0;
for i = 1:length(noise.comp)
  L = L + noiseLogLikelihood(noise.comp{i}, ...
			     mu(:, i), ...
			     varsigma(:, i), ...
			     y(:, i));
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [n,m] = size(ad,dim)
% SIZE for adiff objects

[n,m] = size(ad.x);
if nargin==2 & dim==2, n=m; end
if nargout<2, n = [n,m]; end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


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

%*****************************************************************************80
%
%% GOMPERTZ_CDF_INV inverts the Gompertz CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Johnson, Kotz, and Balakrishnan,
%    Continuous Univariate Distributions, Volume 2, second edition,
%    Wiley, 1994, pages 25-26.
%
%  Parameters:
%
%    Input, real CDF, the value of the CDF.
%
%    Input, real A, B, the parameters of the PDF.
%    1 < A, 0 < B.
%
%    Output, real X, the corresponding argument.
%
  if ( cdf < 0.0 )
    x = 0.0;
  elseif ( cdf < 1.0 )
    x = log ( 1.0 - log ( 1.0 - cdf ) * log ( a ) / b  ) / log ( a );
  else
    x = r8_huge ( );
  end

  return
end


text:
function [t0, x0, u0] = shift(T, t0, x0, u,f)
% add noise to the control actions before applying it
con_cov = diag([0.005 deg2rad(2)]).^2;
con = u(1,:)' + sqrt(con_cov)*randn(2,1);
st = x0;

f_value = f(st,con);
st = st+ (T*f_value);

x0 = full(st);
t0 = t0 + T;

u0 = [u(2:size(u,1),:);u(size(u,1),:)]; % shift the control action
end

author: MMehrez


text:
function [J, grad] = lrCostFunction(theta, X, y, lambda)
%LRCOSTFUNCTION Compute cost and gradient for logistic regression with
%regularization
%   J = LRCOSTFUNCTION(theta, X, y, lambda) computes the cost of using
%   theta as the parameter for regularized logistic regression and the
%   gradient of the cost w.r.t. to the parameters.

% Initialize some useful values
m = length(y); % number of training examples

% You need to return the following variables correctly
J = 0;
grad = zeros(size(theta));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost of a particular choice of theta.
%               You should set J to the cost.
%               Compute the partial derivatives and set grad to the partial
%               derivatives of the cost w.r.t. each parameter in theta
%
% Hint: The computation of the cost function and gradients can be
%       efficiently vectorized. For example, consider the computation
%
%           sigmoid(X * theta)
%
%       Each row of the resulting matrix will contain the value of the
%       prediction for that example. You can make use of this to vectorize
%       the cost function and gradient computations.
%
% Hint: When computing the gradient of the regularized cost function,
%       there're many possible vectorized solutions, but one solution
%       looks like:
%           grad = (unregularized gradient for logistic regression)
%           temp = theta;
%           temp(1) = 0;   % because we don't add anything for j = 0
%           grad = grad + YOUR_CODE_HERE (using the temp variable)
%

tempTheta = theta;
tempTheta(1) = 0;

J = (-1 / m) * sum(y.*log(sigmoid(X * theta)) + (1 - y).*log(1 - sigmoid(X * theta))) + (lambda / (2 * m))*sum(tempTheta.^2);
temp = sigmoid (X * theta);
error = temp - y;
grad = (1 / m) * (X' * error) + (lambda/m)*tempTheta;

% =============================================================

grad = grad(:);

end

author: AvaisP


text:
classdef anisotropicdiffusion < TTeMPS_op_laplace
    % Class for anisotropic diffusion operator with tridiagonal diffusion matrix
    %
    %       [ 1 a 0     ...0 ]
    %       [ a 1 a 0    ..0 ]
    %   D = [ 0 a 1 a 0  ..0 ]
    %       [  ..   .. . .   ]
    %       [ 0 ... .. 0 a 1 ]
    %
    %

    %   TTeMPS Toolbox.
    %   Michael Steinlechner, 2013-2016
    %   Questions and contact: michael.steinlechner@epfl.ch
    %   BSD 2-clause license, see LICENSE.txt

    properties
        L
        D
        % precomputed spectral decomp of 1D Laplace:

    end

    methods

        function A = update_properties( A );

            A.rank = [1,  3*ones(1, length(A.U)-1), 1];  % the TT rank is always three for such Laplace-like tensors
            size_col_ = cellfun( @(y) size(y,1), A.U);
            A.size_col = size_col_ ./ (A.rank(1:end-1).*A.rank(2:end));
            A.size_row = cellfun( @(y) size(y,2), A.U);
            A.order = length( A.size_row );
        end
    end

    methods( Access = public )

        function A = anisotropicdiffusion( n, d, alpha )

            if ~exist('alpha', 'var')
                alpha = 0.25;
            end

            one = ones(n,1);
            q = linspace( -10, 10, n)';
            h = abs(q(2) - q(1));
            L = -spdiags( [one, -2*one, one], [-1 0 1], n, n) / (h^2);

            % superclass constructor
            A = A@TTeMPS_op_laplace( L, d );
            % precompute eigenvalue information and exponential for use in local
            A = initialize_precond( A );
            % preconditioner:
            A.L = L;

			[A.V_L, A.E_L] = eig(full(A.L));
            A.E_L = diag(A.E_L);

            A.D = spdiags( [-one,one], [-1,1], n, n ) / (2*h);
            I = speye( n, n );

            e1 = sparse( 1, 1, 1, 3, 1 );
            e2 = sparse( 2, 1, 1, 3, 1 );
            e3 = sparse( 3, 1, 1, 3, 1 );

            l_mid = sparse( 3, 1, 1, 9, 1 );                % e_3
            b_mid = sparse( 6, 1, 1, 9, 1 );                % e_6
            m_mid = sparse( [1;9], [1;1], [1;1], 9, 1 );    % e_1 + e_9
            c_mid = sparse( 2, 1, 1, 9, 1 );                % e_2

            A.U = cell( 1, d );
            A.U{1} = kron( A.L, e1 ) + kron( 2*alpha*A.D, e2 ) + kron( I, e3);
            A_mid = kron( A.L, l_mid ) + kron( 2*alpha*A.D, b_mid ) + kron( I, m_mid) + kron( A.D, c_mid );
            for i=2:d-1
                A.U{i} = A_mid;
            end
            A.U{d} = kron( I, e1 ) + kron( A.D, e2 ) + kron( A.L, e3);

            A = update_properties( A );

        end

        function expB = constr_precond_inner( A, X, mu )

            n = size(A.L, 1);
            sz = [X.rank(mu), X.size(mu), X.rank(mu+1)]

            B1 = zeros( X.rank(mu) );
            % calculate B1 part:
            for i = 1:mu-1
                % apply L to the i'th core
                tmp = X;
                Xi = matricize( tmp.U{i}, 2 );
                Xi = A.L*Xi;
                tmp.U{i} = tensorize( Xi, 2, [X.rank(i), n, X.rank(i+1)] );
                B1 = B1 + innerprod( X, tmp, 'LR', mu-1);
            end

            B3 = zeros( X.rank(mu+1) );
            % calculate B3 part:
            for i = mu+1:A.order
                tmp = X;
                Xi = matricize( tmp.U{i}, 2 );
                Xi = A.L*Xi;
                tmp.U{i} = tensorize( Xi, 2, [X.rank(i), n, X.rank(i+1)] );
                B3 = B3 + innerprod( X, tmp, 'RL', mu+1);
            end

            [V1,e1] = eig(B1);
            e1 = diag(e1);
            [V3,e3] = eig(B3);
            e3 = diag(e3);

            lmin = min(e1) + min(A.E_L) + min(e3);
            lmax = max(e1) + max(A.E_L) + max(e3);

            R = lmax/lmin

            [omega, alpha] = load_coefficients( R );

            k = 3;
            omega = omega/lmin;
            alpha = alpha/lmin;

            expB = cell(3,k);

            for i = 1:k
                expB{1,i} = omega(i) * V1*diag( exp( -alpha(i)*e1 ))*V1';    % include omega in first part
                expB{2,i} = A.V_L*diag( exp( -alpha(i)*A.E_L ))*A.V_L';
                expB{3,i} = V3*diag( exp( -alpha(i)*e3 ))*V3';
            end
        end

        function expB = constr_precond_outer( A, X, mu1, mu2 )

            n = size(A.L, 1);

            B1 = zeros( X.rank(mu1) );
            % calculate B1 part:
            for i = 1:mu1-1
                % apply L to the i'th core
                tmp = X;
                Xi = matricize( tmp.U{i}, 2 );
                Xi = A.L*Xi;
                tmp.U{i} = tensorize( Xi, 2, [X.rank(i), n, X.rank(i+1)] );
                B1 = B1 + innerprod( X, tmp, 'LR', mu1-1);
            end

            B3 = zeros( X.rank(mu2+1) );
            % calculate B3 part:
            for i = mu2+1:A.order
                tmp = X;
                Xi = matricize( tmp.U{i}, 2 );
                Xi = A.L*Xi;
                tmp.U{i} = tensorize( Xi, 2, [X.rank(i), n, X.rank(i+1)] );
                B3 = B3 + innerprod( X, tmp, 'RL', mu2+1);
            end

            [V1,e1] = eig(B1);
            e1 = diag(e1);
            [V3,e3] = eig(B3);
            e3 = diag(e3);

            lmin = min(e1) + 2*min(A.E_L) + min(e3);
            lmax = max(e1) + 2*max(A.E_L) + max(e3);

            R = lmax/lmin

            [omega, alpha] = load_coefficients( R );

            k = 3;
            omega = omega/lmin;
            alpha = alpha/lmin;

            expB = cell(4,k);

            for i = 1:k
                expB{1,i} = omega(i) * V1*diag( exp( -alpha(i)*e1 ))*V1';    % include omega in first part
                expB{2,i} = A.V_L*diag( exp( -alpha(i)*A.E_L ))*A.V_L';
                expB{3,i} = A.V_L*diag( exp( -alpha(i)*A.E_L ))*A.V_L';
                expB{4,i} = V3*diag( exp( -alpha(i)*e3 ))*V3';
            end
        end

        function P = constr_precond( A, k )

            d = A.order;

            lmin = d*min(A.E_L);
            lmax = d*max(A.E_L);

            R = lmax/lmin

            %  http://www.mis.mpg.de/scicomp/EXP_SUM/1_x/1_xk07_2E2
            %  0.0133615547183825570028305575534521842940   {omega[1]}
            %  0.0429728469424360175410925952177443321034   {omega[2]}
            %  0.1143029399081515586560726591147663100401   {omega[3]}
            %  0.2838881266934189482611071431161775535656   {omega[4]}
            %  0.6622322841999484042811198458711174907876   {omega[5]}
            %  1.4847175320092703810050463464342840325116   {omega[6]}
            %  3.4859753729916252771962870138366952232900   {omega[7]}
            %  0.0050213411684266507485648978019454613531   {alpha[1]}
            %  0.0312546410994290844202411500801774835168   {alpha[2]}
            %  0.1045970270084145620410366606112262388706   {alpha[3]}
            %  0.2920522758702768403556507270657505159761   {alpha[4]}
            %  0.7407504784499061527671195936939341208927   {alpha[5]}
            %  1.7609744335543204401530945069076494746696   {alpha[6]}
            %  4.0759036969145123916954953635638503328664   {alpha[7]}

            if k == 3
                [omega, alpha] = load_coefficients( R );

            elseif k == 7
                omega = [0.0133615547183825570028305575534521842940 0.0429728469424360175410925952177443321034 0.1143029399081515586560726591147663100401,...
                         0.2838881266934189482611071431161775535656 0.6622322841999484042811198458711174907876 1.4847175320092703810050463464342840325116,...
                         3.4859753729916252771962870138366952232900];
                alpha = [0.0050213411684266507485648978019454613531 0.0312546410994290844202411500801774835168 0.1045970270084145620410366606112262388706,...
                         0.2920522758702768403556507270657505159761 0.7407504784499061527671195936939341208927 1.7609744335543204401530945069076494746696,...
                         4.0759036969145123916954953635638503328664];
            else
                error('Unknown rank specified. Choose either k=3 or k=7');
            end

            omega = omega/lmin;
            alpha = alpha/lmin;

            % careful: all cores assumed to be of same size
            E = reshape( expm( -alpha(1) * A.L), [1, A.size_row(2), A.size_col(2), 1]);
            P = omega(1)*TTeMPS_op( repmat({E},1,d) );
            for i = 2:k
                E = reshape( expm( -alpha(i) * A.L), [1, A.size_row(2), A.size_col(2), 1]);
                P = P + omega(i)*TTeMPS_op( repmat({E},1,d) );
            end

        end

    end

end

author: NicolasBoumal

repo: manopt


text:
function X = multi_transpose (X)
% Multi-task cell array transpose.
%
%% INPUT
%   X: {n1 * n2} * t - input matrix
%
%% OUTPUT
%   X: {n2 * n1} * t - transpose every matrix in the cells.
%
%% LICENSE
%   This program is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%   Copyright (C) 2011 - 2012 Jiayu Zhou and Jieping Ye
%
%   You are suggested to first read the Manual.
%   For any problem, please contact with Jiayu Zhou via jiayu.zhou@asu.edu
%
%   Last modified on June 3, 2012.
%

for i = 1:length(X)
    X{i} = X{i}';
end

end

author: jiayuzhou

repo: MALSAR

save_path: github-repos/MATLAB/jiayuzhou-MALSAR


text:
function test_suite = test_clipRay
% Tests the function 'clipRay'
%   output = testClipRay(input)
%
%   Example
%   testClipRay
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2009-04-22,    using Matlab 7.7.0.471 (R2008b)
% Copyright 2009 INRA - Cepia Software Platform.
% Licensed under the terms of the LGPL, see the file "license.txt"

test_suite = functiontests(localfunctions);

function testRightSide(testCase) %#ok<*DEFNU>
% test edges totally inside window, possibly touching edges

% a basic bounding box
box = [0 100 0 100];

% inside
origin      = [30 40];
direction   = [10 0];
ray         = [origin direction];
expected    = [30 40 100 40];
testCase.assertEqual(expected, clipRay(ray, box), 'AbsTol', .01);

% outside
origin      = [30 140];
direction   = [10 0];
ray         = [origin direction];
testCase.assertTrue(sum(isnan(clipRay(ray, box)))==4);

% line inside, but ray outside
origin      = [130 40];
direction   = [10 0];
ray         = [origin direction];
testCase.assertTrue(sum(isnan(clipRay(ray, box)))==4);

function testLeftSide(testCase)
% test edges totally inside window, possibly touching edges

% a basic bounding box
box = [0 100 0 100];

% inside
origin      = [30 40];
direction   = [-10 0];
ray         = [origin direction];
expected    = [30 40 0 40];
testCase.assertEqual(expected, clipRay(ray, box), 'AbsTol', .01);

% outside
origin      = [30 140];
direction   = [-10 0];
ray         = [origin direction];
testCase.assertTrue(sum(isnan(clipRay(ray, box)))==4);

% line inside, but ray outside
origin      = [-30 40];
direction   = [-10 0];
ray         = [origin direction];
testCase.assertTrue(sum(isnan(clipRay(ray, box)))==4);

function testUpSide(testCase)
% test edges totally inside window, possibly touching edges

box = [0 100 0 100];

% inside
origin      = [30 40];
direction   = [0 10];
ray         = [origin direction];
expected    = [30 40 30 100];
testCase.assertEqual(expected, clipRay(ray, box), 'AbsTol', .01);

% outside
origin      = [130 40];
direction   = [0 10];
ray         = [origin direction];
testCase.assertTrue(sum(isnan(clipRay(ray, box)))==4);

% line inside, but ray outside
origin      = [30 140];
direction   = [0 10];
ray         = [origin direction];
testCase.assertTrue(sum(isnan(clipRay(ray, box)))==4);

function testDownSide(testCase)
% test edges totally inside window, possibly touching edges

box = [0 100 0 100];

% inside
origin      = [30 40];
direction   = [0 -10];
ray         = [origin direction];
expected    = [30 40 30 0];
testCase.assertEqual(expected, clipRay(ray, box), 'AbsTol', .01);

% outside
origin      = [130 40];
direction   = [0 -10];
ray         = [origin direction];
testCase.assertTrue(sum(isnan(clipRay(ray, box)))==4);

% line inside, but ray outside
origin      = [30 -40];
direction   = [0 -10];
ray         = [origin direction];
testCase.assertTrue(sum(isnan(clipRay(ray, box)))==4);

function testArray(testCase)
% test with an array of rays and a box

box = [0 100 0 100];

origins     = [30 40;30 40;30 140;130 40];
directions  = [10 0;0 10;10 0;0 10];
rays        = [origins directions];
expected    = [30 40 100 40;30 40 30 100;NaN NaN NaN NaN;NaN NaN NaN NaN];
clipped     = clipRay(rays, box);
testCase.assertEqual(expected, clipped, 'AbsTol', .01);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [l,u,vBest,currBest,lowerBound,isRefined,newInd]=bnbBranch(l,u,...
  vBest,currBest,lowerBound,isRefined)
% Perform the simple branching in brach an d bound
%
% Used in affine and metric Upgrade
%
% USAGE
%   [ l, ul, vBest, currBest, lowerBound, goodInd ]=bnbBranch( l, ...
%    u, vBest, currBest, lowerBound)
%
% INPUTS
%  l             - [ dim x nInterval ] lower bounds
%  u             - [ dim x nInterval ] upper bounds
%  vBest         - [ dim x nInterval ] best value of the parameter in the
%                                      interval
%  currBest      - [ 1 x nInterval ] best value of the criterion on the
%                                      interval
%  lowerBound    - [ 1 x nInterval ] underestimator of the criterion on the
%                                      interval
%  isRefined     - [ 1 x nInterval ] if true, the result in there has been
%                                      refined
%
% OUTPUTS
%  l             - [ dim x nInterval ] updated l
%  u             - [ dim x nInterval ] updated u
%  vBest         - [ dim x nInterval ] updated vBest
%  currBest      - [ 1 x nInterval ] updated currBest
%  lowerBound    - [ 1 x nInterval ] updated lowerBound
%  isRefined     - [ 1 x nInterval ] updated isRefined
%  newInd        - array of the new intervals in l,u,vBest ...
%
% EXAMPLE
%
% See also
%
% Vincent's Structure From Motion Toolbox      Version 3.1
% Copyright (C) 2008-2011 Vincent Rabaud.  [vrabaud-at-cs.ucsd.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the GPL [see external/gpl.txt]

% split the rectangle with the lowest lower bound and the one with the
% lowest current best: both have the best chances to be the right
% interval
[ disc, ind1 ]=min(lowerBound);
[ disc, ind2 ]=min(currBest);
% only do the one with the current best once in a while: if it is indeed
% the best, we are going to refine that interval for no reason
if rand()>0.9; goodInd=unique([ind1,ind1]);
else; goodInd=unique([ind1,ind1]); end

% also add the biggest interval
[ disc, ind3 ]=max(prod(u-l,1));
%  goodInd(end+1)=ind3;
%  goodInd(end+1)=1;

goodInd=unique(goodInd);

for ind=goodInd(end:-1:1)
  % define the boundaries of the new interval
  % one way of doing it is by simply splitting in the middle of the largest
  % dimension

  % you can also split where there is already a minimum, if it is not too
  % far from the center
  [ disc, dim ]=max(u(:,ind)-l(:,ind));
  if 0 && abs(vBest(dim,ind)-(u(dim,ind)+l(dim,ind))/2) < 0.4*(u(dim,ind)-l(dim,ind))/2
    l(:,end+1)=l(:,ind);
    u(:,end+1)=u(:,ind); u(dim,end)=vBest(dim,ind);
    l(:,end+1)=l(:,ind); l(dim,end)=vBest(dim,ind);
    u(:,end+1)=u(:,ind);
  else
    [ disc, dim ]=max(u(:,ind)-l(:,ind));
    l=l(:,[1:end,ind,ind]); u=u(:,[1:end,ind,ind]);
    u(dim,end-1)=(u(dim,ind)+l(dim,ind))/2;
    l(dim,end)=(u(dim,ind)+l(dim,ind))/2;
  end

  vBest(:,end+1:end+2)=[vBest(:,ind),vBest(:,ind)];

  % keep track of the current best and lower bound
  % those can be inherited from the interval we are splitting
  if vBest(dim,ind)>=l(dim,end)
    currBest(end+1:end+2)=[Inf,currBest(ind)];
    lowerBound(end+1:end+2)=[Inf,lowerBound(ind)];
    isRefined(end+1:end+2)=[false,isRefined(ind)];
  else
    currBest(end+1:end+2)=[currBest(ind),Inf];
    lowerBound(end+1:end+2)=[lowerBound(ind),Inf];
    isRefined(end+1:end+2)=[isRefined(ind),false];
  end
  % remove the interval we have split
  l(:,ind)=[]; u(:,ind)=[]; vBest(:,ind)=[]; currBest(ind)=[];
  lowerBound(ind)=[]; isRefined(ind)=[];
end

newInd=[size(l,2)-2*length(goodInd)+1:size(l,2)];

author: vrabaud


text:
function Anew = h_trid(A)
%  H_TRID(A) uses Householder method to form a tridiagonal matrix from A.
%  Must have a SQUARE SYMMETRIC matrix as the input.
%
%
%  Example:
%
%             B=[0 1 1;1 2 1;1 1 1];
%             h_trid(B)
%
%
% Author:  Matt Fig
% Contact: popkenai@yahoo.com

[M N] = size(A);
if M~=N || ~isequal(A,A')  %This just screens matricies that can't work.
   error('Matrix must be square symmetric only, see help.');
end

lngth = length(A);  % Preallocations.
v = zeros(lngth,1);
I = eye(lngth);
Aold = A;

for jj=1:lngth-2  % Build each vector j and run the whole procedure.
    v(1:jj) = 0;
    S = ss(Aold,jj);
    v(jj+1) = sqrt(.5*(1+abs(Aold(jj+1,jj))/(S+2*eps)));
    v(jj+2:lngth) = Aold(jj+2:lngth,jj)*sign(Aold(jj+1,jj))...
                   /(2*v(jj+1)*S+2*eps);
    P = I-2*v*v';
    Anew = P*Aold*P;
    Aold = Anew;
end

Anew(abs(Anew(:))<5e-14)=0; % Tolerence.

function anss = ss(A,jj)
% Subfunction for h_trid.
anss = sqrt(sum(A(jj+1:end,jj).^2));


text:
function gui_l (x, y)
%
% This file is needed by FASTICAG

% The load dialog for loading new data
% and new initial guess.

% @(#)$Id: gui_l.m,v 1.4 2004/07/27 13:09:26 jarmo Exp $

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Global variables

% Handle to the window
global hf_FastICA_Load;

% Handles to some of the controls in window
global he_FastICA_file;

% What is the load type of load dialog
global g_FastICA_loadType;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Configuration options
FIGURENAME = 'FastICA: Load';
FIGURETAG = 'f_FastICALoad';
FIGURESIZE = [x y 450 150];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Check to see if this figure is already open - it should not!
% Can't have more than one copy - otherwise the global
% variables and handles can get mixed up.
if ~isempty(findobj('Tag',FIGURETAG))
  error('Error: load dialog already open!');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize some of the controls' values

% What are we loading - who is calling?
caller = get(gcf, 'CurrentObject');

switch get(caller, 'Tag')
 case 'b_LoadData'                             % Do we load new data...
  loadString = 'Load data from variable in Matlab.';
  g_FastICA_loadType = 'data';
  FIGURENAME = 'FastICA: Load data';

 case 'b_LoadGuess'                            % ... or new initial guess?
  loadString = 'Load initial guess for mixing matrix A from variable in Matlab.';
  g_FastICA_loadType = 'guess';
  FIGURENAME = 'FastICA: Load initial guess';
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create the figure
a = figure('Color',[0.8 0.8 0.8], ...
	   'PaperType','a4letter', ...
	   'Name', FIGURENAME, ...
	   'NumberTitle', 'off', ...
	   'Tag', FIGURETAG, ...
	   'Position', FIGURESIZE, ...
	   'MenuBar', 'none');
set (a, 'Resize', 'off');

hf_FastICA_Load = a;

set(hf_FastICA_Load, 'HandleVisibility', 'callback');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% From here on it get's ugly as I have not had time to clean it up

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create the frames
pos_l=2;
pos_w=FIGURESIZE(3)-4;
pos_h=FIGURESIZE(4)-4;
pos_t=2;
h_f_load_background = uicontrol('Parent',a, ...
  'BackgroundColor',[0.701961 0.701961 0.701961], ...
  'Position',[pos_l pos_t pos_w pos_h], ...
  'Style','frame', ...
  'Tag','f_load_background');

pos_w=120;
pos_l=FIGURESIZE(3)-(pos_w+2+2);
pos_h=FIGURESIZE(4)-2*4;
pos_t=4;
h_f_load_side = uicontrol('Parent',a, ...
  'BackgroundColor',[0.701961 0.701961 0.701961], ...
  'Position',[pos_l pos_t pos_w pos_h], ...
  'Style','frame', ...
  'Tag','f_load_side');

pos_l=4;
pos_w=FIGURESIZE(3)-8-pos_w-2;
pos_h=FIGURESIZE(4)-8;
pos_t=4;
h_f_load = uicontrol('Parent',a, ...
  'BackgroundColor',[0.701961 0.701961 0.701961], ...
  'Position',[pos_l pos_t pos_w pos_h], ...
  'Style','frame', ...
  'Tag','f_load');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Controls in f_load
bgc = get(h_f_load, 'BackgroundColor');

pos_w=230;

pos_frame=get(h_f_load, 'Position');
pos_h = 40;
pos_t = pos_frame(2) + pos_frame(4) - pos_h - 6;
pos_l = pos_frame(1) + 6;

b = uicontrol('Parent',a, ...
  'BackgroundColor',bgc, ...
  'HorizontalAlignment','left', ...
  'Position',[pos_l pos_t pos_w pos_h], ...
  'String',loadString, ...
  'Style','text', ...
  'Tag','t_93');

pos_h = 20;
pos_t = pos_t - pos_h - 10;
pos_l = pos_frame(1) + 6;

b = uicontrol('Parent',a, ...
  'BackgroundColor',bgc, ...
  'HorizontalAlignment','left', ...
  'Position',[pos_l pos_t pos_w pos_h], ...
  'String','Name of the variable:', ...
  'Style','text', ...
  'Tag','t_92');

pos_w = 200;
pos_l = pos_l + 30;
pos_t = pos_t - pos_h;
he_FastICA_file = uicontrol('Parent',a, ...
  'BackgroundColor',[1 1 1], ...
  'HorizontalAlignment','left', ...
  'Position',[pos_l pos_t pos_w pos_h], ...
  'String','', ...
  'Style','edit', ...
  'Tag','e_file');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Controls in f_load_side
pos_vspace = 6;
pos_hspace = 10;
pos_frame = get(h_f_load_side, 'Position');
pos_w = 100;
pos_h = 30;
pos_l = pos_frame(1) + pos_hspace;
pos_t = pos_frame(2) + pos_frame(4) - pos_h - pos_vspace;
b = uicontrol('Parent',a, ...
  'BackgroundColor',[0.701961 0.701961 0.701961], ...
  'Callback','gui_lc Load', ...
  'Position',[pos_l pos_t pos_w pos_h], ...
  'String','Load', ...
  'Tag','b_lLoad');

pos_t=pos_t-pos_h-pos_vspace;
b = uicontrol('Parent',a, ...
  'BackgroundColor',[0.701961 0.701961 0.701961], ...
  'Callback','gui_lc Cancel', ...
  'Position',[pos_l pos_t pos_w pos_h], ...
  'String','Cancel', ...
  'Tag','b_lCancel');

pos_t = pos_frame(2) + pos_vspace;
b = uicontrol('Parent',a, ...
  'BackgroundColor',[0.701961 0.701961 0.701961], ...
  'Callback','gui_lc Help', ...
  'Position',[pos_l pos_t pos_w pos_h], ...
  'String','Help', ...
  'Tag','b_lHelp');

author: angeloyeo

repo: gongdols

save_path: github-repos/MATLAB/angeloyeo-gongdols


text:
function [FEVD, VAR] = VARfevd(VAR,VARopt)
% =======================================================================
% Compute FEVDs for a VAR model estimated with VARmodel. Three
% identification schemes can be specified: zero short-run restrictions,
% zero long run restrictions, and sign restrictions
% =======================================================================
% [FEVD, VAR] = VARfevd(VAR,VARopt)
% -----------------------------------------------------------------------
% INPUT
%   - VAR: structure, result of VARmodel function
%   - VARopt: options of the VAR (see VARopt from VARmodel)
% -----------------------------------------------------------------------
% OUTPUT
%	- FEVD(t,j,k): matrix with 't' steps, the FEVD due to 'j' shock for
%       'k' variable
%   - VAR: structure including VAR estimation results
%       * VAR.invA: identified contemporaneous A matrix
% =======================================================================
% Ambrogio Cesa Bianchi, March 2015
% ambrogio.cesabianchi@gmail.com

% I thank Dora Xia for pointing out a typo in the above description.

%% Check inputs
%===============================================
if ~exist('VAR','var')
    error('You need to provide VAR structure, result of VARmodel');
end
if ~exist('VARopt','var')
    error('You need to provide VAR options (VARopt from VARmodel)');
end

%% Retrieve and initialize variables
%===============================================
nsteps = VARopt.nsteps;
ident  = VARopt.ident;

S     = VAR.S;
Fcomp = VAR.Fcomp;
nvar  = VAR.nvar;
sigma = VAR.sigma;

MSE   = zeros(nvar,nvar,nsteps);
MSE_j = zeros(nvar,nvar,nsteps);
PSI   = zeros(nvar,nvar,nsteps);
FEVD  = zeros(nsteps,nvar,nvar);
SE    = zeros(nsteps,nvar);

%% Compute the multipliers
%===============================================
VARjunk = VAR;
VARjunk.sigma = eye(nvar);
IRFjunk = VARir(VARjunk,VARopt);

% this loop is to save the multipliers for each period
for mm = 1:nvar
    PSI(:,mm,:) = reshape(IRFjunk(:,:,mm)',1,nvar,nsteps);
end

%% Calculate the contribution to the MSE for each shock (i.e, FEVD)
%===============================================

for mm = 1:nvar % loop for the shocks

    % Calculate Total Mean Squared Error
    MSE(:,:,1) = sigma;

    for kk = 2:nsteps;
       MSE(:,:,kk) = MSE(:,:,kk-1) + PSI(:,:,kk)*sigma*PSI(:,:,kk)';
    end;

    % Get the matrix invA containing the structural impulses
    if strcmp(ident,'oir')
        [out, chol_flag] = chol(sigma);
        if chol_flag~=0; error('VCV is not positive definite'); end
        invA = out';
    elseif strcmp(ident,'bq')
        Finf_big = inv(eye(length(Fcomp))-Fcomp);   % from the companion
        Finf     = Finf_big(1:nvar,1:nvar);
        D        = chol(Finf*sigma*Finf')'; % identification: u2 has no effect on y1 in the long run
        invA = Finf\D;
    elseif strcmp(ident,'sr')
        [out, chol_flag] = chol(sigma);
        if chol_flag~=0; error('VCV is not positive definite'); end
        if isempty(S); error('Rotation matrix is not provided'); end
        invA = (out')*(S');
    end

    % Get the column of invA corresponding to the mm_th shock
    column = invA(:,mm);

    % Compute the mean square error
    MSE_j(:,:,1) = column*column';
    for kk = 2:nsteps
        MSE_j(:,:,kk) = MSE_j(:,:,kk-1) + PSI(:,:,kk)*(column*column')*PSI(:,:,kk)';
    end

    % Compute the Forecast Error Covariance Decomposition
    FECD = MSE_j./MSE;

    % Select only the variance terms
    for nn = 1:nsteps
        for ii = 1:nvar
            FEVD(nn,mm,ii) = FECD(ii,ii,nn);
            SE(nn,:) = sqrt( diag(MSE(:,:,nn))' );
        end
    end
end

% Update VARopt
VAR.invA = invA;

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox


text:
function exact = p38_exact ( )

%*****************************************************************************80
%
%% P38_EXACT returns the exact integral for problem 38.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real EXACT, the value of the integral.
%
  alpha = p38_param_get ();

  x = 2.0^alpha;

  exact = pi * besselj ( 0, x );

  return
end

author: johannesgerer


text:
function [influenceM] = getGlobalInfluenceM(IM, structsV)
%function [influenceM] = getGlobalInfluenceM(IM, structsV)
%   Get the global influence matrix (dose(:) = influenceM * weights(:)),
%   for all given structures.  The output influenceM is sparse, and has
%   rowNums equal to the number of voxels in the plan.
%   The row coord is the same as would be gotten by retrieving the
%   structure mask:
%          [mask3D, planC] =getUniformStr(Str_Num, planC, optS);
%
%JOD, 17 Nov 03.
%JRA, 27 Feb 04, cleanup.
%JRA, 24 Mar 04, Supports multiple structs, preallocate inflM, kicked up
%                the speed some.
%JRA, 14 Apr 04, Added monteCarlo MCflag.
%JOD, 19 Dec 05, comment mods.
%JJW, 20 Jun 06, MCflag removed; added call of getSingleGlobalInfluenceM
%APA, 10 Oct 06, updated to make compatible with changed location of beamlets field
%
%Usage:
%   function [gInfluenceM] = getGlobalInfluenceM(IM, structsV);
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

if length(structsV) == 1

    % this function is fast, but it works only for one structure at a time
    [influenceM] = getSingleGlobalInfluenceM(IM, structsV);

else

   global planC;

   beamlets = [IM.beams(:).beamlets];

   %get indices of structures stored under beamlets
   structIndV = getAssociatedStr({beamlets(:,1).strUID});

   if ~all(ismember(structsV,structIndV))
       influenceM = [];
       error('Dose not computed on some of the input structures. Please compute dose on and try again')
       return;
   end

   numPBs = size(beamlets,2);

   %Find minimum number of nonzero elements that need to be put in influenceM.
   for structNum = structsV;
       strBmletInd = find(structNum==structIndV);
        count = 0;
        for i=1:length(beamlets(strBmletInd,:))
            count = count + length(beamlets(strBmletInd,i).influence);
        end
        nnzV(structNum) = count;
   end
    maxnnz = max(nnzV);

    %Pre-initalize influence matrix, greatly speeds things up.
    indexS = planC{end};
    numVoxels = prod(getUniformScanSize(planC{indexS.scan}(getStructureAssociatedScan(structsV(1)))));
    influenceM = spalloc(numVoxels, numPBs, maxnnz);

    %Make sure structsV is a row vector, used in below for loop.
    if size(structsV, 1) ~= 1
        structsV = structsV';
    end

    %Loop over beamlets.
    for PBNum = 1 : numPBs

        %For each requested structure, add the effect of this beamlet to inflM.
        %*** Loops are in this order to cut down on out of order inserts into
        %*** sparse influence matrix. Greatly increases speed. Leave it!
        for structNum = structsV

            strBmletInd = find(structNum==structIndV);

            if ~isempty(beamlets(strBmletInd,PBNum).influence)

                doseV     = double(beamlets(strBmletInd,PBNum).influence);
                indV      = beamlets(strBmletInd,PBNum).indexV;
                maxVal    = beamlets(strBmletInd,PBNum).maxInfluenceVal;
                sizeParam = beamlets(strBmletInd,PBNum).fullLength;

                if isfield(beamlets, 'lowDosePoints')
                    lowDosePoints = unpackLogicals(beamlets(strBmletInd,PBNum).lowDosePoints, size(indV));
                    doseScaledV(~lowDosePoints) = doseV(~lowDosePoints) * (maxVal / (2^8 -1));
                    doseScaledV(lowDosePoints) = doseV(lowDosePoints) * (maxVal / (2^8 -1) / (2^8 -1));
                else
                    doseScaledV = doseV * (maxVal / (2^8 -1));
                end

                influenceM(indV,PBNum) = doseScaledV(:);
                doseScaledV = [];

            end

        end

    end

end % if length(structsV)==1

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni
%
%
% Problem 5- System response to sinusoidal input

w0=2*pi;
Hw0=(3*(j*w0)^2+4*j*w0+2)/((j*w0)^2+2*j*w0+1);
mag=abs(Hw0);
phas=angle(Hw0);
t=0:.1:5;
y=3*mag*cos(w0*t+pi/2+phas);
plot(t,y)
legend('Output signal y(t)')

author: Sable


text:
function bratu ( )

%*****************************************************************************80
%
%% BRATU uses BVP4C to solve the BRATU problem.
%
%  Discussion:
%
%    The Bratu equation includes a parameter lambda.  Depending on the
%    value of lambda, there may be 2, 1, or no solutions to the BVP.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 September 2013
%
%  Author:
%
%    Original MATLAB version by Shampine, Kierzenka, Reichelt.
%    This version by John Burkardt.
%
%  Reference:
%
%    Lawrence Shampine, Jacek Kierzenka, Mark Reichelt,
%    Solving boundary value problems for ordinary differential equations
%    in MATLAB with bvp4c.
%
  global lambda

  lambda_test = [ 0.45, 1.00, 3.50 ];

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BRATU:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Use BVP4C to solve the following boundary value problem:\n' );
  fprintf ( 1, '  y" + lambda * exp ( y ) = 0\n' );
  fprintf ( 1, '  y(0) = 0, y(1) = 0\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  When lambda = 1, there are two solutions.\n' );
  fprintf ( 1, '  Try\n' );
  fprintf ( 1, '    y(x) = 0.1, y''(x) = 0.\n' );
  fprintf ( 1, '  and\n' );
  fprintf ( 1, '    y(x) = 3.0, y''(x) = 0.0\n' );

  figure_num = 0;

  for test = 1 : 3
    lambda = lambda_test(test);
    figure_num = bratu_solver ( lambda, figure_num );
  end
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BRATU:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function figure_num = bratu_solver ( lambda, figure_num )

%*****************************************************************************80
%
%% BRATU_SOLVER calls BVP4C to solve the Bratu problem for a particular lambda.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 September 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real LAMBDA, the value of the parameter lambda.
%
%    Input/output, integer FIGURE_NUM, the number of figures being displayed.
%

%
%  Compute SOL1 with the first initial guess.
%
  x_init = linspace ( 0.0, 1.0, 5 );
  y_init = [ 0.1, 0.0 ];
  solinit = bvpinit ( x_init, y_init );
  sol1 = bvp4c ( @bratu_ode, @bratu_bc, solinit );
%
%  Compute SOL2 with the second initial guess.
%
  x_init = linspace ( 0.0, 1.0, 5 );
  y_init = [ 3.0, 0.0 ];
  solinit = bvpinit ( x_init, y_init );
  sol2 = bvp4c ( @bratu_ode, @bratu_bc, solinit );

  x = linspace ( 0.0, 1.0, 101 );
  y1 = deval ( sol1, x );
  y2 = deval ( sol2, x );
%
%  Display a plot of the two solutions.
%
  figure_num = figure_num + 1;
  figure ( figure_num )
  plot ( x, y1(1,:), 'r-', ...
         x, y2(1,:), 'g-', 'Linewidth', 2 );
  xlabel ( '<--- X --->' )
  ylabel ( '<--- Y --->' )
  title ( sprintf ( 'Bratu''s equation for \\lambda = %g\n', lambda ) );
  grid on
  filename = sprintf ( 'bratu_%f.png', lambda );
  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Saving plot file as "%s"\n', filename );

  return
end
function dydx = bratu_ode ( x, y )

%*****************************************************************************80
%
%% BRATU_ODE evaluates the right hand side of the ODE.
%
%  Discussion:
%
%    We assume that the differential equation has been rewritten as a
%    system of first order equations of the form
%
%      dydx = f(x,y)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 September 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the point at which the ODE is to be evaluated.
%
%    Input, real Y(M), the value of the solution at X.
%
%    Output, real DYDX(M), the value of the right hand side given X and Y.
%
  global lambda

  dydx(1) = y(2);
  dydx(2) = - lambda * exp ( y(1) );

  return
end
function bc = bratu_bc ( ya, yb )

%*****************************************************************************80
%
%% SAMPLE1_BC evaluates the boundary conditions.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 September 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real YA(M), YB(M), the solution value at the left and right endpoints.
%
%    Output, real BC(2), the value of the boundary conditions.
%
  bc(1) = ya(1);
  bc(2) = yb(1);

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end

repo: jburkardt-m


text:
function varargout = brainstorm( varargin )
% BRAINSTORM Brainstorm startup function.
%
% USAGE: brainstorm                 : Start Brainstorm
%        brainstorm start           : Start Brainstorm
%        brainstorm nogui           : Start Brainstorm with hidden interface (for scripts)
%        brainstorm server          : Start Brainstorm on a distant server (completely headless)
%        brainstorm [script] [args] : Start Brainstorm in server mode and execute the input script
%        brainstorm ... local       : Start Brainstorm with a local database (in .brainstorm folder)
%        brainstorm stop            : Quit Brainstorm
%        brainstorm reset           : Re-inialize Brainstorm (delete preferences and database)
%        brainstorm digitize        : Digitize points using a Polhemus system
%        brainstorm update          : Download and install latest Brainstorm update
%        brainstorm autopilot ...   : Call bst_autopilot with the following arguments
%        brainstorm setpath         : Add Brainstorm subdirectories to current path
%        brainstorm startjava       : Add Brainstorm Java classes to dynamic classpath
%        brainstorm info            : Open Brainstorm website
%        brainstorm license         : Displays license agreement window
%        brainstorm tutorial name   : Run the validation script attached to a tutorial (ctf, neuromag, raw, resting, yokogawa
%        brainstorm tutorial all    : Run all the validation scripts
%        brainstorm test            : Run a coverage test
%        brainstorm deploy          : Cleanup files and copy to git repository
%        brainstorm compile         : Compile Brainstorm with Matlab mcc compiler, including all plugins
%        brainstorm compile noplugs : Compile Brainstorm with Matlab mcc compiler, without the plugins
%        brainstorm workshop        : Download OpenMEEG and the SPM atlases, and run some small tests
%  res = brainstorm('status')       : Return brainstorm status (1=running, 0=stopped)

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

% Make sure that "more" is off
more off

% Compiled version
isCompiled = exist('isdeployed', 'builtin') && isdeployed;
if isCompiled
    BrainstormHomeDir = fileparts(fileparts(which(mfilename)));
    % Some versions of the compiler use different subfolders: search for the doc folder
    if ~exist(bst_fullfile(BrainstormHomeDir, 'doc'), 'file')
        docDir = file_find(BrainstormHomeDir, 'doc', 4, 1);
        if ~isempty(docDir)
            BrainstormHomeDir = fileparts(docDir);
        end
    end
else
    % Assume we are in the Brainstorm folder
    BrainstormHomeDir = fileparts(which(mfilename));
    % Add path to the core subfolders, used in this function
    corePath = fullfile(BrainstormHomeDir, 'toolbox', 'core');
    if ~exist(corePath, 'dir')
        error(['Unable to find ' corePath ' directory.']);
    end
    % Test if the folders are already in the path first
    p = path();
    if isempty(strfind(p, corePath))
        addpath(corePath, '-BEGIN');
        addpath(fullfile(BrainstormHomeDir, 'toolbox', 'io'), '-BEGIN');
        addpath(fullfile(BrainstormHomeDir, 'toolbox', 'misc'), '-BEGIN');
    end
end

% Set dynamic JAVA CLASS PATH
if ~exist('org.brainstorm.tree.BstNode', 'class')
    % Brainstorm jar file
    BstJar = fullfile(BrainstormHomeDir, 'java', 'brainstorm.jar');
    UpdateFile = fullfile(BrainstormHomeDir, 'java', 'outdated_jar.txt');
    % If there is a flag file indicating that brainstorm.jar should be deleted (see bst_update.m)
    if exist(UpdateFile, 'file')
        delete(BstJar);
        delete(UpdateFile);
        if exist(UpdateFile, 'file') || exist(BstJar, 'file')
            errordlg(['The following files could not be deleted automatically:' 10 BstJar 10 UpdateFile 10 10 ...
                      'An error occurred during the Brainstorm update.' 10 ...
                      'Delete this brainstorm3 folder and download manually' 10 ...
                      'a new one from the Brainstorm website.']);
        end
    end
    % Download brainstorm.jar if missing
    if ~exist(BstJar, 'file')
        disp('BST> Downloading brainstorm.jar...');
        bst_webread('https://github.com/brainstorm-tools/bst-java/raw/master/brainstorm/dist/brainstorm.jar', BstJar);
    end
    % Add Brainstorm JARs to classpath
    javaaddpath(fullfile(BrainstormHomeDir, 'java', 'RiverLayout.jar'));
    javaaddpath(BstJar);
end

% Default action : start
if (nargin == 0)
    action = 'start';
    BrainstormDbDir = [];
else
    action = lower(varargin{1});
    % Local start
    if ismember(action, {'start', 'nogui', 'server'}) && (nargin == 2) && strcmpi(varargin{2}, 'local')
        BrainstormDbDir = 'local';
    else
        BrainstormDbDir = [];
    end
end

res = 1;
switch action
    case 'start'
        bst_set_path(BrainstormHomeDir);
        bst_startup(BrainstormHomeDir, 1, BrainstormDbDir);
    case 'nogui'
        bst_set_path(BrainstormHomeDir);
        bst_startup(BrainstormHomeDir, 0, BrainstormDbDir);
    case 'server'
        bst_set_path(BrainstormHomeDir);
        bst_startup(BrainstormHomeDir, -1, BrainstormDbDir);
    case 'autopilot'
        if ~isappdata(0, 'BrainstormRunning')
            bst_set_path(BrainstormHomeDir);
            bst_startup(BrainstormHomeDir, 2, BrainstormDbDir);
        end
        res = bst_autopilot(varargin{2:end});
    case 'digitize'
        brainstorm nogui
        panel_digitize('Start');
    case {'status', 'isstarted', 'isrunning'}
        res = isappdata(0, 'BrainstormRunning');
    case {'exit', 'stop', 'quit'}
        bst_exit();
    case 'reset'
        bst_reset();
    case 'setpath'
        disp('Adding all Brainstorm directories to local path...');
        bst_set_path(BrainstormHomeDir);
    case 'startjava'
        disp('Starting Java...');
    case {'info', 'website'}
        web('https://neuroimage.usc.edu/brainstorm/', '-browser');
    case 'forum'
        web('https://neuroimage.usc.edu/forums/', '-browser');
    case 'license'
        bst_set_path(BrainstormHomeDir);
        bst_set('BrainstormHomeDir', BrainstormHomeDir);
        bst_license();
    case 'update'
        % Add path to java_dialog function
        addpath(fullfile(BrainstormHomeDir, 'toolbox', 'gui'));
        % Update
        bst_update(0);
    case 'tutorial'
        bst_set_path(BrainstormHomeDir);
        tutonames = varargin{2};
        % Tutorial folder
        if (nargin < 3)
            tutorial_dir = 'C:\Work\RawData\Tutorials';
        else
            tutorial_dir = varargin{3};
        end
        % Run all the tutorial scripts
        if isequal(tutonames, 'all')
            tutonames = {'ctf', 'raw', 'epilepsy', 'resting', 'neuromag', 'yokogawa', 'auditory'};
        elseif ischar(tutonames)
            tutonames = {tutonames};
        elseif ~iscell(tutonames)
            error('Invalid call.');
        end
        % Run all the requested validation scripts
        for i = 1:length(tutonames)
            disp([10 '===== TUTORIAL: ' upper(tutonames{i}) ' =====']);
            startTime = tic;
            if (length(varargin) == 4)
                eval(['tutorial_' lower(tutonames{i}) '(''' tutorial_dir ''', ''' varargin{4} ''');']);
            else
                eval(['tutorial_' lower(tutonames{i}) '(''' tutorial_dir ''');']);
            end
            % Done
            stopTime = toc(startTime);
            if (stopTime > 60)
                disp(sprintf('BST> Done in %dmin\n', round(stopTime/60)));
            else
                disp(sprintf('BST> Done in %ds\n', round(stopTime)));
            end
            % Close report viewer if it is not the last tutorial to run
            if (i < length(tutonames))
                bst_report('Close');
            end
        end

    case 'test'
        bst_set_path(BrainstormHomeDir);
        if (nargin < 2)
            error(['You must specify an empty test folder.' 10 'Usage: brainstorm test test_dir']);
        end
        test_dir = varargin{2};
        test_all(test_dir);

    case 'workshop'
        % Runs Brainstorm normally (asks for brainstorm_db)
        if ~isappdata(0, 'BrainstormRunning')
            bst_set_path(BrainstormHomeDir);
            bst_startup(BrainstormHomeDir, 1, BrainstormDbDir);
        end
        % Message
        java_dialog('msgbox', 'Brainstorm will now download additional files needed for the workshop.', 'Workshop');
        % Downloads OpenMEEG
        bst_plugin('Install', 'openmeeg', 1);
        % Downloads the TMP.nii SPM atlas
        bst_normalize_mni('install');
        % Message
        java_dialog('msgbox', ['Brainstorm will now test your display and open a 3D figure:' 10 10 ...
                               ' - You should see two surfaces: a brain surface and a transparent head.' 10 ...
                               ' - Make sure you can rotate the brain with your mouse, ' 10 ...
                               ' - Then close the figure.' 10 10], 'Workshop');
        % Creates an empty test protocol
        ProtocolName = 'TestWorkshop';
        gui_brainstorm('DeleteProtocol', ProtocolName);
        gui_brainstorm('CreateProtocol', ProtocolName, 0, 0);
        % Display the default anatomy cortex and head
        hFig = view_surface('@default_subject/tess_cortex_pial_low.mat');
        hFig = view_surface('@default_subject/tess_head.mat', [], [], hFig);
        waitfor(hFig);
        % Delete test protocol
        gui_brainstorm('DeleteProtocol', ProtocolName);
        % Confirmation message
        java_dialog('msgbox', 'You computer is ready for the workshop.', 'Workshop');

    case 'deploy'
        % Add path to deploy function
        bst_set_path(BrainstormHomeDir);
        addpath(fullfile(BrainstormHomeDir, 'deploy'));
        bst_set('BrainstormHomeDir', BrainstormHomeDir);
        % Deploy Braintorm
        bst_deploy();

    case 'compile'
        % Add path to deploy function
        bst_set_path(BrainstormHomeDir);
        addpath(fullfile(BrainstormHomeDir, 'deploy'));
        % Options
        if (nargin > 1)
            if strcmpi(varargin{2}, 'noplugs')
                isPlugs = 0;
            else
                error('Usage: brainstorm compile [noplugs]');
            end
        else
            isPlugs = 1;
        end
        % Matlab < 2020a: Old compilation function using deploytool
        if (bst_get('MatlabVersion') < 908)
            addpath(fullfile(BrainstormHomeDir, 'deploy', 'deprecated'));
            if isPlugs
                bst_deploy_java('2');
            else
                bst_deploy_java('1');
            end
        else
            bst_compile(isPlugs);
        end

    otherwise
        % Check if trying to execute a script
        if file_exist(varargin{1})
            ScriptFile = varargin{1};
        elseif file_exist(fullfile(pwd, varargin{1}))
            ScriptFile = fullfile(pwd, varargin{1});
        else
            ScriptFile = [];
        end
        % Execute script
        if ~isempty(ScriptFile)
            % Start brainstorm in server mode (local database or not)
            % With extra parameters
            if (length(varargin) > 1)
                params = varargin(2:end);
                if any(cellfun(@(c)isequal(c,'local'), varargin(2:end)))
                    brainstorm server local;
                    params(cellfun(@(c)isequal(c,'local'), params)) = [];
                else
                    brainstorm server;
                end
            % Without extra parameters
            else
                brainstorm server;
                params = [];
            end
            % Execute script
            if ~isempty(params)
                bst_call(@panel_command, 'ExecuteScript', ScriptFile, params{:});
            else
                bst_call(@panel_command, 'ExecuteScript', ScriptFile);
            end
            % Quit
            brainstorm stop;

        % Display usage
        else
            disp(' ');
            disp('Usage : brainstorm start           : Start Brainstorm');
            disp('        brainstorm nogui           : Start Brainstorm with hidden interface (for scripts)');
            disp('        brainstorm server          : Start Brainstorm on a distant server (completely headless)');
            disp('        brainstorm <script> <args> : Start Brainstorm in server mode, execute the input script and quit');
            disp('        brainstorm ... local       : Start Brainstorm with a local database (in .brainstorm folder)');
            disp('        brainstorm stop            : Quit Brainstorm');
            disp('        brainstorm update          : Download and install latest Brainstorm update (see bst_update)');
            disp('        brainstorm reset           : Re-initialize Brainstorm database and preferences');
            disp('        brainstorm digitize        : Digitize electrodes positions and head shape using a Polhemus system');
            disp('        brainstorm setpath         : Add Brainstorm subdirectories to current path');
            disp('        brainstorm startjava       : Add Brainstorm Java classes to dynamic classpath');
            disp('        brainstorm info            : Open Brainstorm website');
            disp('        brainstorm forum           : Open Brainstorm forum');
            disp('        brainstorm license         : Display license');
            disp('        brainstorm tutorial name   : Run the validation script attached to a tutorial (ctf, neuromag, raw, resting, yokogawa)');
            disp('        brainstorm tutorial all    : Run all the validation scripts');
            disp('        brainstorm packagebin      : Create separate zip files for all the currently available binary distributions');
            disp('  res = brainstorm(''status'')     : Return brainstorm status (1=running, 0=stopped)');
            disp(' ');
        end
end

% Return value
if (nargout >= 1)
    varargout{1} = res;
end

end

%% ===== SET PATH =====
function bst_set_path(BrainstormHomeDir)
    % Cancel add path in case of deployed application
    if bst_iscompiled()
        return
    end
    % Brainstorm folder itself
    addpath(BrainstormHomeDir, '-BEGIN'); % make sure the main brainstorm folder is in the path
    % List of folders to add
    NEXTDIR = {'external','toolbox'}; % in reverse order of priority
    for i = 1:length(NEXTDIR)
        nextdir = fullfile(BrainstormHomeDir,NEXTDIR{i});
        % Reset the last warning to blank
        lastwarn('');
        % Check that directory exist
        if ~isdir(nextdir)
            error(['Directory "' NEXTDIR{i} '" does not exist in Brainstorm path.' 10 ...
                   'Please re-install Brainstorm.']);
        end
        % Recursive search for subfolders in each main folder
        P = genpath(nextdir);
        % Add directory and subdirectories
        addpath(P, '-BEGIN');
    end
    % Adding user's mex path
    userMexDir = bst_get('UserMexDir');
    addpath(userMexDir, '-BEGIN');
    % Adding user's custom process path
    userProcessDir = bst_get('UserProcessDir');
    addpath(userProcessDir, '-BEGIN');
end

author: brainstorm-tools


text:
function nrgf = rgf_enum ( m )

%*****************************************************************************80
%
%% RGF_ENUM enumerates the restricted growth functions on M.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 January 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Donald Kreher, Douglas Simpson,
%    Combinatorial Algorithms,
%    CRC Press, 1998,
%    ISBN: 0-8493-3988-X,
%    LC: QA164.K73.
%
%  Parameters:
%
%    Input, integer M, the domain of the RGF is the integers
%    from 1 to M.  M must be positive.  However, for the enumeration routine
%    only, it is legal to call with any value of M.
%
%    Output, integer NRGF, the number of restricted growth
%    functions.
%
  if ( m < 0 )

    nrgf = 0;

  elseif ( m == 0 )

    nrgf = 1;

  else

    b = zeros ( m + 1, 1 );
    offset = 1;
    b(0+offset) = 1;
    for j = 1 : m
      b(j+offset) = 0;
      for i = 0 : j - 1
        b(j+offset) = b(j+offset) + i4_choose ( j - 1, i ) * b(i+offset);
      end
    end

    nrgf = b(m+offset);

  end

  return
end


text:
function imBW = segFloodFill(im, row, col)
%% Segment the image using flood fill, the seed locations is (row, col)

%% License
% ACADEMIC OR NON-PROFIT ORGANIZATION NONCOMMERCIAL RESEARCH USE ONLY
% Copyright (c) 2018 Bingyao Huang
% All rights reserved.

% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:

% The above copyright notice and this permission notice shall be included in all
% copies or substantial portions of the Software.

% If you publish results obtained using this software, please cite our paper.

% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
% AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
% OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
% SOFTWARE.

%%
imBW = zeros(size(im));
tol = 0.05;
thresh = 500;

% iterate if not enough flood fill was performed to avoid erasing the image
while(nnz(imBW) < thresh)
    % convert RGB to LAB
    imLab = rgb2lab(im);

    % normalize
    imLabNorm = sum((imLab - imLab(row,col,:)).^2, 3);
    imLabNorm = mat2gray(imLabNorm);

    tol = tol * 1.5; % tolerance
    imBW = grayconnected(imLabNorm, row, col, tol);
end
end

author: BingyaoHuang

repo: single-shot-pro-cam-calib


text:
function display(d)

%tstoolbox/@description/display
%   description/display
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

disp(['  Name : ' d.name]);
disp(['  Type : ' d.type]);
disp(' ')
disp(['  Attributes of data values : '])
disp(['  ' label(d.yunit) ' | '  d.yname])
disp(' ')
disp('  Comment : ')
display(d.comment);
disp(' ')
disp('  History : ')
display(d.history);

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [mv1, mv2, mv3] = specificIntMeanCurvDetails(img, varargin)
%SPECIFICINTMEANCURVDETAILS Ohser's Integral of Mean Curvature with details
%
%   this version is just for debugging.
%   It allows to extact contribution from each directions.
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 16/02/2005.
%

%   HISTORY

img = img~=0;

delta = [1 1 1];

% square faces of the unit cell
kr = [...
    1 2 3 4;...
    1 2 5 6;...
    1 3 5 7;...
    1 2 7 8;...
    3 5 4 6;...
    1 6 3 8;...
    2 4 5 7;...
    2 3 6 7;...
    1 5 3 8];

% triangle faces of the unit cell
kt = [...
    1 7 6;...
    2 5 8;...
    4 7 6;...
    3 5 8;...
    2 3 8;...
    4 1 7;...
    2 3 5;...
    4 1 6];

% unit surface for each cell
c1 = 0.045778;
c2 = 0.036981;
c3 = 0.035196;
c = [c1 c1 c1   c2 c2 c2   c2 c2 c2   c3 c3 c3 c3];

% distances between pixels
d1 = delta(1);
d2 = delta(2);
d3 = delta(3);
d12  = sqrt(delta(1)*delta(1) + delta(2)*delta(2));
d13  = sqrt(delta(1)*delta(1) + delta(3)*delta(3));
d23  = sqrt(delta(2)*delta(2) + delta(3)*delta(3));
s = (d12 + d23 + d13)/2;
a123 = 2*sqrt(s*(s-d12)*(s-d13)*(s-d23));

a = [d1*d2 d1*d3 d2*d3   d3*d12 d3*d12 d2*d13 d2*d13 d1*d23 d1*d23  a123 a123 a123 a123];

% compute gray-tone histogram of the image
h = grayHist(img);

mv1 = 0;
mv2 = 0;
mv3 = 0;

% for each type of configuration
for l=1:256

    v = l-1;
    b(1) = bitand(v,1)~=0;
    b(2) = bitand(v,2)~=0;
    b(3) = bitand(v,4)~=0;
    b(4) = bitand(v,8)~=0;
    b(5) = bitand(v,16)~=0;
    b(6) = bitand(v,32)~=0;
    b(7) = bitand(v,64)~=0;
    b(8) = bitand(v,128)~=0;

    % for each square face
    for nu=1:3

        b1 = b(kr(nu, 1));
        b2 = b(kr(nu, 2));
        b3 = b(kr(nu, 3));
        b4 = b(kr(nu, 4));

        s = sum([b1 b2 b3 b4]);
        if s==1
            mv1 = mv1 + h(l)*c(nu)/4/a(nu);
        elseif s==3
            mv1 = mv1 - h(l)*c(nu)/4/a(nu);
        end
    end

    % for each square diagonal face
    for nu = 4:9

        b1 = b(kr(nu, 1));
        b2 = b(kr(nu, 2));
        b3 = b(kr(nu, 3));
        b4 = b(kr(nu, 4));

        s = sum([b1 b2 b3 b4]);
        if s==1
            mv2 = mv2 + h(l)*c(nu)/4/a(nu);
        elseif s==3
            mv2 = mv2 - h(l)*c(nu)/4/a(nu);
        end
    end

    % for each triangular face
    for nu=10:13

        b1 = b(kt(nu-9, 1));
        b2 = b(kt(nu-9, 2));
        b3 = b(kt(nu-9, 3));
        b4 = b(kt(nu-5, 1));
        b5 = b(kt(nu-5, 2));
        b6 = b(kt(nu-5, 3));

        s1 = sum([b1 b2 b3]);
        s2 = sum([b4 b5 b6]);
        if s1==1
            mv3 = mv3 + h(l)*c(nu)/3/a(nu);
        end
        if s2==2
            mv3 = mv3 - h(l)*c(nu)/3/a(nu);
        end
    end
end

%imc = 4*pi*mv/sum(h(:));

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage