text:
function [collect] = ft_channelcombination(channelcmb, datachannel, includeauto, dirflag)

% FT_CHANNELCOMBINATION creates a cell-array with combinations of EEG/MEG
% channels for subsequent cross-spectral-density and coherence analysis
%
% You should specify channel combinations as a two-column cell-array,
%   cfg.channelcmb = {  'EMG' 'MLF31'
%                       'EMG' 'MLF32'
%                       'EMG' 'MLF33' };
% to compare EMG with these three sensors, or
%   cfg.channelcmb = { 'MEG' 'MEG' };
% to make all MEG combinations, or
%   cfg.channelcmb = { 'EMG' 'MEG' };
% to make all combinations between the EMG and all MEG channels.
%
% For each column, you can specify a mixture of real channel labels
% and of special strings that will be replaced by the corresponding
% channel labels. Channels that are not present in the raw datafile
% are automatically removed from the channel list.
%
% When directional connectivity measures will subsequently be computed, the
% interpretation of each channel-combination is that the direction of the
% interaction is from the first column to the second column.
%
% Note that the default behavior is to exclude symmetric pairs and
% auto-combinations.
%
% See also FT_CHANNELSELECTION

% Undocumented options:
%   includeauto = 0 (or 1), include auto-combinations
%   dirflag     = 0 (or 1, or 2) specifies the treatment of the order in
%                   the columns of the output. If dirflag = 0, the order is
%                   preserved, if dirflag = 1, the order is reversed,
%                   second input column is put first. If dirflag = 2, both
%                   directions are added to the list of combinations.

% Copyright (C) 2003-2015, Jan-Mathijs Schoffelen & Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

if nargin<3 || isempty(includeauto), includeauto = 0; end
if nargin<4 || isempty(dirflag),     dirflag     = 0; end

% ensure that datachannel is a column vector
datachannel = datachannel(:);

if ischar(channelcmb) && strcmp(channelcmb, 'all')
  % make all possible combinations of all channels
  channelcmb = {'all' 'all'};
end

% it should have a selection of two channels or channelgroups in each row
if size(channelcmb,1)==2 && size(channelcmb,2)~=2
  ft_warning('transposing channelcombination matrix');
  channelcmb = channelcmb';
end

% allow for channelcmb to be a 1x2 cell-array containing cells
if numel(channelcmb)==2 && iscell(channelcmb{1}) && iscell(channelcmb{2})
  channelcmb{1} = ft_channelselection(channelcmb{1}, datachannel);
  channelcmb{2} = ft_channelselection(channelcmb{2}, datachannel);
  n1  = numel(channelcmb{1});
  n2  = numel(channelcmb{2});
  tmp = cell(n1*n2+n1+n2,2);
  for k = 1:n1
    tmp((k-1)*n2+(1:n2), 1) = channelcmb{1}(k);
    tmp((k-1)*n2+(1:n2), 2) = channelcmb{2};
    tmp(n2*k+(1:n1),     1) = channelcmb{1};
    tmp(n2*k+(1:n1),     2) = channelcmb{1};
    tmp(n2*k+n1+(1:n2),  1) = channelcmb{2};
    tmp(n2*k+n1+(1:n2),  2) = channelcmb{2};
  end
  collect = tmp;
  return;
end

if isempty(setdiff(channelcmb(:), datachannel))
  % there is not much to do, since there are no channelgroups with special names
  % each element of the input therefore already contains a proper channel name

  switch dirflag
    case 0
      % nothing to do
    case 1
      % switch the order
      channelcmb = channelcmb(:,[2 1]);
    case 2
      channelcmb = [channelcmb; channelcmb(:,[2 1])];
    otherwise
      ft_error('unknown value for input argument ''dirflag''');
  end
  collect = channelcmb;

  if includeauto
    autochannel = unique(channelcmb(:));
    for ch=1:numel(autochannel)
      collect = cat(1, collect, [autochannel(ch) autochannel(ch)]);
    end
  end

else
  % a combination is made for each row of the input selection after
  % translating the channel group (such as 'all') to the proper channel names
  % and within each set, double occurences and autocombinations are removed

  selmat = false(numel(datachannel));
  for sel=1:size(channelcmb,1)
    % translate both columns and subsequently make all combinations
    channelcmb1 = ft_channelselection(channelcmb(sel,1), datachannel);
    channelcmb2 = ft_channelselection(channelcmb(sel,2), datachannel);

    % translate both columns and subsequently make all combinations
    list1 = match_str(datachannel, channelcmb1);
    list2 = match_str(datachannel, channelcmb2);

    selmat(list1,list2) = true;

    if ~includeauto,
      % exclude the auto-combinations
      selmat = selmat & ~eye(size(selmat));
    else
      % ensure that the appropriate diagonal entries are filled
      autovec = false(numel(datachannel),1);
      autovec([list1(:);list2(:)]) = true;
      selmat = selmat | diag(autovec);
    end

    if dirflag==2,
      % also fill the other 'direction'
      selmat(list2,list1) = true;
    end
  end

  if dirflag<2
    % remove double occurrences
    selmat   = selmat & ~tril(selmat, -1)';
  end
  [i1, i2] = find(selmat);

  switch dirflag
    case 0
      % original behavior,
      % row-to-column, i.e. outflow according to FT's convention
      indx = [i1 i2];
    case 1
      % column-to-row, i.e. inflow according to FT's convention
      indx = [i2 i1];
    case 2
      % both
      indx = [i1 i2];

    otherwise
      ft_error('unknown value for input argument ''dirflag''');
  end

  collect = [datachannel(indx(:,1)) datachannel(indx(:,2))];
end

author: spm


text:
function mono_total_enum_test ( )

%*****************************************************************************80
%
%% MONO_TOTAL_ENUM_TEST tests MONO_TOTAL_ENUM.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 September 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MONO_TOTAL_ENUM_TEST\n' );
  fprintf ( 1, '  MONO_TOTAL_ENUM can enumerate the number of monomials\n' );
  fprintf ( 1, '  in M variables, of total degree N.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '    N:' );
  for n = 0 : 8
    fprintf ( 1, '  %4d', n );
  end
  fprintf ( 1, '\n' );
  fprintf ( 1, '   M +------------------------------------------------------\n' );
  for m = 1 : 8
    fprintf ( 1, '  %2d |', m );
    for n = 0 : 8
      v = mono_total_enum ( m, n );
      fprintf ( 1, '  %4d', v );
    end
    fprintf ( 1, '\n' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function A = warmUpExercise()
%WARMUPEXERCISE Example function in octave
%   A = WARMUPEXERCISE() is an example function that returns the 5x5 identity matrix

A = [];
% ============= YOUR CODE HERE ==============
% Instructions: Return the 5x5 identity matrix
%               In octave, we return values by defining which variables
%               represent the return values (at the top of the file)
%               and then set them accordingly.
A = eye(5,5);

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

end

author: imLogM


text:
% std_topo() - uses topoplot() to get the interpolated Cartesian grid of the
%               specified component topo maps. The topo map grids are saved
%               into a (.icatopo) file and a pointer to the file is stored
%               in the EEG structure. If such a file already exists,
%               loads the information from it.
%
%               Returns the topo map grids of all the requested components. Also
%               returns the EEG sub-structure etc (i.e EEG.etc), which is modified
%               with a pointer to the float file and some information about the file.
% Usage:
%               >> X = std_topo(EEG, components, option);
%
%                  % Returns the ICA topo map grid for a dataset.
%                  % Updates the EEG structure in the Matlab environment and re-saves
% Inputs:
%   EEG        - an EEG dataset structure.
%   components - [numeric vector] components in the EEG structure to compute topo maps
%                      {default|[] -> all}
%   option     - ['gradient'|'laplacian'|'none'] compute gradient or laplacian of
%                the scale topography. This does not acffect the saved file which is
%                always 'none' {default is 'none' = the interpolated topo map}
% Optional inputs
%   'recompute'  - ['on'|'off'] force recomputing topo file even if it is
%                  already on disk.
%   'fileout'    - [string] Path of the folder to save output. The default
%                  is EEG.filepath
% Outputs:
%   X          - the topo map grid of the requested ICA components, each grid is
%                     one ROW of X.
%
% File output: [dataset_name].icatopo
%
% Authors:  Hilit Serby, Arnaud Delorme, SCCN, INC, UCSD, January, 2005
%
% See also  topoplot(), std_erp(), std_ersp(), std_spec(), std_preclust()

% Copyright (C) Hilit Serby, SCCN, INC, UCSD, October 11, 2004, hilit@sccn.ucsd.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function [X] = std_topo(EEG, comps, option, varargin)

if nargin < 1
    help std_topo;
    return;
end;
if isfield(EEG,'icaweights')
   numc = size(EEG.icaweights,1);
else
   error('EEG.icaweights not found');
end
if nargin < 2
   comps = 1:numc;
elseif isempty(comps)
   comps = 1:numc;
end

if nargin < 3
    option = 'none';
end;

g = finputcheck( varargin, { 'recompute'   'string'   { 'on','off' }   'off' ;...
                             'fileout'     'string'   []                EEG.filepath},...
                             'std_topo');
% if isstr(g), error(g); end;

% figure; toporeplot(grid,'style', 'both','plotrad', 0.5, 'intrad', 0.5, 'xsurface' ,Xi, 'ysurface',Yi );

% Topo information found in dataset
% ---------------------------------
if exist(fullfile(g.fileout, [ EEG.filename(1:end-3) 'icatopo' ])) && strcmpi(g.recompute, 'off')
    for k = 1:length(comps)
        tmp = std_readtopo( EEG, 1, comps(k));
        if strcmpi(option, 'gradient')
            [tmpx, tmpy]  = gradient(tmp); %Gradient
            tmp = [tmpx(:); tmpy(:)]';
        elseif strcmpi(option, 'laplacian')
            tmp = del2(tmp); %Laplacian
            tmp = tmp(:)';
        else
            tmp = tmp(:)';
        end;

        tmp = tmp(find(~isnan(tmp)));
        if k == 1
            X = zeros(length(comps),length(tmp)) ;
        end
        X(k,:) =  tmp;
    end
    return
end

all_topos = [];
for k = 1:numc

    % compute topo map grid (topoimage)
    % ---------------------------------
    chanlocs = EEG.chanlocs(EEG.icachansind);
    if isempty( [ chanlocs.theta ] )
        error('Channel locations are required for computing scalp topographies');
    end;
    [hfig grid plotrad Xi Yi] = topoplot( EEG.icawinv(:,k), chanlocs, ...
                                          'verbose', 'off',...
                                           'electrodes', 'on' ,'style','both',...
                                           'plotrad',0.55,'intrad',0.55,...
                                           'noplot', 'on', 'chaninfo', EEG.chaninfo);

    all_topos = setfield(all_topos, [ 'comp' int2str(k) '_grid' ], grid);
    all_topos = setfield(all_topos, [ 'comp' int2str(k) '_x' ]   , Xi(:,1));
    all_topos = setfield(all_topos, [ 'comp' int2str(k) '_y' ]   , Yi(:,1));

end

% Save topos in file
% ------------------
all_topos.datatype = 'TOPO';
tmpfile = fullfile( g.fileout, [ EEG.filename(1:end-3) 'icatopo' ]);
std_savedat(tmpfile, all_topos);

for k = 1:length(comps)
    tmp =  getfield(all_topos, [ 'comp' int2str(comps(k)) '_grid' ]);

    if strcmpi(option, 'gradient')
        [tmpx, tmpy]  = gradient(tmp); % Gradient
        tmp = [tmpx(:); tmpy(:)]';
    elseif strcmpi(option, 'laplacian')
        tmp = del2(tmp); % Laplacian
        tmp = tmp(:)';
    else
        tmp = tmp(:)';
    end;

    tmp = tmp(find(~isnan(tmp)));
    if k == 1
        X = zeros(length(comps),length(tmp)) ;
    end
    X(k,:) =  tmp;
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function net = alexnet(labelNum,lossWeight,opts)
channel_num=256;
mag=0.1;
netname='alexnet';

partRate=1;
textureRate=0;
output_num=labelNum;
net=load('../nets/imagenet-caffe-alex.mat');

net.layers=net.layers(1:end-7);
net.layers{end-1}.learningRate=[1,5];
net=add_block_mask(net,lossWeight, opts,'6',3,3,256,channel_num,1,1,partRate,textureRate,mag); %%%%%%%%%%%%%%%%%
net.layers{end+1} = struct('type', 'pool', 'name', 'pool6', ...
                           'method', 'max', ...
                           'pool', [3 3], ...
                           'stride', 2, ...
                           'pad', 0) ;
net=add_block_mask(net,lossWeight,opts,'7',6,6,channel_num,4096,1,0,partRate,textureRate,mag); %%%%%%%%%%%%%%%%%

net = add_dropout(net, opts, 'fc-7') ;
net = add_block(net, opts, 'fc-8', 1, 1, 4096, 4096, 1, 0) ;
net = add_dropout(net, opts, 'fc-8') ;

net = add_block(net, opts, '9', 1, 1, 4096, output_num, 1, 0) ;
net.layers(end) = [] ;
net.meta.classes.name={'pos'};
net.meta.classes.description={'pos'};
net.meta.normalization.imageSize=net.meta.normalization.imageSize(1:3);
net.meta.netname=netname;
end

repo: interpretableCNN


text:
function [varargout] = isequaln(varargin)

% ISEQUALN is a drop-in replacement for the same function that was
% introduced in MATLAB R2012a.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% try to automatically remove incorrect compat folders from the path, see https://github.com/fieldtrip/fieldtrip/issues/899

if isempty(strfind(mfilename('fullpath'), matlabroot))
  % automatic cleanup does not work when FieldTrip is installed inside the MATLAB path, see https://github.com/fieldtrip/fieldtrip/issues/1527

  alternatives = which(mfilename, '-all');
  if ~iscell(alternatives)
    % this is needed for octave, see https://github.com/fieldtrip/fieldtrip/pull/1171
    alternatives = {alternatives};
  end

  keep = true(size(alternatives));
  for i=1:numel(alternatives)
    keep(i) = keep(i) && ~any(alternatives{i}=='@');  % exclude methods from classes
    keep(i) = keep(i) && alternatives{i}(end)~='p';   % exclude precompiled files
  end
  alternatives = alternatives(keep);

  if exist(mfilename, 'builtin') || any(strncmp(alternatives, matlabroot, length(matlabroot)) & cellfun(@isempty, strfind(alternatives, fullfile('private', mfilename))))
    % remove this directory from the path
    p = fileparts(mfilename('fullpath'));
    warning('removing "%s" from your path, see http://bit.ly/2SPPjUS', p);
    rmpath(p);
    % call the original MATLAB function
    if exist(mfilename, 'builtin')
      [varargout{1:nargout}] = builtin(mfilename, varargin{:});
    else
      [varargout{1:nargout}] = feval(mfilename, varargin{:});
    end
    return
  end

end % automatic cleanup of compat directories

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% this is where the actual replacement code starts
% function tf = isequaln(varargin)

tf = isequalwithequalnans(varargin{:});

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% deal with the output arguments

varargout = {tf};

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function c = typeadj(a,TYPE)
%TYPEADJ      typecast of  a  to type TYPE
%
%   c = typeadj(a,TYPE)
%
%For details, see intval\typeof and intval\typeadj.
%

% written  10/16/98     S.M. Rump
% written  12/06/98     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
%

% gradient\@gradient\typeadj:  a  must be gradient (superior to intval)
  switch TYPE
    case 'gradient',       c = mid(a);
    case 'gradientintval', c = a;
  otherwise
    error('invalid type in call of typeadj')
  end

repo: OpenMAS


text:
% Figure 8.10: Approximate linear discrimination via linear programming
% Section 8.6.1, Boyd & Vandenberghe "Convex Optimization"
% Original by Lieven Vandenberghe
% Adapted for CVX by Joelle Skaf - 10/16/05
% (a figure is generated)
%
% The goal is to find a function f(x) = a'*x - b that classifies the non-
% separable points {x_1,...,x_N} and {y_1,...,y_M} by allowing some
% misclassification. a and b can be obtained by solving the following
% problem:
%           minimize    1'*u + 1'*v
%               s.t.    a'*x_i - b >= 1 - u_i        for i = 1,...,N
%                       a'*y_i - b <= -(1 - v_i)     for i = 1,...,M
%                       u >= 0 and v >= 0

% data generation
n = 2;
randn('state',2);
N = 50; M = 50;
Y = [1.5+0.9*randn(1,0.6*N), 1.5+0.7*randn(1,0.4*N);
     2*(randn(1,0.6*N)+1), 2*(randn(1,0.4*N)-1)];
X = [-1.5+0.9*randn(1,0.6*M),  -1.5+0.7*randn(1,0.4*M);
      2*(randn(1,0.6*M)-1), 2*(randn(1,0.4*M)+1)];
T = [-1 1; 1 1];
Y = T*Y;  X = T*X;

% Solution via CVX
cvx_begin
    variables a(n) b(1) u(N) v(M)
    minimize (ones(1,N)*u + ones(1,M)*v)
    X'*a - b >= 1 - u; %#ok
    Y'*a - b <= -(1 - v); %#ok
    u >= 0; %#ok
    v >= 0; %#ok
cvx_end

% Displaying results
linewidth = 0.5;  % for the squares and circles
t_min = min([X(1,:),Y(1,:)]);
t_max = max([X(1,:),Y(1,:)]);
tt = linspace(t_min-1,t_max+1,100);
p = -a(1)*tt/a(2) + b/a(2);
p1 = -a(1)*tt/a(2) + (b+1)/a(2);
p2 = -a(1)*tt/a(2) + (b-1)/a(2);

graph = plot(X(1,:),X(2,:), 'o', Y(1,:), Y(2,:), 'o');
set(graph(1),'LineWidth',linewidth);
set(graph(2),'LineWidth',linewidth);
set(graph(2),'MarkerFaceColor',[0 0.5 0]);
hold on;
plot(tt,p, '-r', tt,p1, '--r', tt,p2, '--r');
axis equal
title('Approximate linear discrimination via linear programming');
% print -deps svc-discr.eps

author: yu-jiang

repo: radpbook

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


text:
n = 2;
clear robot
robot.NB = n;
robot.parent = [0:n-1];
robot.pitch = zeros(1,n);  % revolute joints

robot.Xtree(1) = Twist(inv(SE3.Rx(pi/2)));  % joint 1 at origin, z in horizontal plane
robot.Xtree(2) = Twist(inv(SE3([1 0 0])));

robot.I(1) = SpatialInertia( 1, [0.5 0 0], diag([0 0 0]) );
robot.I(2) = SpatialInertia( 1, [0.5 0 0], diag([0 0 0]) );

q = [0 0]; qd = [0 0]; qdd = [0 0];

tau = ID( robot, q, qd, qdd)'

author: petercorke

repo: spatialmath-matlab


text:
function om_constants

% astrodynamic and utility constants

% Orbital Mechanics with MATLAB

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

global dtr rtd mu mmu smu omega req flat j2 aunit

dtr = pi / 180.0;

rtd = 180.0 / pi;

% earth gravitational constant (km**3/sec**2)

mu = 398600.436233;

% moon gravitational constant (km**3/sec**2)

mmu = 4902.800076;

% sun gravitational constant (km**3/sec**2)

smu = 132712440040.944;

% earth inertial rotation rate (radians/second)

omega = 7.292115486e-5;

% earth equatorial radius (kilometers)

req = 6378.1363;

% earth flattening factor (non-dimensional)

flat = 1.0 / 298.257;

% earth oblateness gravity coefficient (non-dimensional)

j2 = 0.00108263;

% astronomical unit (kilometers)

aunit = 149597870.691;

author: Sable

repo: mcbench-benchmarks


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: Setting up advection-diffusion solver
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [C,laplacian_C] = please_Update_Adv_Diff_Concentration_Flux_Limiter_FV(C,dt,dx,dy,uX,uY,k)

% C:     concentration
% dt:    time-step
% dx,dy: spatial steps in x and y, respectively
% uX:    x-Component of Velocity
% uY:    y-Component of Velocity
% k:     diffusion coefficient

% Compute Fluxes (Note: these calculations could be parallalized)
selection = 'superbee';
Fx = give_Necessary_Fluxes(C,dx,uX,'x',selection,dt); % Fluxes in x
Fy = give_Necessary_Fluxes(C,dy,uY,'y',selection,dt); % Fluxes in y

% "forward difference of fluxes in x"
F2x = [Fx(:,2:end) Fx(:,1)];
F1x = [Fx(:,end) Fx(:,1:end-1)];
diffX = 0.5/dx*( F2x - F1x );

% "forward differences of fluxes in y"
F2y = [Fy(2:end,:); Fy(1,:)];
F1y = [Fy(end,:); Fy(1:end-1,:)];
diffY = 0.5/dy*( F2y - F1y );

% Compute 2nd Derivative Terms
Cxx = DD(C,dx,'x');
Cyy = DD(C,dy,'y');

% Forms Laplacian
laplacian_C = Cxx+Cyy;

% UPWIND
C = C - dt * ( diffX + diffY ) + dt*( k*laplacian_C );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: Computes derivative based on sign of Velocity, u, using UPWIND
% approach
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function C_z = give_Necessary_Fluxes(C,dz,uZ,string,selection,dt)

C_z = zeros(size(C));
len = length(uZ(:,1));
signs = sign(uZ);

if strcmp(string,'x')

    %For periodicity on ends w/ UPWIND
    for i=1:len

        %left side of grid
        if signs(i,1) <= 0
            r = ( C(i,2) - C(i,1) ) / ( C(i,1) - C(i,len) );
            phi = please_Give_Flux_Limiter(r,selection);
            C_z(i,1) =  uZ(i,1)*C(i,2)   + 0.5*abs( uZ(i,1) )*( 1 - abs( 0.5*dt*uZ(i,1)/dz ) )*phi*( C(i,2) - C(i,len) );
        else
            r = ( C(i,1) - C(i,len) ) / ( C(i,2) - C(i,1) );
            phi = please_Give_Flux_Limiter(r,selection);
            C_z(i,1) =  uZ(i,1)*C(i,len) + 0.5*abs( uZ(i,1) )*( 1 - abs( 0.5*dt*uZ(i,1)/dz ) )*phi*( C(i,2) - C(i,len) );
        end

        %right side of grid
        if signs(len,1) <= 0
            r = ( C(i,1) - C(i,len) ) / ( C(i,len) - C(i,len-1) );
            phi = please_Give_Flux_Limiter(r,selection);
            C_z(i,len) =  uZ(i,len)*C(i,1)     + 0.5*abs( uZ(i,len) )*( 1 - 0.5*abs( dt*uZ(i,len)/dz ) )*phi*( C(i,1) - C(i,len-1) );
        else
            r = ( C(i,len) - C(i,len-1) ) / ( C(i,1) - C(i,len) );
            phi = please_Give_Flux_Limiter(r,selection);
            C_z(i,len) =  uZ(i,len)*C(i,len-1) + 0.5*abs( uZ(i,len) )*( 1 - 0.5*abs( dt*uZ(i,len)/dz ) )*phi*( C(i,1) - C(i,len-1) );
        end

    end
    %Standard Upwind
    for i=1:len
        for j=2:len-1
            if signs(i,j) <= 0
                r = ( C(i,j+1) - C(i,j) ) / ( C(i,j) - C(i,j-1) );
                phi = please_Give_Flux_Limiter(r,selection);
                C_z(i,j) = uZ(i,j)*C(i,j+1) + 0.5*abs( uZ(i,j) )*( 1 - 0.5*abs( dt*uZ(i,j)/dz ) )*phi*( C(i,j+1) - C(i,j-1) );
            else
                r = ( C(i,j) - C(i,j-1) ) / ( C(i,j+1) - C(i,j) );
                phi = please_Give_Flux_Limiter(r,selection);
                C_z(i,j) = uZ(i,j)*C(i,j-1) + 0.5*abs( uZ(i,j) )*( 1 - 0.5*abs( dt*uZ(i,j)/dz ) )*phi*( C(i,j+1) - C(i,j-1) );
            end
        end
    end

    % Ends x-Direction calculation %

elseif strcmp(string,'y')

    %For periodicity on ends w/ UPWIND
    for i=1:len

        %bottom of grid
        if signs(1,i) <= 0
            r = ( C(2,i) - C(1,i) ) / ( C(1,i) - C(len,i) );
            phi = please_Give_Flux_Limiter(r,selection);
            C_z(1,i) =  uZ(1,i)*C(2,i)   + 0.5*abs( uZ(1,i) )*( 1 - abs( 0.5*dt*uZ(1,i)/dz ) )*phi*( C(2,i) - C(len,i) );
        else
            r = ( C(1,i) - C(len,i) ) / ( C(2,i) - C(1,i) );
            phi = please_Give_Flux_Limiter(r,selection);
            C_z(1,i) =  uZ(1,i)*C(len,i) + 0.5*abs( uZ(1,i) )*( 1 - abs( 0.5*dt*uZ(1,i)/dz ) )*phi*( C(2,i) - C(len,i) );
        end

        %top of grid
        if signs(len,1) <= 0
            r = ( C(1,i) - C(len,i) ) / ( C(len,i) - C(len-1,i) );
            phi = please_Give_Flux_Limiter(r,selection);
            C_z(len,i) =  uZ(len,i)*C(1,i)       + 0.5*abs( uZ(len,i) )*( 1 - abs( 0.5*dt*uZ(len,i)/dz ) )*phi*( C(1,i) - C(len-1,i) );
        else
            r = ( C(len,i) - C(len-1,i) ) / ( C(1,i) - C(len,i) );
            phi = please_Give_Flux_Limiter(r,selection);
            C_z(len,i) =  uZ(len,i)*C(len-1,i)   + 0.5*abs( uZ(len,i) )*( 1 - abs( 0.5*dt*uZ(len,i)/dz ) )*phi*( C(1,i) - C(len-1,i) );
        end

    end

    %Standard Upwind
    for i=1:len
        for j=2:len-1
            if signs(j,i) <= 0
                r = ( C(j+1,i) - C(j,i) ) / ( C(j,i) - C(j-1,i) );
                phi = please_Give_Flux_Limiter(r,selection);
                C_z(j,i) = uZ(j,i)*C(j+1,i) + 0.5*abs( uZ(j,i) )*( 1 - abs( 0.5*dt*uZ(j,i)/dz ) )*phi*( C(j+1,i) - C(j-1,i) );
            else
                r = ( C(j,i) - C(j-1,i) ) / ( C(j+1,i) - C(j,i) );
                phi = please_Give_Flux_Limiter(r,selection);
                C_z(j,i) = uZ(j,i)*C(j-1,i) + 0.5*abs( uZ(j,i) )*( 1 - abs( 0.5*dt*uZ(j,i)/dz ) )*phi*( C(j+1,i) - C(j-1,i) );
            end
        end
    end

    % Ends y-Direction calculation %

else

    fprintf('\n\n\n ERROR IN FUNCTION FOR COMPUTING FLUX LIMITERS\n');

end

clear signs; clear len;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: computes flux limiter with choice of which one
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function phi = please_Give_Flux_Limiter(r,selection)

if strcmp(selection,'superbee')
    max1 = max( min(1,2*r),min(2,r) );
    phi = max(0,max1);
else
    fprintf('\n\n');
    error('NEED TO CHOOSE AN APPROPRIATE FLUX LIMITER');
end

author: nickabattista


text:
function y=logsum(x,d)
%LOGSUM logsum(x,d)=log(sum(exp(x),d))
%  d gives dimension to sum along

%      Copyright (C) Mike Brookes 1998
%
%      Last modified Mon Oct 12 15:47:25 1998
%
%   VOICEBOX is a MATLAB toolbox for speech processing. Home page is at
%   http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%

if nargin==1
   d=[find(size(x)-1) 1];
   d=d(1);
end
n=size(x,d);
if n<=1, y=x; return; end
s=size(x);
p=[d:ndims(x) 1:d-1];
z=reshape(permute(x,p),n,prod(s)/n);

y=max(z);
y=y+log(sum(exp(z-y(ones(n,1),:))));

s(d)=1;
y=ipermute(reshape(y,s(p)),p);


text:
function sphere_grid_test07 ( )

%*****************************************************************************80
%
%% SPHERE_GRID_TEST07 tests SPHERE_GRID_Q4.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 August 2010
%
%  Author:
%
%    John Burkardt
%
  lat_num = 3;
  long_num = 4;
  rectangle_num = lat_num * long_num;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPHERE_GRID_TEST07\n' );
  fprintf ( 1, '  SPHERE_GRID_Q4 computes a grid\n' );
  fprintf ( 1, '  of Q4 rectangular elements on a sphere in 3D.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of latitudes is      %d\n', lat_num );
  fprintf ( 1, '  Number of longitudes is     %d\n', long_num );
  fprintf ( 1, '  The number of rectangles is %d\n', rectangle_num );

  rectangle_node = sphere_grid_q4 ( lat_num, long_num );

  i4mat_transpose_print ( 4, rectangle_num, rectangle_node, ...
    '  Rectangle vertices:' );

  return
end


text:
function [] = visBinaryCRBMLearning(RBM);
%------------------------------------------------------------------
%  visBinaryCRBMLearning(RBM);
%------------------------------------------------------------------
% DES

figure(99)
colormap gray
set(99,'name','Learning Convolutional RBM (Binary Visible)');
subplot(221);
imagesc(RBM.auxVars.batchX'); axis image; axis off;
title('Visible Data');

subplot(222);
imagesc(RBM.eVis'); axis image; axis off;
title('Reconstruction');

[c,r,k]=size(RBM.eHid);
subplot(223);
data = reshape(RBM.eHid,r*c,k);
visWeights(data,1);
title('Feature Maps');

[c,r,k]=size(RBM.W);
subplot(224);
data = reshape(RBM.W,r*c,k);
visWeights(data,1);
title('Learned Filters');
drawnow

figure(98)
colormap gray;
set(98,'name','Learning Convolutional RBM (Binary Visible)');
subplot(221);
bar(1:RBM.nFM,RBM.c);  axis square
xlim([1 RBM.nFM]);axis square;
xlabel('Feature Index')
title('Hidden Biases')

subplot(222);
data = squeeze(mean(mean(RBM.eHid0)));
bar(1:RBM.nFM,data);
xlim([1 RBM.nFM]);axis square;
xlabel('Feature Index')
title(sprintf('Mean Hidden Activation\nactual=%g \n target = %g\n',mean(data),RBM.sparsity))
drawnow

subplot(223);
data = reshape(RBM.dW,r*c,k);
visWeights(data); axis image; axis off
title('Weight Gradients')

subplot(224);
bar(1:RBM.nFM,RBM.auxVars.dcSparse); axis square;
xlim([1 RBM.nFM]);axis square;
xlabel('Hidden Feature Index')
title('Sparsenss Offset')
%  title('Hidden Bias Gradients')

author: dustinstansbury

repo: medal


text:
% Copyright (C) 1993-2013, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for MATLAB (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

%%begin

% Create a white square
im = testpattern('squares', 256, 256, 128);
% and rotate it
im = irotate(im, -0.3);
idisp(im)

% find the edges
edges = icanny(im);
idisp(edges)

% compute the Hough transform of the edge image
h = Hough(edges)
about h
% where we have a Hough object

% We can view the Hough accumulator array
figure(2); h.show()
% where each point represents a line and the bright spots represent
% dominant lines in the scene

% We can find those dominant lines
lines = h.lines();
about lines
% which is an array of LineFeature objects

% Each LineFeature object has a number of properties
lines(1)
% which we can access individually by
lines(1).theta
lines(1).rho

% If we display these
lines
% we see that

axis([-1.4 -1.1 -190 -110])

% We can avoid this problem by using non-local minima suppression
h = Hough(edges, 'suppress', 5)
% in this case with a radius of 5 Hough cells

% Repeating the process above we see fewer lines
lines = h.lines()

figure(1); idisp(im)
h.plot('b')
im = iread('church.png', 'grey', 'double');
idisp(im)
edges = icanny(im);
idisp(edges)
h = Hough(edges, 'suppress', 5)
h.show()
h.lines()
lines = h.lines()
idisp(im)
lines(1:10).plot();

repo: machinevision-toolbox-matlab


text:
function blkStruct = slblocks
% This function specifies that the library 'mylib'
% should appear in the Library Browser with the
% name 'My Library'

    Browser.Library = 'FMT_Model_Lib';
    % 'mylib' is the name of the library

    Browser.Name = 'FMT Toolbox';
    % 'My Library' is the library name that appears
    % in the Library Browser

    blkStruct.Browser = Browser;

author: Firmament-Autopilot

repo: FMT-Model


text:
% ---------------------------------------------------------
% Copyright (c) 2017, Gurkirt Singh
% This code and is available
% under the terms of MIT License provided in LICENSE.
% Please retain this notice and LICENSE if you use
% this file (or any portion of it) in your project.
% ---------------------------------------------------------

function opts = initDatasetOpts(data_root,baseDir,dataset,imgType,model_type,listid,iteration_num,iouthresh,costtype,gap)

opts = struct();
opts.imgType = imgType;
opts.costtype = costtype;
opts.gap = gap;
opts.baseDir = baseDir;
opts.imgType = imgType;
opts.dataset = dataset;
opts.iouThresh = iouthresh;
opts.weight = iteration_num;
opts.listid = listid;

testlist = ['display',listid];
% opts.vidList = sprintf('%s/%s/splitfiles/%s.txt',data_root,dataset,testlist);
opts.vidList = sprintf('%s.txt',testlist);

if strcmp(dataset,'ucf24')
    opts.actions = {'Basketball','BasketballDunk','Biking','CliffDiving','CricketBowling',...
        'Diving','Fencing','FloorGymnastics','GolfSwing','HorseRiding','IceDancing',...
        'LongJump','PoleVault','RopeClimbing','SalsaSpin','SkateBoarding','Skiing',...
        'Skijet','SoccerJuggling','Surfing','TennisSwing','TrampolineJumping',...
        'VolleyballSpiking','WalkingWithDog'};
elseif strcmp(dataset,'JHMDB')
    opts.actions = {'brush_hair','catch','clap','climb_stairs','golf','jump',...
        'kick_ball','pick','pour','pullup','push','run','shoot_ball','shoot_bow',...
        'shoot_gun','sit','stand','swing_baseball','throw','walk','wave'};
elseif strcmp(dataset,'LIRIS')
    opts.actions = {'discussion', 'give_object_to_person','put_take_obj_into_from_box_desk',...
        'enter_leave_room_no_unlocking','try_enter_room_unsuccessfully','unlock_enter_leave_room',...
        'leave_baggage_unattended','handshaking','typing_on_keyboard','telephone_conversation'};
end

opts.imgDir = sprintf('%s/%s/%s-images/',data_root,dataset,imgType);

opts.detDir = sprintf('%s/%s/detections/%s-%s-%s-%06d/',baseDir,dataset,model_type,imgType,listid,iteration_num);
opts.annotFile = sprintf('%s/%s/splitfiles/annots.mat',data_root,dataset);

opts.actPathDir = sprintf('%s/%s/actionPaths/%s-%s-%s-%06d-%s-%d-%04d/',baseDir,dataset,model_type,imgType,listid,iteration_num,costtype,gap,iouthresh*100);
opts.tubeDir = sprintf('%s/%s/actionTubes/%s-%s-%s-%06d-%s-%d-%04d/',baseDir,dataset,model_type,imgType,listid,iteration_num,costtype,gap,iouthresh*100);

if exist(opts.detDir,'dir')
    if ~isdir(opts.actPathDir)
        fprintf('Creating %s\n',opts.actPathDir);
        mkdir(opts.actPathDir)
    end
    if ~isdir(opts.tubeDir)
        mkdir(opts.tubeDir)
    end
    if strcmp(dataset,'ucf24') || strcmp(dataset,'JHMDB')
        createdires({opts.actPathDir},opts.actions)
    end
end

author: gurkirt

repo: realtime-action-detection


text:
function Index = SVM_AFdetection_withTrainingModel(training_feature,training_class,test_feature,test_class)
%******************************************************
% $ This function is used for calculating the SVM-based AF detection
% indices. You need to train the SVM model first using the training_feature
% and training_class, and then using the test_feature and test_class for calculating the AF
% detection incices
%
% $ Reference:
% Q. Li, C. Y. Liu, J. Oster and G. D. Clifford. Chapter: Signal processing
% and feature selection preprocessing for classification in noisy healthcare data. In book: Machine Learning for Healthcare Technologies, Edition: 1st, Publisher: IET, Editors: David A. Clifton, 2016.
%
% $ Variable declaration:
% Input:
% training_feature: the feature metrix from AF detection, you need to run
% AF_feature function first to obtain the AF features, each row is a
% feature vector from an RR interval time series
% training_class: the labels of AF for the trained RR ineterval time series, N*1
% metrix, 1 for AF and 0 for non-AF
% test_feature: the feature metrix from AF detection, you need to run
% AF_feature function first to obtain the AF features, each row is a
% feature vector from an RR interval time series
% test_class: the labels of AF for the tested RR ineterval time series, N*1
% metrix, 1 for AF and 0 for non-AF
% Output:
% Index: AF detection indices output, with the defined order of indices of TP FN FP TN Se Sp Acc PPV NPV J in turn.
%
%   LICENSE:
%       This software is offered freely and without warranty under
%       the GNU (v3 or later) public license. See license file for
%       more information

addpath('E:\My algorithms\Matlab_SVM_tool\libsvm-3.20\matlab');

meanx=nanmean(training_feature);
stdx=nanstd(training_feature);

xtrain=bsxfun(@minus, training_feature, meanx);
xtrain=bsxfun(@rdivide, xtrain,stdx);

xtest=bsxfun(@minus, test_feature, meanx);
xtest=bsxfun(@rdivide, xtest,stdx);
model = svmtrain(training_class,xtrain, '-b 1');

[predict_all, accuracy_test, ytest_out] = svmpredict(test_class, xtest, model, '-b 1');

TP=length(find(test_class==1 & predict_all==1));
FP=length(find(test_class==0 & predict_all==1));
FN=length(find(test_class==1 & predict_all==0));
TN=length(find(test_class==0 & predict_all==0));
Se=TP/(TP+FN);
Sp=TN/(TN+FP);
Acc=(TP+TN)/length(test_class);
PPV=TP/(TP+FP);
NPV=TN/(TN+FN);
J=Se+Sp-1;
Index=[TP FN FP TN Se Sp Acc PPV NPV J];
end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
 % evaluate simulated data - ICA and NMF
% V ~ WH
% V -> N_pix x N_t
% W -> N_pix x N_comp - xth pixel of the ith components
% H -> N_copm x N_t - contribution of the i-th component in the time t

function separNMFICA(sep0, offset0, path_data, path_res, prename, niter, savethis, sep_how)

if ~exist('initval', 'var')
    initval = 0;
end

if ~exist('sep_how', 'var')
    sep_how = 'in';
end
ncomp = 3; %number of components to be separated + background

for rr = 1: length(offset0)
    for ll=1 : length(sep0)
        namedir = [prename num2str(100*sep0(ll)) 'offset_' num2str(offset0(rr))];
        cd ([path_data namedir])
        for mm=1:niter
            %reads the first...
            namefile = [namedir '-iter_' num2str(mm)];
            load ([namefile '.mat'])
            %cat the folowing,,,
            p.catitervec=[1];
            [dpixc, dveccr, blinkmat, p] = catsimul(namedir, p.catitervec);
            if sum(sep_how == 'i')>0 %ICA
                [icasig{mm}, A{mm}, W{mm}] = fastica (dveccr, 'numOfIC', ncomp, 'g', 'tanh');
                icapixICA{mm} = reshape(A{mm},p.nx, p.ny, ncomp);
            end

            if sum(sep_how == 'n')>0 %NMF
                % background estimation:
                %[out, bg(mm), bg_im]=backgroundoffset(dpixc);
                [out, bg(mm), bg_im]=backgroundoffset(dpixc, 'no', 5, 20, 8); %empirical values...

                dvec_bg = ones(p.nx*p.ny, 1);
                %                     dvec_bg = p.offset*ones(1, p.nx*p.ny); %changed for
                %                     offset 10...

                dvec_ind = squeeze(reshape(double(array2im(dpixc_ind)), p.nx*p.ny, 1, 2)); % vectors of resized images
                %                     sum_dvec_ind = sum(dvec_ind, 1);
                %                     dvec_ind = dvec_ind./repmat(sum_dvec_ind, p.nx*p.ny,1); %normlaized

                %                     winit = [f*dvec_ind'; dvec_bg];       %original 'true' points + background
                % % %                     winittmp = [dvec_ind, dvec_bg];
                winittmp = [rand(size(dvec_ind)), dvec_bg];
                sumw = sum(winittmp,1);
                winit = winittmp./repmat(sumw, p.nx*p.ny, 1); %normalized to 1
                f = mean(dveccr(:)-bg(mm))/mean(mean(winit(:, 1:2))); %ration of the data/psf
                %                     winit = [rand(ncomp, p.nx*p.ny); dvec_bg];
                blinkmatrand = rand(ncomp-1, p.Nt); %uniform random;
                hinit = [f*blinkmatrand; bg(mm)*sumw(ncomp)*ones(1, p.Nt)];             %random weights will be assigned to firts two and bg fixed
                % % %                     hinit = [f*blinkmat./repmat(mean(blinkmat,2),1,size(blinkmat,2)); bg(mm)*sumw(ncomp)*ones(1, p.Nt)];             %random weights will be assigned to firts two and bg fixed

                %                    [w{mm},h{mm}, wtrace{mm},wtrace{mm}]=nmf_test(double(dveccr'),ncomp+1,1,hinit,winit, [3], [3]);
                %                     [w{mm},h{mm}, wtrace,htrace,ddiv{mm}]=nmf_testconvD(double(dveccr'),ncomp+1,1,hinit,winit, [3], [3]);

                dvec_bg = ones(p.nx*p.ny, 1);
                winittmp = [dvec_ind, dvec_bg];
                sumw = sum(winittmp,1);
                winit = winittmp./repmat(sumw, p.nx*p.ny, 1); %normalized to 1
                f = mean(dveccr(:)-bg(mm))/mean(mean(winit(:, 1:2))); %ration of the data/psf

                p.Nt=1;
                p.meanblinkmat=mean(mean(blinkmat,2));
                p.htrue=[0.3; 0.7];
                htruef=p.meanblinkmat*p.htrue;

                %hinit = [p.meanblinkmat*rand(2,1); p.offset*sumw(ncomp)*ones(1, p.Nt)];

                    hinit = [p.meanblinkmat*rand(2,1); p.offset*sumw(ncomp)*ones(1, p.Nt)];
%                     hinit = [htruef; p.offset*sumw(ncomp)*ones(1, p.Nt)];
                    dveccr=dvec_ind*htruef+p.offset;
                    dpixcd = dip_image(reshape(dveccr,p.nx,p.ny));
                    dpixcdn = noise(dpixcd,'poisson');
                    dveccr=reshape(double(dpixcdn),p.nx*p.ny,1);
                    %[c,w{mm},h{mm}, X1,X2, dhr, minXr, miXvalr, mhdr, htrace, p]
                    [c,w{mm},h{mm}, X1,X2, dh, minX, miXval, mhd, htr, p]=nmf_S43(double(dveccr),ncomp,1,winit,hinit, [3], [3],p);
                    res=struct('c',c,'w',w,'h',h,'X1',X1,'X2',X2,'dh',dh,'minX',minX,'miXval',miXval, 'mhd',mhd, 'htr',htr, 'p',p);

                %hinit = [f*blinkmat./repmat(mean(blinkmat,2),1,size(blinkmat,2)); bg(mm)*sumw(ncomp)*ones(1, p.Nt)];             %random weights will be assigned to firts two and bg fixed
                %[c,w{mm},h{mm}, X1,X2, dht, minXt, miXvalt, mhdt, htrace,p]=nmf_S41(double(dveccr),ncomp,1,winit,hinit, [3], [3],p);
                qqq=[];
                icapixNMF{mm} = reshape(w{mm},p.nx,p.ny,ncomp);
            end

        end

        p.path_data = path_data;
        p.path_res = path_res;

        if savethis == 1
            fprintf('saving data \n');
            if ~(strcmp(path_res, path_data)) %not identical
                mkdir ([path_res namedir]);
                cd ([path_res namedir]);
            end
            %             save (p.namedir)
            %save ([namedir '_separMulti'],save ([namedir '_separMulti'],'p', 'X1','X2','dhr','minXr', 'miXvalr','dht'))
            save ([namedir '_DdviMap'],'res')
            writedata([],[],p,[namedir '_param'])
        end
    end
end

fprintf('\n')

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
% A SYMBOLIC DERIVATION OF THE JACOBIAN MANIPULATOR OF A KUKA LBR IIWA 14
%
% Copyright (C) 2016, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.
function jacobian_symbolic_kuka_iiwa
% link lengths

syms q1 q2 q3 q4 q5 q6 q7
robot = load_robot('KUKA', 'LBR_IIWA_R820_COP')

d = eval(robot.DH.d);
a = eval(robot.DH.a);
alpha = eval(robot.DH.alpha);

% matrices DH
A01 = dh_sym(q1, d(1), a(1), alpha(1));
A12 = dh_sym(q2, d(2), a(2), alpha(2));
A23 = dh_sym(q3, d(3), a(3), alpha(3));
A34 = dh_sym(q4, d(4), a(4), alpha(4));
A45 = dh_sym(q5, d(5), a(5), alpha(5));
A56 = dh_sym(q6, d(6), a(6), alpha(6));
A67 = dh_sym(q7, d(7), a(7), alpha(7));

A02 = A01*A12;
A03 = A02*A23;
A04 = A03*A34;
A05 = A04*A45;
A06 = A05*A56;
A07 = A06*A67;
A07 = simplify(A07)

z0 = [0 0 1]';
z1 = A01(1:3,3);
z2 = A02(1:3,3);
z3 = A03(1:3,3);
z4 = A04(1:3,3);
z5 = A05(1:3,3);
z6 = A06(1:3,3);

% simplify expressions
z2 = simplify(z2);
z3 = simplify(z3);
z4 = simplify(z4);
z5 = simplify(z5);
z6 = simplify(z6);

p07=A07(1:3,4);
p17=A07(1:3,4)-A01(1:3,4);
p27=A07(1:3,4)-A02(1:3,4);
p37=A07(1:3,4)-A03(1:3,4);
p47=A07(1:3,4)-A04(1:3,4);
p57=A07(1:3,4)-A05(1:3,4);
p67=A07(1:3,4)-A06(1:3,4);

% Jacobian in linear speed
Jv = [cross(z0, p07) cross(z1, p17) cross(z2, p27) cross(z3, p37) cross(z4, p47) cross(z5, p57) cross(z6, p67)];
% Jacobian in angular speed
Jw = [z0 z1 z2 z3 z4 z5 z6];

Jw = simplify(Jw)
Jv = simplify(Jv)
% singularities = det(Jv)
Js = [Jv; Jw];

q1 = 0.1
q2 = 0.1
q3 = 0.1
q4 = 0.1
q5 = 0.1
q6 = 0.1
q7 = 0.1

Js = eval(Js)

q = [q1 q2 q3 q4 q5 q6 q7]
J = manipulator_jacobian(robot, q)

J-Js

Ts = eval(A07)
T = directkinematic(robot, q)

T-Ts

function A = dh_sym(theta, d, a, alpha)
syms q1 q2 q3 q4 q5 q6 q7
% avoid almost zero elements in cos(alpha) and sin(alpha)
ca = cos(alpha);
sa = sin(alpha);
if abs(ca) < 1e-6
    ca = 0;
end
if abs(sa) < 1e-6
    sa = 0;
end

A=[cos(theta)  -ca*sin(theta)   sa*sin(theta)   a*cos(theta);
   sin(theta)   ca*cos(theta)  -sa*cos(theta)   a*sin(theta);
            0         sa             ca             d;
            0         0               0             1];

repo: ARTE


text:
function submit(partId, webSubmit)
%SUBMIT Submit your code and output to the pgm-class servers
%   SUBMIT() will connect to the pgm-class server and submit your solution
%   There is no penalty for submitting, so go ahead and try this!
%
%   If this function does not work for you, use the web-submission mechanism.
%   Call the submitWeb function, which will produce a file for the part you
%   wish to submit. Then, submit the file to the class servers using the
%   "Web Submission" button on the Programming Assignments page on the course
%   website.
%
%   webSubmit is a boolean variable that specifies whether to prepare
%   a file for web-submission (if webSubmit = 1) or to submit directly
%   to the server (if webSubmit = 0, or is not specified); you should call
%   submitWeb if you want to do a web-submission.
%

  fprintf('==\n== [pgm-class] Submitting Solutions | Programming Assignment %s\n==\n', ...
	  homework_id());

  if ~exist('partId', 'var') || isempty(partId)
    partId = promptPart();
  end

  if ~exist('webSubmit', 'var') || isempty(webSubmit)
    webSubmit = 0; % submit directly by default
  end

  % Check valid partId
  partNames = validParts();
  if ~isValidPartId(partId)
    fprintf('!! Invalid assignment part selected.\n');
    fprintf('!! Expected an integer from 1 to %d.\n', numel(partNames) + 1);
    fprintf('!! Submission Cancelled\n');
    return
  end

  if ~exist('228_login_data.mat','file')
    [login password] = loginPrompt();
    save('228_login_data.mat','login','password');
  else
    load('228_login_data.mat');
    [login password] = quickLogin(login,password);
    save('228_login_data.mat','login','password');
  end

  if isempty(login)
    fprintf('!! Submission Cancelled\n');
    return
  end

  fprintf('\n== Connecting to pgm-class ... ');
  if exist('OCTAVE_VERSION')
    fflush(stdout);
  end

  % Setup submit list
  if partId == numel(partNames) + 1
    submitParts = 1:numel(partNames);
  else
    submitParts = [partId];
  end

  for s = 1:numel(submitParts)
    % Submit this part
    partId = submitParts(s);

    for thisPartId = subParts(partId)
      if (~webSubmit) % submit directly to server

        [login, ch, signature, auxstring] = getChallenge(login, thisPartId);
        if isempty(login) || isempty(ch) || isempty(signature)
          % Some error occured, error string in first return element.
          fprintf('\n!! Error: %s\n\n', login);
          return
        end

        % Attempt Submission with Challenge
        ch_resp = challengeResponse(login, password, ch);

        [result, str] = submitSolution(login, ch_resp, thisPartId, ...
               output(thisPartId, auxstring), source(partId), signature);

        if (~isTest(thisPartId))
          partName = partNames{partId};
        else
          partName = [partNames{partId} ' (test)'];
        end

        fprintf('\n== [pgm-class] Submitted Assignment %s - Part %d - %s\n', ...
          homework_id(), partId, partName);
        fprintf('== %s\n', strtrim(str));

        if exist('OCTAVE_VERSION')
          fflush(stdout);
        end

      else % make web submission files

        [result] = submitSolutionWeb(login, thisPartId, output(thisPartId), ...
                            source(partId));
        result = base64encode(result);

        if (~isTest(thisPartId))
          partType = 'sample';
        else
          partType = 'test';
        end

        fprintf('\nSave as submission file [submit_pa%s_part%d_%s.txt (enter to accept default)]:', ...
          homework_id(), partId, partType);
        saveAsFile = input('', 's');
        if (isempty(saveAsFile))
          saveAsFile = sprintf('submit_pa%s_part%d_%s.txt', homework_id(), partId, partType);
        end

        fid = fopen(saveAsFile, 'w');
        if (fid)
          fwrite(fid, result);
          fclose(fid);
          fprintf('\nSaved your solutions to %s.\n\n', saveAsFile);
          fprintf(['You can now submit your solutions using the ' ...
             'Web Submission button\non the Assignments page\n']);
        else
          fprintf('Unable to save to %s\n\n', saveAsFile);
          fprintf(['You can create a submission file by saving the \n' ...
             'following text in a file: (press enter to continue)\n\n']);
          pause;
          fprintf(result);
        end

      end
    end
  end
end

% ====================== SERVER CONFIGURATION ===========================

% ***************** REMOVE -staging WHEN YOU DEPLOY *********************
function url = site_url()
  url = 'http://www.coursera.org/pgm';
end

function url = challenge_url()

  url = [site_url() '/assignment/challenge'];
end

function url = submit_url()
  url = [site_url() '/assignment/submit'];
end

% ================== CONFIGURABLES FOR EACH HOMEWORK ==================

function id = homework_id()
  id = '1';
end

function [partNames] = validParts()
  partNames = { 'Constructing a Credit Network', ...
                'Factor Product', ...
                'Factor Marginalization', ...
                'Observing Evidence', ...
                'Computing the Joint Distribution', ...
                'Computing Marginals'
              };
end

function srcs = sources()
  % Separated by part
  srcs = { { 'Credit_net.net'}, ...
           { 'FactorProduct.m'}, ...
           { 'FactorMarginalization.m'}, ...
           { 'ObserveEvidence.m'}, ...
           { 'ComputeJointDistribution.m'}, ...
           { 'ComputeMarginal.m'}
         };
end

% defines the shown part to back-end part mappings
function parts = subParts(partId)
  if (partId == 1)
      parts = [1];
  else
      parts = [2 * (partId - 1), 2 * (partId - 1) + 1];
  end
end

% specifies which parts are test parts
function result = isTest(partId)
  if (bitand(partId, 1))
    result = true;
  else
    result = false;
  end
end

function out = output(partId, auxstring)
  load submit_input;
  if partId == 1
    [F, names, assignments] = ConvertNetwork('Credit_net.net');
    [F, names, assignments] = RelabelVars(F, names, assignments);
    if (ValidateNetwork(F, names, assignments))
      out = SerializeFactorsFg(F);
    else
      out = '';
    end
  elseif partId == 2
    fp = FactorProduct(PART2.SAMPLEINPUT{:});
    out = SerializeFactorsFg(fp);
  elseif partId == 3
    fp = [FactorProduct(PART2.INPUT1{:}), FactorProduct(PART2.INPUT2{:})];
    out = SerializeFactorsFg(fp);
  elseif partId == 4
    fm = [FactorMarginalization(PART3.SAMPLEINPUT{:})];
    out = SerializeFactorsFg(fm);
  elseif partId == 5
    fm = [FactorMarginalization(PART3.INPUT1{:}), FactorMarginalization(PART3.INPUT2{:})];
    out = SerializeFactorsFg(fm);
  elseif partId == 6
    oe = [ObserveEvidence(PART4.SAMPLEINPUT{:})];
    out = SerializeFactorsFg(oe);
  elseif partId == 7
    oe = [ObserveEvidence(PART4.INPUT1{:})];
    out = SerializeFactorsFg(oe);
  elseif partId == 8
    jd = [ComputeJointDistribution(PART5.SAMPLEINPUT{:})];
    out = SerializeFactorsFg(jd);
  elseif partId == 9
    jd = [ComputeJointDistribution(PART5.INPUT1{:})];
    out = SerializeFactorsFg(jd);
  elseif partId == 10
    cm = [ComputeMarginal(PART6.SAMPLEINPUT{:})];
    out = SerializeFactorsFg(cm);
  elseif partId == 11
    cm = [ComputeMarginal(PART6.INPUT1{:}), ComputeMarginal(PART6.INPUT2{:}), ComputeMarginal(PART6.INPUT3{:}), ComputeMarginal(PART6.INPUT4{:})];
    out = SerializeFactorsFg(cm);
  end
end

function ok = ValidateNetwork(F, names, assignments)
  % Check number of variables in network
  assert(numel(names) == 8, sprintf('Error: your network should have 8 variables but it only has %d. Did you delete some variables?', numel(names)));

  % Check number of CPDs in network
  % If this happens, it's most likely a malformed .net file, or else it is
  % a bug in ConvertNetwork.m
  assert(numel(F) == 8, sprintf('Error: your network should have 8 CPDs but it only has %d.', numel(F)));

  % Check that the assignment of variables is the same as in the original
  valNames = cell(8, 1);
  valNames{1} = {'Low','High'};
  valNames{2} = {'High','Medium','Low'};
  valNames{3} = {'Positive','Negative'};
  valNames{4} = {'High','Medium','Low'};
  valNames{5} = {'Excellent','Acceptable','Unacceptable'};
  valNames{6} = {'Promising','Not_promising'};
  valNames{7} = {'Reliable','Unreliable'};
  valNames{8} = {'Between16and21','Between22and64','Over65'};

  for i = 1:numel(valNames)
    assert(all(strcmp(valNames{i}, assignments{i})), ...
      sprintf('Error: variable values/value ordering for variable %s does not match the original. Did you change the order/name of some variable values?', names{i}));
  end

  ok = true;
end

% Relabel variable ids according to the grader ordering
function [newF, newNames, newAssignments] = RelabelVars(F, names, assignments)
  varNames = cell(8, 1);
  varNames{1} = 'DebtIncomeRatio';
  varNames{2} = 'Assets';
  varNames{3} = 'CreditWorthiness';
  varNames{4} = 'Income';
  varNames{5} = 'PaymentHistory';
  varNames{6} = 'FutureIncome';
  varNames{7} = 'Reliability';
  varNames{8} = 'Age';

  newToOrig = zeros(8, 1);
  origToNew = zeros(8, 1);
  for i = 1:8
    result = find(strcmp(varNames{i}, names), 1);
    if (isempty(result)) % no match!
      error('Variable %s not found in network! Did you change the variable names/identifiers?', varNames{i});
    else
      newToOrig(result) = i;
      origToNew(i) = result;
    end
  end

  newNames = names(origToNew(:));
  newAssignments = assignments(origToNew(:));

  % relabel the vars
  newF = F;
  for i = 1:numel(F)
    newF(i).var = newToOrig(F(i).var);
  end
end

function out = SerializeFactorsFg(F)
% Serializes a factor struct array into the .fg format for libDAI
% http://cs.ru.nl/~jorism/libDAI/doc/fileformats.html
%
% To avoid incompatibilities with EOL markers, make sure you write the
% string to a file using the appropriate file type ('wt' for windows, 'w'
% for unix)

lines = cell(5*numel(F) + 1, 1);

lines{1} = sprintf('%d\n', numel(F));
lineIdx = 2;
for i = 1:numel(F)
  lines{lineIdx} = sprintf('\n%d\n', numel(F(i).var));
  lineIdx = lineIdx + 1;

  lines{lineIdx} = sprintf('%s\n', num2str(F(i).var(:)')); % ensure that we put in a row vector
  lineIdx = lineIdx + 1;

  lines{lineIdx} = sprintf('%s\n', num2str(F(i).card(:)')); % ensure that we put in a row vector
  lineIdx = lineIdx + 1;

  lines{lineIdx} = sprintf('%d\n', numel(F(i).val));
  lineIdx = lineIdx + 1;

  % Internal storage of factor vals is already in the same indexing order
  % as what libDAI expects, so we don't need to convert the indices.
  vals = [0:(numel(F(i).val) - 1); F(i).val(:)'];
  lines{lineIdx} = sprintf('%d %0.8g\n', vals);
  lineIdx = lineIdx + 1;
end

out = sprintf('%s', lines{:});

end

% ========================= CHALLENGE HELPERS =========================

function src = source(partId)
  src = '';
  src_files = sources();
  if partId <= numel(src_files)
      flist = src_files{partId};
      for i = 1:numel(flist)
          fid = fopen(flist{i});
          if (fid == -1)
            error('Error opening %s (is it missing?)', flist{i});
          end
          line = fgets(fid);
          while ischar(line)
            src = [src line];
            line = fgets(fid);
          end
          fclose(fid);
          src = [src '||||||||'];
      end
  end
end

function ret = isValidPartId(partId)
  partNames = validParts();
  ret = (~isempty(partId)) && (partId >= 1) && (partId <= numel(partNames) + 1);
end

function partId = promptPart()
  fprintf('== Select which part(s) to submit:\n', ...
          homework_id());
  partNames = validParts();
  srcFiles = sources();
  for i = 1:numel(partNames)
    fprintf('==   %d) %s [', i, partNames{i});
    fprintf(' %s ', srcFiles{i}{:});
    fprintf(']\n');
  end
  fprintf('==   %d) All of the above \n==\nEnter your choice [1-%d]: ', ...
          numel(partNames) + 1, numel(partNames) + 1);
  selPart = input('', 's');
  partId = str2num(selPart);
  if ~isValidPartId(partId)
    partId = -1;
  end
end

function [email,ch,signature,auxstring] = getChallenge(email, part)
  str = urlread(challenge_url(), 'post', {'email_address', email, 'assignment_part_sid', [homework_id() '-' num2str(part)], 'response_encoding', 'delim'});

  str = strtrim(str);
  r = struct;
  while(numel(str) > 0)
    [f, str] = strtok (str, '|');
    [v, str] = strtok (str, '|');
    r = setfield(r, f, v);
  end

  email = getfield(r, 'email_address');
  ch = getfield(r, 'challenge_key');
  signature = getfield(r, 'state');
  auxstring = getfield(r, 'challenge_aux_data');
end

function [result, str] = submitSolutionWeb(email, part, output, source)

  result = ['{"assignment_part_sid":"' base64encode([homework_id() '-' num2str(part)], '') '",' ...
            '"email_address":"' base64encode(email, '') '",' ...
            '"submission":"' base64encode(output, '') '",' ...
            '"submission_aux":"' base64encode(source, '') '"' ...
            '}'];
  str = 'Web-submission';
end

function [result, str] = submitSolution(email, ch_resp, part, output, ...
                                        source, signature)

  params = {'assignment_part_sid', [homework_id() '-' num2str(part)], ...
            'email_address', email, ...
            'submission', base64encode(output, ''), ...
            'submission_aux', base64encode(source, ''), ...
            'challenge_response', ch_resp, ...
            'state', signature};

  str = urlread(submit_url(), 'post', params);

  % Parse str to read for success / failure
  result = 0;

end

% =========================== LOGIN HELPERS ===========================

function [login password] = loginPrompt()
  % Prompt for password
  [login password] = basicPrompt();

  if isempty(login) || isempty(password)
    login = []; password = [];
  end
end

function [login password] = basicPrompt()
  login = input('Login (Email address): ', 's');
  password = input('Submission Password (from Assignments page): ', 's');
end

function [login password] = quickLogin(login,password)
  cont_token = input(['You are currently logged in as ' login '.\nIs this you? (y/n - type n to reenter password)'],'s');
  if(isempty(cont_token) || cont_token(1)=='Y'||cont_token(1)=='y')
    return;
  else
    [login password] = loginPrompt();
  end
end

function [str] = challengeResponse(email, passwd, challenge)
  str = sha1([challenge passwd]);
end

% =============================== SHA-1 ================================

function hash = sha1(str)

  % Initialize variables
  h0 = uint32(1732584193);
  h1 = uint32(4023233417);
  h2 = uint32(2562383102);
  h3 = uint32(271733878);
  h4 = uint32(3285377520);

  % Convert to word array
  strlen = numel(str);

  % Break string into chars and append the bit 1 to the message
  mC = [double(str) 128];
  mC = [mC zeros(1, 4-mod(numel(mC), 4), 'uint8')];

  numB = strlen * 8;
  if exist('idivide')
    numC = idivide(uint32(numB + 65), 512, 'ceil');
  else
    numC = ceil(double(numB + 65)/512);
  end
  numW = numC * 16;
  mW = zeros(numW, 1, 'uint32');

  idx = 1;
  for i = 1:4:strlen + 1
    mW(idx) = bitor(bitor(bitor( ...
                  bitshift(uint32(mC(i)), 24), ...
                  bitshift(uint32(mC(i+1)), 16)), ...
                  bitshift(uint32(mC(i+2)), 8)), ...
                  uint32(mC(i+3)));
    idx = idx + 1;
  end

  % Append length of message
  mW(numW - 1) = uint32(bitshift(uint64(numB), -32));
  mW(numW) = uint32(bitshift(bitshift(uint64(numB), 32), -32));

  % Process the message in successive 512-bit chs
  for cId = 1 : double(numC)
    cSt = (cId - 1) * 16 + 1;
    cEnd = cId * 16;
    ch = mW(cSt : cEnd);

    % Extend the sixteen 32-bit words into eighty 32-bit words
    for j = 17 : 80
      ch(j) = ch(j - 3);
      ch(j) = bitxor(ch(j), ch(j - 8));
      ch(j) = bitxor(ch(j), ch(j - 14));
      ch(j) = bitxor(ch(j), ch(j - 16));
      ch(j) = bitrotate(ch(j), 1);
    end

    % Initialize hash value for this ch
    a = h0;
    b = h1;
    c = h2;
    d = h3;
    e = h4;

    % Main loop
    for i = 1 : 80
      if(i >= 1 && i <= 20)
        f = bitor(bitand(b, c), bitand(bitcmp(b), d));
        k = uint32(1518500249);
      elseif(i >= 21 && i <= 40)
        f = bitxor(bitxor(b, c), d);
        k = uint32(1859775393);
      elseif(i >= 41 && i <= 60)
        f = bitor(bitor(bitand(b, c), bitand(b, d)), bitand(c, d));
        k = uint32(2400959708);
      elseif(i >= 61 && i <= 80)
        f = bitxor(bitxor(b, c), d);
        k = uint32(3395469782);
      end

      t = bitrotate(a, 5);
      t = bitadd(t, f);
      t = bitadd(t, e);
      t = bitadd(t, k);
      t = bitadd(t, ch(i));
      e = d;
      d = c;
      c = bitrotate(b, 30);
      b = a;
      a = t;

    end
    h0 = bitadd(h0, a);
    h1 = bitadd(h1, b);
    h2 = bitadd(h2, c);
    h3 = bitadd(h3, d);
    h4 = bitadd(h4, e);

  end

  hash = reshape(dec2hex(double([h0 h1 h2 h3 h4]), 8)', [1 40]);

  hash = lower(hash);

end

function ret = bitadd(iA, iB)
  ret = double(iA) + double(iB);
  ret = bitset(ret, 33, 0);
  ret = uint32(ret);
end

function ret = bitrotate(iA, places)
  t = bitshift(iA, places - 32);
  ret = bitshift(iA, places);
  ret = bitor(ret, t);
end

% =========================== Base64 Encoder ============================
% Thanks to Peter John Acklam
%

function y = base64encode(x, eol)
%BASE64ENCODE Perform base64 encoding on a string.
%
%   BASE64ENCODE(STR, EOL) encode the given string STR.  EOL is the line ending
%   sequence to use; it is optional and defaults to '\n' (ASCII decimal 10).
%   The returned encoded string is broken into lines of no more than 76
%   characters each, and each line will end with EOL unless it is empty.  Let
%   EOL be empty if you do not want the encoded string broken into lines.
%
%   STR and EOL don't have to be strings (i.e., char arrays).  The only
%   requirement is that they are vectors containing values in the range 0-255.
%
%   This function may be used to encode strings into the Base64 encoding
%   specified in RFC 2045 - MIME (Multipurpose Internet Mail Extensions).  The
%   Base64 encoding is designed to represent arbitrary sequences of octets in a
%   form that need not be humanly readable.  A 65-character subset
%   ([A-Za-z0-9+/=]) of US-ASCII is used, enabling 6 bits to be represented per
%   printable character.
%
%   Examples
%   --------
%
%   If you want to encode a large file, you should encode it in chunks that are
%   a multiple of 57 bytes.  This ensures that the base64 lines line up and
%   that you do not end up with padding in the middle.  57 bytes of data fills
%   one complete base64 line (76 == 57*4/3):
%
%   If ifid and ofid are two file identifiers opened for reading and writing,
%   respectively, then you can base64 encode the data with
%
%      while ~feof(ifid)
%         fwrite(ofid, base64encode(fread(ifid, 60*57)));
%      end
%
%   or, if you have enough memory,
%
%      fwrite(ofid, base64encode(fread(ifid)));
%
%   See also BASE64DECODE.

%   Author:      Peter John Acklam
%   Time-stamp:  2004-02-03 21:36:56 +0100
%   E-mail:      pjacklam@online.no
%   URL:         http://home.online.no/~pjacklam

   if isnumeric(x)
      x = num2str(x);
   end

   % make sure we have the EOL value
   if nargin < 2
      eol = sprintf('\n');
   else
      if sum(size(eol) > 1) > 1
         error('EOL must be a vector.');
      end
      if any(eol(:) > 255)
         error('EOL can not contain values larger than 255.');
      end
   end

   if sum(size(x) > 1) > 1
      error('STR must be a vector.');
   end

   x   = uint8(x);
   eol = uint8(eol);

   ndbytes = length(x);                 % number of decoded bytes
   nchunks = ceil(ndbytes / 3);         % number of chunks/groups
   nebytes = 4 * nchunks;               % number of encoded bytes

   % add padding if necessary, to make the length of x a multiple of 3
   if rem(ndbytes, 3)
      x(end+1 : 3*nchunks) = 0;
   end

   x = reshape(x, [3, nchunks]);        % reshape the data
   y = repmat(uint8(0), 4, nchunks);    % for the encoded data

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Split up every 3 bytes into 4 pieces
   %
   %    aaaaaabb bbbbcccc ccdddddd
   %
   % to form
   %
   %    00aaaaaa 00bbbbbb 00cccccc 00dddddd
   %
   y(1,:) = bitshift(x(1,:), -2);                  % 6 highest bits of x(1,:)

   y(2,:) = bitshift(bitand(x(1,:), 3), 4);        % 2 lowest bits of x(1,:)
   y(2,:) = bitor(y(2,:), bitshift(x(2,:), -4));   % 4 highest bits of x(2,:)

   y(3,:) = bitshift(bitand(x(2,:), 15), 2);       % 4 lowest bits of x(2,:)
   y(3,:) = bitor(y(3,:), bitshift(x(3,:), -6));   % 2 highest bits of x(3,:)

   y(4,:) = bitand(x(3,:), 63);                    % 6 lowest bits of x(3,:)

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Now perform the following mapping
   %
   %   0  - 25  ->  A-Z
   %   26 - 51  ->  a-z
   %   52 - 61  ->  0-9
   %   62       ->  +
   %   63       ->  /
   %
   % We could use a mapping vector like
   %
   %   ['A':'Z', 'a':'z', '0':'9', '+/']
   %
   % but that would require an index vector of class double.
   %
   z = repmat(uint8(0), size(y));
   i =           y <= 25;  z(i) = 'A'      + double(y(i));
   i = 26 <= y & y <= 51;  z(i) = 'a' - 26 + double(y(i));
   i = 52 <= y & y <= 61;  z(i) = '0' - 52 + double(y(i));
   i =           y == 62;  z(i) = '+';
   i =           y == 63;  z(i) = '/';
   y = z;

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Add padding if necessary.
   %
   npbytes = 3 * nchunks - ndbytes;     % number of padding bytes
   if npbytes
      y(end-npbytes+1 : end) = '=';     % '=' is used for padding
   end

   if isempty(eol)

      % reshape to a row vector
      y = reshape(y, [1, nebytes]);

   else

      nlines = ceil(nebytes / 76);      % number of lines
      neolbytes = length(eol);          % number of bytes in eol string

      % pad data so it becomes a multiple of 76 elements
      y = [y(:) ; zeros(76 * nlines - numel(y), 1)];
      y(nebytes + 1 : 76 * nlines) = 0;
      y = reshape(y, 76, nlines);

      % insert eol strings
      eol = eol(:);
      y(end + 1 : end + neolbytes, :) = eol(:, ones(1, nlines));

      % remove padding, but keep the last eol string
      m = nebytes + neolbytes * (nlines - 1);
      n = (76+neolbytes)*nlines - neolbytes;
      y(m+1 : n) = '';

      % extract and reshape to row vector
      y = reshape(y, 1, m+neolbytes);

   end

   % output is a character array
   y = char(y);

end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function bool = istequal(t1,t2)

%ISTEQUAL: Compare 2 Matlab time values for equality
%
%USAGE: bool = istequal(t1,t2)
%
%INPUTS: t1,t2 - Matlab formatted date numbers
%
%OUTPUTS: bool - 1 or 0
%
% See also ADD_SST, CHK_T, COMPARE_SST, DELETE_SST, EXTRACT_SST, IS_SST,
%          ISTEQUAL, MERGE_SST, SEARCH_SST, SORT_SST, NAN2SST, SSD2SST,
%          SST2NAN, SST2SSD, SST2VAL, SST2WFA, WFA2SST
%
% Author: Dane Ketner, Alaska Volcano Observatory
% $Date$
% $Revision$

t1 = round(t1*24*60*60*1000); % Units: seconds/1000
t2 = round(t2*24*60*60*1000); % Units: seconds/1000
bool = isequal(t1,t2);

author: geoscience-community-codes

repo: GISMO


text:
% VBRFA2011_TESTBED_PLOT - Make the Testbed plots for the VBRFA journal
% paper.

% Last modified 2011-06-15
% Copyright (c) Jaakko Luttinen (jaakko.luttinen@aalto.fi)

function data = vbrfa2011_testbed_plot(model)

% Result files
folder = '/share/climate/jluttine/testbed';
%files = {'testbed_results_ind-t_D=30_20110616'};
files = {'testbed_results_gaussian_D=30_20110617', ...
         'testbed_results_multi-t_D=30_20110617', ...
         'testbed_results_ind-t_D=30_20110617', ...
         'testbed_results_laplace_D=30_20110617'};

% Pick stations from: 12 13 20 21 28 48 57 59
corrupted = [12 13 20 21 28 48 57 59];
stations = [5 corrupted([8 7 6 1 3])];

%
% PLOT DATA OBSERVATIONS
%

data = testbed_loaddata();
data = testbed_preprocess(data);
return
plot_timeseries(data.time, ...
                data.observations(stations,:));
%return

print('-depsc2', '-loose', ...
      '/home/jluttine/papers/rpca_journal/fig_testbed_timeseries_data');
%
% PLOT RESULTS FOR EACH METHOD
%

suffix = {'gaussian', 'multi-t', 'ind', 'laplace'};

% Plot
if nargin < 1
  model = 1:length(files);
end
for n=model
  fprintf('Loading results from file:\n%s\n', files{n});
  Q = load(files{n}, 'W', 'X', 'W_struct');
% $$$   plot_loadings(Q,suffix{n});
  plot_reconstructions(Q,stations,suffix{n});
  %figure, semilogy(Q.W_struct.alpha)
  %fprintf('alpha(%d): %.4e\n', [1:length(Q.w); Q.w(:)']);
end

% $$$ %
% $$$ % PLOT THE STATIONS
% $$$ %
% $$$
% $$$ fig = figure();
% $$$
% $$$ % Plot the research area on the map of Finland
% $$$ hax(1) = subplot(1,2,1);
% $$$ testbed_plot_finland();
% $$$
% $$$ % Plot the stations on a topographic map
% $$$ hax(2) = subplot(1,2,2);
% $$$ set(hax(2), 'FontSize', 9);
% $$$ testbed_plot_topography();
% $$$ hold on;
% $$$ testbed_coast();
% $$$ map_grid();
% $$$ testbed_plot_stations(data.coordinates);
% $$$ hcb = colorbar('peer',hax(2),'FontSize',9);
% $$$
% $$$ % Set layout
% $$$ bottom = 0.07;
% $$$ top = 0.92;
% $$$ set(hax(1), 'Units','Normalized', 'Position',[0.0 bottom 0.22 (top-bottom)]);
% $$$ set(hax(2), 'Units','Normalized', 'Position',[0.27 bottom 0.63 (top-bottom)]);
% $$$ set(hcb, 'Units','Normalized', 'Position',[0.91 bottom 0.02 (top-bottom)]);
% $$$ set_figure_size(15,6.3,fig);
% $$$
% $$$ % $$$ % Print
% $$$ % $$$ print('-depsc2', '-loose', ...
% $$$ % $$$       '/home/jluttine/papers/vbrfa/figures_journal/fig_testbed_stations');
% $$$
% $$$

function plot_timeseries(time,Y)

% Plot timeseries
hax = tsplot(time, Y, 'k');

% Set the range of the x-axis
N = length(hax);
set(hax, 'xlim', [min(time), max(time)]);

% Put date ticks on x-axis
dateticks(hax, 'x', 'mmmyy', 'keeplimits');
set(hax(1:(N-1)), 'xticklabel', []);
set_ticks_fontsize(hax,9);
xtick = get(hax(end), 'xtick');
set(hax, 'xtick', xtick(2:end));

% Remove ticks from y-axis
set(hax, 'ylim', [-50 50], 'yticklabel', [], 'ytick', []);

% Set the layout of the subplots
set_subplot_positions(hax, N, 1, [0.01 0.01 0.01 0.07], [0.02 0.02]);
set_figure_size(15,9)

function plot_reconstructions(Q,stations,suffix)

% Reconstructions
Yh = Q.W(:,stations)' * Q.X;

% Remove time instances with no observations to make the plot more clear
data = load(['/share/climate/jluttine/testbed/' ...
             'testbed_vbrfa2011_traindata']);
% $$$ data = testbed_loaddata();
% $$$ data = testbed_preprocess(data);
ind = sum(~isnan(data.observations),1) == 0;
Yh(:,ind) = nan;

% Plotplotplot
plot_timeseries(data.time, Yh);

% Print figure
print('-depsc2', '-loose', ...
      ['/home/jluttine/papers/rpca_journal/fig_testbed_timeseries_', ...
       suffix]);

function plot_loadings(Q,suffix)

%
% Plot spatial components
%

data = testbed_loaddata();
data = testbed_preprocess(data);

components = 1:4;
hax = testbed_loadings(Q.W(components,:)', ...
                       data.coordinates, ...
                       'resolution', [200 150], ...
                       'method', 'rbf');
% Subplot layout
set_subplot_positions(hax, 2, 2, [0.01 0.01 0.08 0.06], [0.08 0.1]);
% Colorbar layout
colorbars(hax, ...
          'Size', 0.02, ...
          'Separation', 0.0, ...
          'Location', 'EastOutside', ...
          'FontSize', 9);
% X-labels (and layout)
for n=1:numel(components)
  xlabels(hax(n),sprintf('(%d)', components(n)), 'Separation',0.03, 'FontSize', 9);
end
% Print figure
set_figure_size(15,8.5)
% $$$ print('-depsc2', '-loose', ...
% $$$       ['/home/jluttine/papers/vbrfa/figures_journal/fig_testbed_loadings_', ...
% $$$        suffix]);

% Plot reconstructions of problematic stations
%testbed_rcplot(Q, 'stations', [12 13 20 21 28 48 57 59]);

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function imf = bwtfilter(im,f)
% function imf = bwtfilter(im,f)
%
% BWT decimated filtering
% Version 1.2
%
% This is faster than doing the complete filtering operation with conv2 and
% then decimating the result
%
% Arguments:
%  im: A square image, with side length a multiple of 3
%  f:  A 3x3 filter
%
% Result:
%  imf: im is filtered with f, and the result is decimated, taking every
%       3rd value in x and y. This code is equivalent to:
%       imf = conv2(im, rot90(rot90(f)),'same');
%       imf = imf(2:3:end,2:3:end);
%
% Citation:
%  Willmore B, Prenger RJ, Wu MC and Gallant JL (2008). The Berkeley
%  Wavelet Transform: A biologically-inspired orthogonal wavelet transform.
%  Neural Computation 20:6, 1537-1564
%
% The article is available at:
%  <http://dx.doi.org/10.1162/neco.2007.05-07-513>
%
% Copyright (c) 2008 Ben Willmore
%
% Permission is hereby granted, free of charge, to any person
% obtaining a copy of this software and associated documentation
% files (the "Software"), to deal in the Software without
% restriction, including without limitation the rights to use,
% copy, modify, merge, publish, distribute, sublicense, and/or sell
% copies of the Software, and to permit persons to whom the
% Software is furnished to do so, subject to the following
% conditions:
%
% The above copyright notice and this permission notice shall be
% included in all copies or substantial portions of the Software.
%
% 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.

sz = size(im);

if (length(sz) ~= 2) || (sz(1) ~= sz(2))
  disp('Input must be square');
  decomp = nan;
  return;
end

sz = sz(1);
ssz = sz/3;

if ( (ssz-floor(ssz)) > abs(ssz)*eps )
  fprintf('Side length must be a multiple of 3');
  decomp = nan;
  return;
end

f_rep = repmat(f,ssz,ssz);

imdot = im .* f_rep;

imf = zeros(ssz);

for yy = 1:3
  for xx = 1:3
    imf = imf + imdot(yy:3:end,xx:3:end);
  end
end


text:
function pers=ordinal_persistence(S)
% ORDINAL_PERSISTENCE computes the persistence of an ordered multilayer partition
%
% Version: 2.2.0
% Date: Thu 11 Jul 2019 12:25:42 CEST
%
% pers = ORDINAL_PERSISTENCE(S) with a single multilayer partition or a
% cell of multilayer partitions S (with S{i} the ith multilayer partition
% of cell S) computes the "persistence" of ordered multilayer partitions. The
% output pers is a vector such that pers(i) is the persistence of multilayer
% partition S{i} (stored as an N*T matrix where N is the number of nodes in
% each layer and T is the number of layers). For a multilayer partition with
% ordered layers, the value of persistence is the number of nodes that do
% not change community assignments between consecutive layers (see Bazzi
% et al. 2016 for more detail). Ordinal persistence varies between 0 and 1.
% Ordinal persistence is related to the multilayer quality function
% developed in Mucha et al. 2010 when one uses ordinal interlayer coupling.
%
%   References:
%
%     Mucha, Peter J., Thomas Richardson, Kevin Macon, Mason A. Porter, and
%     Jukka-Pekka Onnela. "Community Structure in Time-Dependent,
%     Multiscale, and Multiplex Networks," Science 328, 876-878 (2010).
%
%     Bazzi, Marya, Mason A. Porter, Stacy Williams, Mark McDonald, Daniel
%     J. Fenn, and Sam D. Howison. "Community Detection in Temporal Multilayer
%     Networks, with an Application to Correlation Networks", MMS: A SIAM
%     Interdisciplinary Journal 14, 1-41 (2016).

if ~iscell(S)
    S={S};
end

[N,T]=size(S{1});
pers=zeros(length(S),1);
for i=1:length(S)
    pers(i)=sum(sum(S{i}(:,1:end-1)==S{i}(:,2:end)))/(N*(T-1));
end

end

author: GenLouvain

repo: GenLouvain


text:
function check_dl_model( model )
% Download model file
if ~exist(model,'file')
  [~, baseModel] = fileparts(model);
  fprintf('Downloading base model file: %s ...\n', baseModel);
  mkdir_if_missing(fileparts(model)) ;
  urlwrite(...
  ['http://ftp.tugraz.at/pub/feichtenhofer/detect-track/models/' baseModel '.caffemodel'], ...
    model) ;
end

end

author: feichtenhofer

repo: Detect-Track


text:
function model = rmSearchFit_oneGaussian(model, data, params, wProcess, t)
% rmSearchFit_oneGaussian - wrapper for 'fine' one Gaussian fit
%
% model = rmSearchFit_oneGaussian(model, data, params, wProcess, t);
%
% 2008/01 SOD: split of from rmSearchFit.
% 2010/02 SOD: cleanup.

% fminsearch options
searchOptions = params.analysis.fmins.options;
expandRange   = params.analysis.fmins.expandRange;

% convert to double just in case
params.analysis.X = double(params.analysis.X);
params.analysis.Y = double(params.analysis.Y);
params.analysis.allstimimages = double(params.analysis.allstimimages);
data = double(data);

% get starting upper and lower range and reset TolFun
% (raw rss computation (similar to norm) and TolFun adjustments)
[range TolFun] = rmSearchFit_range(params,model,data);

% amount of negative fits
nNegFit  = 0;
vethresh = params.analysis.fmins.vethresh;
trends   = t.trends;
t_id     = t.dcid+1;

%-----------------------------------
% Go for each voxel
%-----------------------------------
progress = 0;tic;
for ii = 1:numel(wProcess),
    % progress monitor (10 dots)
    if floor(ii./numel(wProcess)*10)>progress,
        % print out estimated time left
        if progress==0,
            esttime = toc.*10;
            if floor(esttime./3600)>0,
                fprintf(1,'[%s]:Estimated processing time: %d voxels: %d hours.\n',...
                    mfilename,numel(wProcess),ceil(esttime./3600));
            else
                fprintf(1,'[%s]:Estimated processing time: %d voxels: %d minutes.\n',...
                    mfilename,numel(wProcess),ceil(esttime./60));
            end;
            fprintf(1,'[%s]:Nonlinear optimization (x,y,sigma):',mfilename);
        end;
        fprintf(1,'.');drawnow;
        progress = progress + 1;
    end;

    % volume index
    vi = wProcess(ii);
    vData = data(:,ii);

    % reset tolFun: Precision of evaluation function.
    % We define RMS improvement relative to the initial raw 'no-fit' data
    % RMS. So, 1 means stop if there is less than 1% improvement on the fit:
    % searchOptions = optimset(searchOptions,'tolFun',optimget(params.analysis.fmins.options,'tolFun')./100.*rawrss);
    % optimset and optimget are a little slow so:
    searchOptions.TolFun = TolFun(ii);

    % actual fitting routine
    if searchOptions.MaxIter>0
        outParams = ...
            fmincon(@(x) rmModelSearchFit_oneGaussian(x,vData,...
            params.analysis.X,...
            params.analysis.Y,...
            params.analysis.allstimimages,trends),...
            range.start(:,vi),[],[],[],[],range.lower(:,vi),range.upper(:,vi),...
            @(x) distanceCon(x,range.start(:,vi),range.step(:,vi).*expandRange),searchOptions);
    else
        outParams = range.start(:,vi);
    end

    %[ outParams range.lower(:,vi) range.start(:,vi) range.upper(:,vi)]

    % make RF, prediction and get rss,b
    Xv = params.analysis.X-outParams(1);
    Yv = params.analysis.Y-outParams(2);
    rf = exp( (Yv.*Yv + Xv.*Xv) ./ (-2.*(outParams(3).^2)) );

    % GLU: checked, they are basically the same
    % {
    % rf=rfGaussian2d(XY{1}, XY{2}, rf.sigmaMajor,rf.sigmaMinor,rf.Theta, ...
    %                            rf.Centerx0,rf.Centery0);
    rf2=rfGaussian2d(params.analysis.X, params.analysis.Y, ...
                    outParams(3),outParams(3),0, ...
                               outParams(1),outParams(2));
    %}

    X  = [params.analysis.allstimimages * rf trends];
    b    = pinv(X)*vData;
    rss  = norm(vData-X*b).^2;

    % store results only if the first beta is positive, somehow fmincon
    % outputs negative fits. If the fit is negative keep old (grid) fit. We
    % do adjust the rss, so it won't be accidentally counted as a 'good'
    % fit.
    if b(1)>0,
        model.x0(vi)   = outParams(1);
        model.y0(vi)   = outParams(2);
        model.s(vi)    = outParams(3);
        model.s_major(vi)    = outParams(3);
        model.s_minor(vi)    = outParams(3);
        model.s_theta(vi)    = 0;
        model.rss(vi)  = rss;
        model.b([1 t_id],vi)  = b;
    else
        % change the percent variance explained to be just under the
        % current vethresh. So it counts as a 'coarse'-fit but can still be
        % included in later 'fine'-fits
        model.rss(vi)  = (1-max((vethresh-0.01),0)).*model.rawrss(vi);
        nNegFit = nNegFit + 1;
    end
end

% end time monitor
et  = toc;
if floor(et/3600)>0,
    fprintf(1,'Done [%d hours].\n',ceil(et/3600));
else
    fprintf(1,'Done [%d minutes].\n',ceil(et/60));
end;
fprintf(1,'[%s]:Removed negative fits: %d (%.1f%%).\n',...
    mfilename,nNegFit,nNegFit./numel(wProcess).*100);
return;

%-----------------------------------
% make sure that the pRF can only be moved "step" away from original
% position "startParams" - for the one Gaussian model
function [C, Ceq]=distanceCon(x,startParams,step)
Ceq = [];
dist = x([1 2])-startParams([1 2]);
C = norm(dist) - step;
return;
%-----------------------------------

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function state = dss_create_state(X_or_state, parameters)
% Initializes DSS state structure based on existing state
% and given parameters.
%   state = dss_create_state(X or state, parameters)
%     X           Mixed source, matrix with row vectors as signals.
%     state       New or updated persistent algorithm state structure
%     parameters  Optional parameters that are included in the state.
%                 Overwrites parameters in the state variable.
%
% Description of parameter and state variables:
%   algorithm       DSS algorithm
%                     'defl' = Deflation approach (default)
%		      'symm' = Symmetric approach
%		      'pca'  = Linear PCA DSS
%   sdim            Number of components to be extracted
%   wdim            Dimension of the data after whitening (preprocessing)
%   verbose         Verbosity
%                     0 = no print output
%                     1 = terse print output (default)
%                     2 = more detailed print output
%                     3 = verbose print output (shows iteration
%                     progress)
%
%   preprocf.h      Preprocessing function (default: pre_sphere)
%   preprocf.params Modifiable parameter structure
%   orthof.h        Orthogonalization function (default: ortho_default)
%   orthof.params   Modifiable parameter structure
%   denf.h          Denoising function (default: denoise_fica_tanh)
%   denf.params     Modifiable parameter structure
%   stopf.h         Stopping criteria function (default: default_stop)
%   stopf.params    Modifiable parameter structure
%
%   reportf         Vector or list of iteration reporting functions.
%   report          Reporting data used by reporting functions
%
%   alphaf.h        Denoise normalization function (default:none)
%   alphaf.params   Modifiable parameter structure
%   betaf.h         Denoise spectral shift function (default: none)
%   betaf.params    Modifiable parameter structure
%   gammaf.h        Adaptive learning rate function (default:none)
%   gammaf.params   Modifiable parameter structure
%
%   alpha           Constant value for alpha (not adaptive, default: 1)
%   beta            Constant value for beta (not adaptive, default: 0)
%   gamma           Constant value for gamma (not adaptive, default: 1)
%
% Internal data and state variables
%   X          Original mixed source
%   Y          Sphered mixed source
%   A          Mixing matrix
%   B          Unmixing matrix (for X)
%   W          Matrix of all/found projection vectors (for Y)
%   s          Estimate of the current source component (deflation)
%   w          Estimate of the current projection vector (deflation)
%   w_old      Projection vector of the previous iteration (deflation)
%   S          Estimate of all the signal components (symmetric)
%   W_old      Projection vector matrix of the previous iteration (symmetric)
%   iteration  Current iteration number
%   component  Index of currently calculated component (deflation)
%   V          Sphering matrix
%   dV         De-sphering matrix
%
% There are three ways to give values for alpha, beta and gamma:
%  1) Constant value (alpha, beta, gamma),
%  2) function handle (alphaf, betaf, gammaf) with
%    a) constant value (adapt_... = 0) or
%    b) adaptive value (adapt_... = 1)

% Copyright (C) 2004, 2005 DSS MATLAB package team (dss@cis.hut.fi).
% Distributed by Laboratory of Computer and Information Science,
% Helsinki University of Technology. http://www.cis.hut.fi/projects/dss/.
% $Id$

if nargin<1
  error('DSS initialization requires state structure with input data.');
end

% -- Determine calling type
if isstruct(X_or_state)
  state = X_or_state;
else
  state.X = X_or_state;
end

% -- default functions
def_denoise=struct('h', @denoise_fica_tanh,'params', '');
def_stop   =struct('h', @default_stop, ...
    'params', struct('epsilon', 1e-1,'maxiters',200));
def_preproc=struct('h', @pre_sphere);
def_ortho  =struct('h', @ortho_default,  'params', '');
def_report =struct('data', {{}});

% -- parameter & state variables
% Strcture definition: name, description, [default value], [allowed values]
param_def = {
  {'verbose',      'Output verbosity', 1, {0,1,2,3}},
  {'algorithm',    'DSS algorithm type', 'defl', {'defl', 'symm'}},
  {'preprocf',     'Preprocessing function', def_preproc},
  {'orthof',       'Orthogonalization function', def_ortho},
  {'denf',         'Denoising function', def_denoise},
  {'stopf',        'Stoppping criterion function', def_stop},
  {'reportf',      'Vector of iteration reporting functions'},
  {'sdim',         'Source signal dimension'},
  {'wdim',         'Projection dimension'},
  {'report',       'Reporting data', def_report},

  {'alphaf',      'Denoise normalization function'},
  {'betaf',       'Denoise spectral shift function'},
  {'gammaf',      'Adaptive gain function'},
  {'alpha',       'Fixed value for alpha if alphaf is not defined,',1},
  {'beta',        'Fixed value for beta if betaf is not defined,',0},
  {'gamma',       'Fixed value for gamma if gammaf is not defined,',1},
  {'W',           'Initial guess for projection matrix'}
  {'S',           'Initial guess for components'}
  {'Y',           'Pre-sphered data'}
};

if nargin<2
  params = struct([]);
else
  if length(parameters)==0
    params = struct([]);
  else
    if isstruct(parameters)
      dss_message(state,1,'Parse parameter struct\n');
      params = parameters;
    else
      dss_message(state,1,'Parse parameter list\n');
      if mod(size(parameters,2),2)~=0
        error('Parameter list must contain parameter and value pairs.');
      end

      for a=1:2:size(parameters,2)
        parameter = parameters{a};
        if ~ischar(parameter)
          error(['Invalid parameter name at ' num2str(a+1)]);
        end
        value = parameters{a+1};
        % for ML65+: params.(parameter) = value;
        params = setfield(params, parameter, value);
      end
    end
  end
end

% -- Initialize given parameters to state variable
for p=1:size(param_def,1)
  paramname = param_def{p}{1};
  description = param_def{p}{2};
  default_value=[];
  if (size(param_def{p},2)>=3); default_value = param_def{p}{3}; end
  allowed_values=[];
  if (size(param_def{p},2)>=4); allowed_values = param_def{p}{4}; end

  % -- copy given parameter value to state
  if isfield(params, paramname)
    % for ML65+: state.(paramname)=params.(paramname);
    state = setfield(state, paramname, getfield(params, paramname));
    %dss_message(state,1,['Parameter ' paramname '\n']);
    % Clear parameter from struct so that invalid parameters are detected
    params = rmfield(params, paramname);
  end

  % -- default value
  if ~isfield(state, paramname) & ~isempty(default_value)
    % for ML65+: state.(paramname) = default_value;
    state = setfield(state, paramname, default_value);
  end

  % TODO: check allowed values

  % -- check function parameters
  if isfield(state, paramname)
    if isfield(getfield(state, paramname), 'h') & ~isfield(getfield(state, paramname), 'params')
      %dss_message(state,1, sprintf('Creating empty parameter structure for parameter %s.\n', paramname));

      % for ML65+: state.(paramname).params=[];
      p = getfield(state, paramname);
      p.params = [];
      state = setfield(state, paramname, p);

    end
  end

  % -- Print info about variable
  if ~isfield(state, paramname) str = '-';
  else str = tostring(getfield(state, paramname));
  end
  dss_message(state, 2, ['  Parameter ' paramname ': ' str '\n']);

end

% -- Check invalid parameters
% TODO: should invalid parameters throw an error?

% for ML65+: fields  = fieldnames(params);
if isa(params, 'struct') fields  = fieldnames(params);
else fields = [];
end

if size(fields,1)>0
  for p = size(fields, 1)
    dss_message(state,1,['IGNORING INVALID PARAMETER: ''' fields{p} '''\n']);
  end
end

% other state initializations
if ~isfield(state, 'report')
  state.report = [];
end

repo: fieldtrip


text:
function colors = prtPlotUtilClassColorsBW(nClasses)
% Internal function, default color specs.
% xxx Need Help xxx

% colors = [55  126 184; % Blue
%           228 26  28;  % Red
%           77  175 74;  % Green
%           255 127 0;   % Orange
%           153 153 153; % Gray
%           166 86  40;  % Brown
%           152 78  163; % Purple
%           255 255 51;  % Yello
%           247 129 191; % Pink
%           ]/255;

%255, 255, 191;
colors = [215,  25,  28;
          253, 174,  97;
          171, 221, 164;
           43, 131, 186;]/255;

% colors = gray(5);
% colors = colors(1:4,:);

colorMapInd = repmat((1:size(colors,1))',ceil(nClasses/size(colors,1)),1);
colorMapInd = colorMapInd(1:nClasses);

colors = colors(colorMapInd,:);

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function classSymbols = prtPlotUtilClassSymbolsBW(nClasses)
% Internal function, default color specs.
% xxx Need Help xxx

classSymbols = 'os^dv';

classSymbolsInd = repmat((1:length(classSymbols))',ceil(nClasses/length(classSymbols)),1);
classSymbolsInd = classSymbolsInd(1:nClasses);

classSymbols = classSymbols(classSymbolsInd);


text:
function consistent = swiftcc(S, rev, varargin)
% swiftcc is an even faster version of fastcc
%
% USAGE:
%
%    consistent = swiftcc(S, rev [, solver])
%
% INPUTS:
%    S:      the associated sparse stoichiometric matrix
%    rev:    the 0-1 vector with 1's corresponding to the reversible reactions
%
% OPTIONAL INPUT:
%    solver:    the LP solver to be used; the currently available options are
%               'gurobi', 'linprog', and 'cplex' with the default value of
%               'linprog'. It fallbacks to the COBRA LP solver interface if
%               another supported solver is called.
%
% OUTPUT:
%    consistent:    the 0-1 indicator vector of the reactions constituting
%                   the maximum flux consistent metabolic subnetwork
%
% .. Authors:
%       - Mojtaba Tefagh, Stephen P. Boyd, 2019, Stanford University

    [m, n] = size(S);
    consistent = true(n, 1);

    %% setting up the LP solver
    if ~isempty(varargin)
        solver = varargin{1};
    else
        solver = 'linprog';
    end

    %% identifying the blocked irreversible reactions
    result = blocked(S, rev, solver);
    consistent(result.x(m+1:end) < -0.5) = false;

    %% setting up the zero-tolerance parameter
    tol = norm(S(:, consistent), 'fro')*eps(class(S));

    %% identifying the blocked reversible reactions
    [Q, R, ~] = qr(transpose(S(:, consistent)));
    Z = Q(rev(consistent) == 1, sum(abs(diag(R)) > tol)+1:end);

    %% finding the consistent reactions of the original metabolic network
    consistent(consistent & rev == 1) = diag(Z*Z.') > tol^2;
    consistent = find(consistent);
end

author: opencobra

repo: cobratoolbox


text:
function C1 = C1f(epsi)
%C1F  Evaluate C_{1,k}
%
%   C1 = C1F(EPSI) evaluates C_{1,l} using Eq. (18).  EPSI is a K x 1
%   array and C1 is a K x 6 array.

  nC1 = 6;
  C1 = zeros(length(epsi), nC1);
  eps2 = epsi.^2;
  d = epsi;
  C1(:,1) = d.*((6-eps2).*eps2-16)/32;
  d = d.*epsi;
  C1(:,2) = d.*((64-9*eps2).*eps2-128)/2048;
  d = d.*epsi;
  C1(:,3) = d.*(9*eps2-16)/768;
  d = d.*epsi;
  C1(:,4) = d.*(3*eps2-5)/512;
  d = d.*epsi;
  C1(:,5) = -7*d/1280;
  d = d.*epsi;
  C1(:,6) = -7*d/2048;
end

repo: mcbench-benchmarks


text:
filename='BridgeCool_Quadrilateral_Bilinear_Structured';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance','perimeter'};
weights = [1 0.1];
constraint = {'volume'};
optimizer = 'MMA'; incrementFactor = 1;
filterType = 'P1';
constraint_case = 'INEQUALITY';
line_search_initiator = 'STANDARD';

nsteps = 15;
Vfrac_final = 0.2;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-3;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
%% FUNCTION Logistic_iMSF
%   Incompletet Multi-Source Learning with Logistic Loss (interface). This
%   code generates structural information for MultiSource_LogisticR and
%   performs the optimization.
%
%% OBJECTIVE
%   see manual.
%
%% INPUT
%   A_Set: {n * d} * t - input. The cell array of data. Each cell should be
%          n by p_i, with missing samples denoted by a whole row of NaNs.
%   Y_set: {n * 1} * t - output. Y_i \in {-1, 1}.
%   z:     group sparsity regularization parameter (a relative value, [0,1])
%
%% OUTPUT
%   Sol: {struct} solution struct array.
%         For task i: model: struct{i}.x, bais: struct{i}.c
%   funVal: objective function values at each iteration.
%
%% 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 Lei Yuan, 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 12, 2012.
%
%% RELATED PAPERS
%   [1] Lei Yuan, Yalin Wang, Paul M. Thompson, Vaibhav A. Narayan and Jieping
%       Ye, Multi-Source Learning for Joint Analysis of Incomplete
%       Multi-Modality Neuroimaging Data, KDD 2012
%   [2] Lei Yuan, Yalin Wang, Paul M. Thompson, Vaibhav A. Narayan and Jieping
%       Ye, for the Alzheimer's Disease Neuroimaging Initiative, Multi-source
%       Feature Learning for Joint Analysis of Incomplete Multiple Heterogeneous
%       Neuroimaging Data, NeuroImage 2012 Jul 2; 61(3):622-632.
%
%% RELATED FUNCTIONS
%   init_opts, Construct_iMSF, MultiSource_LogisticR

function [Sol, funVal] = Logistic_iMSF(X_Set, Y, lambda, opts)

if nargin< 4
    opts = [];
end

[A_Set, Y_Set, W, G, ind] = Construct_iMSF(X_Set, Y);
[Sol, funVal] = MultiSource_LogisticR(A_Set, Y_Set, lambda, G, W, ind, opts);

author: jiayuzhou

repo: MALSAR

save_path: github-repos/MATLAB/jiayuzhou-MALSAR


text:
function y = BasisQuadPhiX(x1,x2)
y = [x1*x1;
    x1*x2;
    x2*x2;
    x1*x1*x1;
    x1*x1*x2;
    x1*x2*x2;
    x2*x2*x2;
    x1*x1*x1*x1;
    x1*x1*x1*x2;
    x1*x1*x2*x2;
    x1*x2*x2*x2;
    x2*x2*x2*x2;];
end

author: yu-jiang

repo: radpbook


text:
%compute tail-inflection direction and correct trx if necessary
function [data,units]=compute_tailinflang(trx,n)

larvae=trx.exp2flies{n};
numlarvae=numel(larvae);
tailinflang=cell(1,numlarvae);
for i=1:numlarvae
    larva=larvae(i);
     tailinflang{1,i}=bsxfun(@atan2,trx(larva).yinflection_mm-trx(larva).ytail_mm,trx(larva).xinflection_mm-trx(larva).xtail_mm);

end
units=parseunits('rad');
data=tailinflang;

author: kristinbranson

repo: JAABA