text:
function [X,info] = IRsirt(A,b,varargin)
%IRsirt Simuletaneous Iterative Reconstruction Technique
%
% options  = IRsirt('defaults')
% [X,info] = IRsirt(A,b,K)
% [X,info] = IRsirt(A,b,K,options)
%
% This function calls the function 'sirt' in AIR Tools II that implements a
% number of Simultaneous Iterative Reconstruction Technique (SIRT) methods:
% cav, cimmino, drop, landweber and sart (see AIR Tools II for more details).
% The default method is sart which, in the CT community, is known as 'sirt'.
%
% With 'defaults' as input prints the default options.  Otherwise outputs
% the iterates specified in K, using max(K) as MaxIter, and using all other
% default options.  With options as input: uses the user- specified options
% and all the other default options.
%
% Inputs:
%  A : either (a) a full or sparse matrix
%             (b) a matrix object that performs the matrix*vector operation
%             (c) user-defined function m-file
%  b : right-hand side vector
%  K : (optional) integer vector that specifies which iterates are returned
%      in X; the maximum number of iterations is assumed to be max(K)
%      [ positive integer | vector of positive components ]
%  options : structure with the following fields (optional)
%      sirt_mthod   - the specific SIRT method to be used;
%                      [ 'cav' | 'cimmino' | 'drop' | 'landweber' | {'sart'} ]
%      x0            - initial guess for the iterations; default = zero vector
%                      [ array | {'none'} ]
%      MaxIter       - maximum allowed number of iterations
%                      [ {100} | positive integer ]
%      x_true        - true solution; allows us to returns error norms
%                      with respect to x_true at each iteration
%                      [ array | {'none'} ]
%      stopCrit      - stopping criterion for the iterations
%                      [ {'none'} | 'discrep' ]
%                      Note: 'discrep' requires NoiseLevel and eta
%      NoiseLevel    - norm of noise in rhs divided by norm of rhs
%                      [ {'none'} | nonnegative scalar]
%      eta           - safety factor for the discrepancy principle
%                      [ {1.01} | scalar greater than (and close to) 1 ]
%      relaxParam    - constant relaxation parameter, bounded above 2 divided
%                      by the spectral radius of the iteration matrix
%                      [ positive scalar | {'none'} ]
%                      If 'none' then a good value is chosed by the function
%      nonnegativity - apply nonnegativity constraints
%                      [ 'on' | {'off'} ]
%      Ubound        - apply box constraints in the interval [0,Ubound]
%                      [ positive scalar | {'off'} ]
%      IterBar       - shows the progress of the iterations
%                      [ {'on'} | 'off' ]
% Note: the options structure can be created using the function IRset.
%
% Outputs:
%   X : computed solutions, stored column-wise (at the iterations listed in K)
%   info : structure with the following fields:
%      its              - number of the last computed iteration
%      saved_iterations - iteration numbers of iterates stored in X
%      StopFlag         - a flag that describes the stopping condition:
%                           1 : reached maximum number of iterations
%                           2 : discrepancy principle satisfied
%      relaxParam       - the used relaxation parameter
%      Enrm             - relative error norms (requires x_true) for the
%                         stored iterations
%
% Note that this function provides a simplified call to the function "sirt"
% in AIR Tools II which has more features than we allow here.  To use the
% full power of the SIRT methods consider using AIR Tools II directly.
%
% See also: IRart, IRget, IRset

% Silvia Gazzola, University of Bath
% Per Christian Hansen, Technical University of Denmark
% James G. Nagy, Emory University
% April, 2018.

% This file is part of the IR Tools package and is distributed under the
% 3-Clause BSD License. A separate license file should be provided as part
% of the package.

% Set default values for options.
defaultopt = struct('x0', 'none', 'IterBar', 'on', 'stopCrit', 'none', ...
    'eta', 1.01, 'nonnegativity', 'off', 'NoiseLevel','none', ...
    'Ubound', 'none', 'relaxParam', 'none', 'sirt_method', 'sart', ...
    'MaxIter', 100, 'x_true', 'none');

% If input is 'defaults,' return the default options in X.
if nargin==1 && nargout <= 1 && isequal(A,'defaults')
    X = defaultopt;
    return;
end

% Check for acceptable number of optional input arguments.
switch length(varargin)
    case 0
        K = []; options = [];
    case 1
        if isa(varargin{1}, 'double')
            K = varargin{1}; options = [];
        else
            K = []; options = varargin{1};
        end
    case 2
        if isa(varargin{1}, 'double')
            K = varargin{1}; options = varargin{2};
        else
            K = varargin{2}; options = varargin{1};
        end
    otherwise
        error('Too many input parameters')
end

if isempty(options)
    options = defaultopt;
end

options = IRset(defaultopt, options);

if min(K) < 1, error('Number of iterations must be positive'), end

x0 = IRget(options, 'x0', [], 'fast');
if strcmpi(x0, 'none')
    x0 = [];
end

MaxIter       = IRget(options, 'MaxIter',       [], 'fast');
method        = IRget(options, 'sirt_method',   [], 'fast');
IterBar       = IRget(options, 'IterBar',       [], 'fast');
relaxParam    = IRget(options, 'relaxParam',    [], 'fast');
stopCrit      = IRget(options, 'stopCrit',      [], 'fast');
nonnegativity = IRget(options, 'nonnegativity', [], 'fast');
Ubound        = IRget(options, 'Ubound',        [], 'fast');
x_true        = IRget(options, 'x_true',        [], 'fast');

if isempty(K)
    K = MaxIter;
end
% Sorting the iteration numbers (in case they are shuffled in input).
K = sort(K,'ascend'); K = unique(K);
if ~((isreal(K) && (all(K > 0)) && all(K == floor(K))))
    error('K must be a vector of positive real integers')
end

% Set options for SIRT.  Always used a fixed lambda, either chosen by the
% user or set by the SIRT function.  Always use no stopping rule or stop
% by the size of the residual (discrepancy principle).

if strcmpi(relaxParam,'none')
    % The field lambda is not present.
else
   sirtoptions.relaxpar = relaxParam;
end

if strcmpi(stopCrit,'discrep')
    sirtoptions.stoprule.type = 'DP';
    if isempty(options.NoiseLevel) || strcmp(options.NoiseLevel,'none')
        error('When using discrepancy principle: options.NoiseLevel must be specified')
    end
    sirtoptions.stoprule.taudelta = options.eta*options.NoiseLevel*norm(b);
else
    sirtoptions.stoprule.type = 'none';
end

if strcmpi(nonnegativity,'on')
    sirtoptions.lbound = 0;
end

if (isreal(Ubound) && isscalar(Ubound) && Ubound > 0)
    sirtoptions.lbound = 0;
    sirtoptions.ubound = Ubound;
end

if strcmp(IterBar,'on')
    sirtoptions.waitbar = true;
end

% Call the requested SIRT method.
[X,sirtinfo] = sirt(method,A,b,K,x0,sirtoptions);

switch sirtinfo.stoprule
    case 0
        info.StopFlag = 1;
    case 2
        info.StopFlag = 2;
        if sirtinfo.finaliter==0
            warning(['No iterations were performed, beause the starting',...
                     ' vector x0 satisfies the discrepancy principle'])
        end
end
info.its = sirtinfo.finaliter;
info.saved_iterations = sirtinfo.itersaved;
info.relaxParam = sirtinfo.relaxpar;

% Compute relative error norms, if requested.
if ~strcmp(x_true,'none')
    Enrm = zeros(length(info.saved_iterations),1);
    for k=1:length(info.saved_iterations)
        Enrm(k) = norm(x_true-X(:,k));
    end
    info.Enrm = Enrm/norm(x_true);
end

repo: IRtools


text:
% F1=figure;
% set(F1,'position',[18 144 560 420]);%Numero de paralelepipedos
lo=2;               %Lado del paralelepipedo base
deltax=5;           %Altura del paralelepipedo
xo=-lo/2;           %Punto de origen en x
yo=-lo/2;           %Punto de origen en y
zo=0;               %Punto de origen en z

line([0 0],[0 0],[0 10],'color','r')
line([0 5],[0 0],[5 5],'color','g')
line([0 0],[0 5],[5 5],'color','b')

x1=xo;
x2=x1+lo;
x3=x1+lo;
x4=x1;
x5=x1;
x6=x1;
x7=x1+lo;
x8=x1+lo;

y1=yo;
y2=y1;
y3=y1+lo;
y4=y1+lo;
y5=y1+lo;
y6=y1;
y7=y1;
y8=y1+lo;

z1=zo;
z2=zo;
z3=zo;
z4=zo;
z5=zo+deltax;
z6=zo+deltax;
z7=zo+deltax;
z8=zo+deltax;

x=[x1 x2 x3 x4 x5 x6 x7 x8];
y=[y1 y2 y3 y4 y5 y6 y7 y8];
z=[z1 z2 z3 z4 z5 z6 z7 z8];

%///////////////////////////////////////////////////
%//////////////////////////////////////////////////
for n=1:8
T(1:4,n)=[cos(teta) -sin(teta) 0 0
   sin(teta) cos(teta) 0 0
   0 0 1 0
   0 0 0 1]*[x(n);y(n);z(n);1];
end

a=0.4;
tcolor(1,1,1:3) = [0.8 0.4 0.3];
A1=patch([T(1,1) T(1,2) T(1,3) T(1,4)],[T(2,1) T(2,2) T(2,3) T(2,4)],[T(3,1) T(3,2) T(3,3) T(3,4)],tcolor);
A2=patch([T(1,1) T(1,2) T(1,7) T(1,6)],[T(2,1) T(2,2) T(2,7) T(2,6)],[T(3,1) T(3,2) T(3,7) T(3,6)],tcolor);
A3=patch([T(1,5) T(1,6) T(1,7) T(1,8)],[T(2,5) T(2,6) T(2,7) T(2,8)],[T(3,5) T(3,6) T(3,7) T(3,8)],tcolor);
A4=patch([T(1,2) T(1,3) T(1,8) T(1,7)],[T(2,2) T(2,3) T(2,8) T(2,7)],[T(3,2) T(3,3) T(3,8) T(3,7)],tcolor);
A5=patch([T(1,1) T(1,4) T(1,5) T(1,6)],[T(2,1) T(2,4) T(2,5) T(2,6)],[T(3,1) T(3,4) T(3,5) T(3,6)],tcolor);
A6=patch([T(1,3) T(1,4) T(1,5) T(1,8)],[T(2,3) T(2,4) T(2,5) T(2,8)],[T(3,3) T(3,4) T(3,5) T(3,8)],tcolor);
%
% light('Position',[1 0 0],'Style','infinite');
% light('Position',[0 1 0],'Style','infinite');
% light('Position',[0 0 1],'Style','infinite');
% light('Position',[-1 0 0],'Style','infinite');
% light('Position',[0 -1 0],'Style','infinite');
% % light('Position',[0 0 -1],'Style','infinite');

deltax2=deltax;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%ESLABON 2
lo=1.5;               %Lado del paralelepipedo base
deltax=1.5;           %Altura del paralelepipedo
l3=2;                 %Longitud del eslabon 3
xo=-1;                %Punto de origen en x
yo=1;                 %Punto de origen en y
zo=-1.5/2;            %Punto de origen en z
l2=7;
x1=xo;
x2=x1+l2;
x3=x1+l2;
x4=x1;
x5=x1;
x6=x1;
x7=x1+l2;
x8=x1+l2;

y1=yo;
y2=y1;
y3=y1+lo;
y4=y1+lo;
y5=y1+lo;
y6=y1;
y7=y1;
y8=y1+lo;

z1=zo;
z2=zo;
z3=zo;
z4=zo;
z5=zo+deltax;
z6=zo+deltax;
z7=zo+deltax;
z8=zo+deltax;
x=[x1 x2 x3 x4 x5 x6 x7 x8];
y=[y1 y2 y3 y4 y5 y6 y7 y8];
z=[z1 z2 z3 z4 z5 z6 z7 z8];

for n=1:8

T(1:4,n)=[cos(teta) -sin(teta) 0 0
   sin(teta) cos(teta) 0 0
   0 0 1 deltax2-deltax/2
   0 0 0 1]*[cos(teta2) 0 sin(teta2) 0
          0 1 0 0
          -sin(teta2) 0 cos(teta2) 0
          0 0 0 1]*[x(n);y(n);z(n);1];

% Ta(1:4,n)=[cos(teta) -sin(teta) 0 0
%    sin(teta) cos(teta) 0 0
%    0 0 1 0
%    0 0 0 1]*[x(n);y(n);z(n);1];
end

a=0.4;
tcolor(1,1,1:3) = [0.4 0.1 0.1 ];
A7=patch([T(1,1) T(1,2) T(1,3) T(1,4)],[T(2,1) T(2,2) T(2,3) T(2,4)],[T(3,1) T(3,2) T(3,3) T(3,4)],tcolor);
A8=patch([T(1,1) T(1,2) T(1,7) T(1,6)],[T(2,1) T(2,2) T(2,7) T(2,6)],[T(3,1) T(3,2) T(3,7) T(3,6)],tcolor);
A9=patch([T(1,5) T(1,6) T(1,7) T(1,8)],[T(2,5) T(2,6) T(2,7) T(2,8)],[T(3,5) T(3,6) T(3,7) T(3,8)],tcolor);
A10=patch([T(1,2) T(1,3) T(1,8) T(1,7)],[T(2,2) T(2,3) T(2,8) T(2,7)],[T(3,2) T(3,3) T(3,8) T(3,7)],tcolor);
A11=patch([T(1,1) T(1,4) T(1,5) T(1,6)],[T(2,1) T(2,4) T(2,5) T(2,6)],[T(3,1) T(3,4) T(3,5) T(3,6)],tcolor);
A12=patch([T(1,3) T(1,4) T(1,5) T(1,8)],[T(2,3) T(2,4) T(2,5) T(2,8)],[T(3,3) T(3,4) T(3,5) T(3,8)],tcolor);
% light('Position',[1 0 0],'Style','infinite');
% light('Position',[0 1 0],'Style','infinite');
% light('Position',[0 0 1],'Style','infinite');
% light('Position',[-1 0 0],'Style','infinite');
% light('Position',[0 -1 0],'Style','infinite');
% light('Position',[0 0 -1],'Style','infinite');

grid on

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%//////Eslabon 3 //////////////////////////////////////////////////

lo=1;                   %Lado del paralelepipedo base
deltax=1;               %Altura del paralelepipedo
xo=-deltax/2;           %Punto de origen en x
yo=0;                   %Punto de origen en y
zo=-.5;                 %Punto de origen en z
l2=5.5;
x1=xo;
x2=x1+l2;
x3=x1+l2;
x4=x1;
x5=x1;
x6=x1;
x7=x1+l2;
x8=x1+l2;

y1=yo;
y2=y1;
y3=y1+lo;
y4=y1+lo;
y5=y1+lo;
y6=y1;
y7=y1;
y8=y1+lo;

z1=zo;
z2=zo;
z3=zo;
z4=zo;
z5=zo+deltax;
z6=zo+deltax;
z7=zo+deltax;
z8=zo+deltax;
x=[x1 x2 x3 x4 x5 x6 x7 x8];
y=[y1 y2 y3 y4 y5 y6 y7 y8];
z=[z1 z2 z3 z4 z5 z6 z7 z8];

for n=1:8
T(1:4,n)=[cos(teta) -sin(teta) 0 0
   sin(teta) cos(teta) 0 0
   0 0 1 deltax2-deltax/2
   0 0 0 1]*[cos(teta2) 0 sin(teta2) 0
          0 1 0 0
          -sin(teta2) 0 cos(teta2) 0
          0 0 0 1]*[cos(teta3) 0 sin(teta3) 5+deltax
          0 1 0 0
          -sin(teta3) 0 cos(teta3) 0
          0 0 0 1]*[x(n);y(n);z(n);1];
end

crgb=[249 193 120]/255;
tcolor(1,1,1:3) = [crgb(1) crgb(2) crgb(3)];

A13=patch([T(1,1) T(1,2) T(1,3) T(1,4)],[T(2,1) T(2,2) T(2,3) T(2,4)],[T(3,1) T(3,2) T(3,3) T(3,4)],tcolor);
A14=patch([T(1,1) T(1,2) T(1,7) T(1,6)],[T(2,1) T(2,2) T(2,7) T(2,6)],[T(3,1) T(3,2) T(3,7) T(3,6)],tcolor);
A15=patch([T(1,5) T(1,6) T(1,7) T(1,8)],[T(2,5) T(2,6) T(2,7) T(2,8)],[T(3,5) T(3,6) T(3,7) T(3,8)],tcolor);
A16=patch([T(1,2) T(1,3) T(1,8) T(1,7)],[T(2,2) T(2,3) T(2,8) T(2,7)],[T(3,2) T(3,3) T(3,8) T(3,7)],tcolor);
A17=patch([T(1,1) T(1,4) T(1,5) T(1,6)],[T(2,1) T(2,4) T(2,5) T(2,6)],[T(3,1) T(3,4) T(3,5) T(3,6)],tcolor);
A18=patch([T(1,3) T(1,4) T(1,5) T(1,8)],[T(2,3) T(2,4) T(2,5) T(2,8)],[T(3,3) T(3,4) T(3,5) T(3,8)],tcolor);
%
% light('Position',[1 0 0],'Style','infinite');
% light('Position',[0 1 0],'Style','infinite');
% light('Position',[0 0 1],'Style','infinite');
% light('Position',[-1 0 0],'Style','infinite');
% light('Position',[0 -1 0],'Style','infinite');
% light('Position',[0 0 -1],'Style','infinite');

lo=10;                %Lado del paralelepipedo base
deltax=1;             %Altura del paralelepipedo
xo=-lo/2;             %Punto de origen en x
yo=-lo/2;             %Punto de origen en y
zo=-deltax;           %Punto de origen en z

x1=xo;
x2=x1+lo;
x3=x1+lo;
x4=x1;
x5=x1;
x6=x1;
x7=x1+lo;
x8=x1+lo;

y1=yo;
y2=y1;
y3=y1+lo;
y4=y1+lo;
y5=y1+lo;
y6=y1;
y7=y1;
y8=y1+lo;

z1=zo;
z2=zo;
z3=zo;
z4=zo;
z5=zo+deltax;
z6=zo+deltax;
z7=zo+deltax;
z8=zo+deltax;
a=0.4;
tcolor(1,1,1:3) = [0.3 0.5 0.4];
patch([x1 x2 x3 x4],[y1 y2 y3 y4],[z1 z2 z3 z4],tcolor);
patch([x1 x2 x7 x6],[y1 y2 y7 y6],[z1 z2 z7 z6],tcolor);
patch([x5 x6 x7 x8],[y5 y6 y7 y8],[z5 z6 z7 z8],tcolor);
patch([x2 x3 x8 x7],[y2 y3 y8 y7],[z2 z3 z8 z7],tcolor);
patch([x1 x4 x5 x6],[y1 y4 y5 y6],[z1 z4 z5 z6],tcolor);
patch([x3 x4 x5 x8],[y3 y4 y5 y8],[z3 z4 z5 z8],tcolor);
%
% light('Position',[1 0 0],'Style','infinite');
% light('Position',[0 1 0],'Style','infinite');
% light('Position',[0 0 1],'Style','infinite');
% light('Position',[-1 0 0],'Style','infinite');
% light('Position',[0 -1 0],'Style','infinite');
% light('Position',[0 0 -1],'Style','infinite');

p=[0;0;0;1];

Tf=[cos(teta) -sin(teta) 0 0
   sin(teta) cos(teta) 0 0
   0 0 1 deltax2-deltax/2
   0 0 0 1]*[cos(teta2) 0 sin(teta2) 0
          0 1 0 0
          -sin(teta2) 0 cos(teta2) 0
          0 0 0 1]*[cos(teta3) 0 sin(teta3) 5+deltax
          0 1 0 0
          -sin(teta3) 0 cos(teta3) 0
          0 0 0 1]*p
Tf(2)=T(2,8)
Tf(3)=T(3,8)
Tf(1)=T(1,8)

E1=gca;
axis([-10 8 -10 8 0 15])
hold on
plot3(E1,p(1),p(2),p(3));
grid on

author: Sable

repo: mcbench-benchmarks


text:
function [ a1, a2 ] = i4vec_swap ( n, a1, a2 )

%*****************************************************************************80
%
%% I4VEC_SWAP swaps the entries of two I4VEC's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries in the arrays.
%
%    Input, integer A1(N), A2(N), the vectors to swap.
%
%    Output, integer A1(N), A2(N), the swapped vectors.
%
  a3(1:n) = a1(1:n);
  a1(1:n) = a2(1:n);
  a2(1:n) = a3(1:n);

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [F,V]=mesh2tri(X,Y,Z,tri_type)

% function [F,V]=mesh2tri(X,Y,Z,tri_type)
% ------------------------------------------------------------------------
%
% This function converts a regular mesh defined by X,Y and Z into a regular
% triangulation. The output is patch data (triangles) in the faces �F� and
% vertices �V� format. The quadrilateral mesh faces are converted to
% triangles by splitting the faces into triangles according to the setting
% tri_type:
%   tri_type ='f' -> forward slash division of quadrilateral
%   tri_type ='b' -> back slash division of quadrilateral
%   tri_type ='x' -> Cross division of quadrilateral
%
% The output coordinates "V" are in the form of V=[X(:),Y(:),Z(:)];
% For forward and back slash subdivision no extra coordinates are
% introduced and therefore the original meshgrid formatted coordinates can
% still be used for plotting, see examples below.
% For cross division extra points are created at the centre of each
% quadrilateral face using the mean of the input coordinates. The extra
% coordinates are the last prod(size(X)-1) points (e.g.
% V((numel(X)+1):end,:) )and can therefore be replaced by interpolated
% coordinates if desired, see example.
%
%
% %% EXAMPLE
% clear all; close all; clc;
%
% [X,Y] = meshgrid(linspace(-10,10,25));
% Z = sinc(sqrt((X/pi).^2+(Y/pi).^2));
%
% figure('units','normalized','Position',[0 0 1 1],'Color','w'); colordef('white');
% subplot(2,2,1);
% surf(X,Y,Z); hold on;
% axis tight; axis square; grid on; axis off; view(3); view(-30,70);
% title('Meshgrid','FontSize',20);
%
% [F,V]=mesh2tri(X,Y,Z,'f');
% C=V(:,3); C=mean(C(F),2);
% subplot(2,2,2);
% patch('Faces',F,'Vertices',V,'FaceColor','flat','CData',C); hold on;
% axis tight; axis square; grid on; axis off; view(3); view(-30,70);
% title('Forward slash','FontSize',20);
%
% [F,V]=mesh2tri(X,Y,Z,'b');
% C=V(:,3); C=mean(C(F),2);
% subplot(2,2,3);
% Example of using original meshgrid coordinates instead
% trisurf(F,X,Y,Z);
% axis tight; axis square; grid on; axis off; axis off; view(3); view(-30,70);
% title('Back slash','FontSize',20);
%
% [F,V]=mesh2tri(X,Y,Z,'x');
% Replace Z-coordinates of added points by interpolated values if desired
% IND=(numel(X)+1):size(V,1);
% ZI = interp2(X,Y,Z,V(IND,1),V(IND,2),'cubic');
% V(IND,3)=ZI;
%
% C=V(:,3); C=mean(C(F),2);
% subplot(2,2,4);
% patch('Faces',F,'Vertices',V,'FaceColor','flat','CData',C); hold on;
% axis tight; axis square; grid on; axis off; view(3); view(-30,70);
% title('Crossed','FontSize',20);
%
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% 15/07/2010
%------------------------------------------------------------------------

[J,I]=meshgrid(1:1:size(X,2)-1,1:1:size(X,1)-1);

switch tri_type
    case 'f'%Forward slash
        TRI_I=[I(:),I(:)+1,I(:)+1;  I(:),I(:),I(:)+1];
        TRI_J=[J(:),J(:)+1,J(:);   J(:),J(:)+1,J(:)+1];
        F = sub2ind(size(X),TRI_I,TRI_J);
    case 'b'%Back slash
        TRI_I=[I(:),I(:)+1,I(:);  I(:)+1,I(:)+1,I(:)];
        TRI_J=[J(:)+1,J(:),J(:);   J(:)+1,J(:),J(:)+1];
        F = sub2ind(size(X),TRI_I,TRI_J);
    case 'x'%Cross
        TRI_I=[I(:)+1,I(:);  I(:)+1,I(:)+1;  I(:),I(:)+1;    I(:),I(:)];
        TRI_J=[J(:),J(:);    J(:)+1,J(:);    J(:)+1,J(:)+1;  J(:),J(:)+1];
        IND=((numel(X)+1):numel(X)+prod(size(X)-1))';
        F = sub2ind(size(X),TRI_I,TRI_J);
        F(:,3)=repmat(IND,[4,1]);
        Fe_I=[I(:),I(:)+1,I(:)+1,I(:)]; Fe_J=[J(:),J(:),J(:)+1,J(:)+1];
        Fe = sub2ind(size(X),Fe_I,Fe_J);
        Xe=mean(X(Fe),2); Ye=mean(Y(Fe),2);  Ze=mean(Z(Fe),2);
        X=[X(:);Xe(:)]; Y=[Y(:);Ye(:)]; Z=[Z(:);Ze(:)];
end

V=[X(:),Y(:),Z(:)];

end

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [x1,x2,nS]=bracket(S,x0,d,problem,stepsize)
%   Bracket the minimum (or maximum) of the objective function
%   in the search direction.
%
%   [x1,x2,nS]=bracket(S,x0,d,problem,stepsize)
%
%   S: objective function
%   x0: initial point
%   d: search direction vector
%   problem: (-1): minimum (default), (1): maximum
%   stepsize: initial stepsize (default = 0.01*norm(d))
%   [x1,x2]: unsorted lower and upper limits
%   nS: number of objective function evaluations

%   Copyright (c) 2001 by LASIM-DEQUI-UFRGS
%   $Revision: 1.0 $  $Date: 2001/07/04 21:45:10 $
%   Argimiro R. Secchi (arge@enq.ufrgs.br)

 if nargin < 3,
   error('bracket requires three input arguments');
 end
 if nargin < 4,
   problem=-1;
 end
 if nargin < 5,
   stepsize=0.01*norm(d);
 end

 d=d(:);
 x0=x0(:);
 j=0; nS=1;
 y0=feval(S,x0)*problem;

 while j < 2,
  x=x0+stepsize*d;
  y=feval(S,x)*problem;
  nS=nS+1;

  if y0 >= y,
    stepsize=-stepsize;
    j=j+1;
  else
    while y0 < y,
      stepsize=2*stepsize;
      y0=y;
      x=x+stepsize*d;
      y=feval(S,x)*problem;
      nS=nS+1;
    end
    j=1;
    break;
  end
 end

 x2=x;
 x1=x0+stepsize*(j-1)*d;


text:
function [params, names] = matern32KernExtractParam(kern)

% MATERN32KERNEXTRACTPARAM Extract parameters from the MATERN32 kernel structure.
% FORMAT
% DESC extracts parameters from the matern kernel with nu=3/2
% kernel structure into a vector of parameters for optimisation.
% ARG kern : the kernel structure containing the parameters to be
% extracted.
% RETURN param : vector of parameters extracted from the kernel. If
% the field 'transforms' is not empty in the kernel structure, the
% parameters will be transformed before optimisation (for example
% positive only parameters could be logged before being returned).
%
% FORMAT
% DESC extracts parameters and parameter names from the matern kernel with nu=3/2
% kernel structure.
% ARG kern : the kernel structure containing the parameters to be
% extracted.
% RETURN param : vector of parameters extracted from the kernel. If
% the field 'transforms' is not empty in the kernel structure, the
% parameters will be transformed before optimisation (for example
% positive only parameters could be logged before being returned).
% RETURN names : cell array of strings containing names for each
% parameter.
%
% SEEALSO matern32KernParamInit, matern32KernExpandParam, kernExtractParam, scg, conjgrad
%
% COPYRIGHT : Neil D. Lawrence, 2006
%
% KERN

params = [kern.lengthScale kern.variance];
if nargout > 1
  names={'length scale', 'variance'};
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
classdef LaunchVehicleTankState < matlab.mixin.SetGet & matlab.mixin.Copyable
    %LaunchVehicle Summary of this class goes here
    %   Detailed explanation goes here

    properties
        %validators commented out to improve performance
        stageState = LaunchVehicleStageState.empty(1,0)% LaunchVehicleStageState
        tank% LaunchVehicleTank
        tankMass = 0%(1,1) double = 0; %mT
    end

    methods
        function obj = LaunchVehicleTankState(stageState)
            obj.stageState = stageState;
        end

        function tankMass = getTankMass(obj)
            tankMass = [obj.tankMass];
        end

        function setTankMass(obj, tankMass)
            obj.tankMass = tankMass;
        end

        function newTankState = deepCopy(obj)
%             newTankState = LaunchVehicleTankState(obj.stageState);
            newTankState = obj.copy();

%             newTankState.tank = obj.tank;
%             newTankState.tankMass = obj.tankMass;
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function spm_dcm_search(P)
% Post hoc optimisation of DCMs (under Laplace approximation)
% FORMAT spm_dcm_search(P)
%
% P         -  character/cell array of DCM filenames
%
%--------------------------------------------------------------------------
% spm_dcm_search operates on different DCMs of the same data to identify
% the best model. It will invert the full model whose free-parameters are
% the union (superset) of all free parameters in each model specified. The
% routine then uses a post hoc selection procedure to evaluate the log-
% evidence and conditional density over free-parameters of each model
% specified.
%
% The DCM specified does not need to be estimated. spm_dcm_search will
% invert the requisite (full DCM) automatically.
%
% The outputs of this routine are graphics reporting the model space search
% (optimisation) and a DCM_optimum (in the first DCMs directory) for the
% best DCM. The structural and function (spectral embedding) graphs are
% based on this DCM.
%
% DCM_optimum  contains the fields:
%        DCM.P   - character/cell array of DCM filenames
%        DCM.PF  - their associated free energies
%        DCM.PP  - and posterior (model) probabilities
%
% In addition, the free energies and posterior estimates of each DCM in P
% are saved for subsequent searches over different partitions of model
% space.
%
% See alos: spm_dcm_post_hoc.m
%
%__________________________________________________________________________
% Copyright (C) 2008-2011 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_dcm_search.m 6615 2015-11-30 12:56:02Z peter $

% get filenames
%--------------------------------------------------------------------------
try
    P;
catch
    [P, sts] = spm_select([2 Inf],'^DCM.*\.mat$','Select DCM*.mat files');
    if ~sts, return; end
end

if ischar(P), P = cellstr(P); end
N = numel(P);

%-Check models are compatible in terms of their data
%==========================================================================
for j = 1:N

    % get prior covariances
    %----------------------------------------------------------------------
    load(P{j});

    % and compare it with the first model
    %----------------------------------------------------------------------
    if j == 1
        Y = DCM.Y.y;
    else
        try
            if any(any(Y - DCM.Y.y))
                fprintf('Please check model %i for compatibility',j)
                return
            end
        catch
            fprintf('Please check model %i for compatibility',j)
            return
        end
    end

    % accumate model in terms of which parameters are free
    %----------------------------------------------------------------------
    A = DCM.a;
    B = DCM.b;
    C = DCM.c;
    D = DCM.d;

    % Get full models free parameters
    %----------------------------------------------------------------------
    if j == 1
        a    = A;
        b    = B;
        c    = C;
        d    = D;
    else
        a    = a | A;
        b    = b | B;
        c    = c | C;
        d    = d | D;
    end
end

%-Estimate full model
%==========================================================================
DCM.a = a;
DCM.b = b;
DCM.c = c;
DCM.d = d;

DCM.name = 'optimum';

% Get full priors and posteriors
% -------------------------------------------------------------------------
FUL   = spm_dcm_estimate(DCM);

qE    = FUL.Ep;
qC    = FUL.Cp;
pE    = FUL.M.pE;
pC    = FUL.M.pC;

%-Loop through models and get log-evidences
%==========================================================================
for j = 1:N

    load(P{j});

    % Fix for endogenous DCM (to match spm_dcm_estimate)
    % ---------------------------------------------------------------------
    if isempty(DCM.c) || isempty(U.u)
        DCM.c  = zeros(DCM.n,1);
        DCM.b  = zeros(DCM.n,DCM.n,1);
    end

    % Get model (priors) and evaluate (reduced) free-energy and posteriors
    % ---------------------------------------------------------------------
    [rE,rC]   = spm_dcm_fmri_priors(DCM.a,DCM.b,DCM.c,DCM.d,DCM.options);
    [F,Ep,Cp] = spm_log_evidence(qE,qC,pE,pC,rE,rC);

    % Put reduced conditional estimates in DCM
    % =====================================================================

    % Bayesian inference and variance
    %----------------------------------------------------------------------
    sw       = warning('off','SPM:negativeVariance');
    Pp       = spm_unvec(1 - spm_Ncdf(0,abs(spm_vec(Ep)),diag(Cp)),Ep);
    Vp       = spm_unvec(diag(Cp),Ep);
    warning(sw);

    % Store parameter estimates
    %----------------------------------------------------------------------
    DCM.M.pC = rC;
    DCM.Ep   = Ep;
    DCM.Cp   = Cp;
    DCM.Pp   = Pp;
    DCM.Vp   = Vp;
    DCM.T    = 0;

    % Store predictions of states from full model for simplicity
    %----------------------------------------------------------------------
    DCM.Ce   = FUL.Ce;
    DCM.H1   = FUL.H1;
    DCM.K1   = FUL.K1;
    DCM.R    = FUL.R;
    DCM.y    = FUL.y;

    % Save approximations to model evidence: negative free energy, AIC, BIC
    %----------------------------------------------------------------------
    evidence = spm_dcm_evidence(DCM);
    DCM.F    = F;
    DCM.AIC  = evidence.aic_overall;
    DCM.BIC  = evidence.bic_overall;

    % Save DCM
    %======================================================================
    save(P{j},'DCM','F','Ep','Cp', spm_get_defaults('mat.format'));

    % Record free-energy
    %----------------------------------------------------------------------
    G(j) = F;

end

% Model evidences and best model
% =========================================================================
G     = G - min(G);
p     = exp(G - max(G));
p     = p/sum(p);

% Get selected model
%--------------------------------------------------------------------------
[q,j] = max(p);
load(P{j});

i   = spm_fieldindices(DCM.Ep,'A','B','C','D');
qE  = spm_vec(FUL.Ep);
Ep  = spm_vec(DCM.Ep);
qC  = DCM.Cp;
Cp  = DCM.Cp;
F   = DCM.F;

% Show results
% -------------------------------------------------------------------------
spm_figure('Getwin','Graphics'); clf

subplot(2,2,1)
if length(P) > 32, plot(G,'k'), else bar(G,'c'), end
title('log-posterior','FontSize',16)
xlabel('model','FontSize',12)
ylabel('log-probability','FontSize',12)
axis square

subplot(2,2,2)
if length(P) > 32, plot(p,'k'), else, bar(p,'r'), end
title('model posterior','FontSize',16)
xlabel('model','FontSize',12)
ylabel('probability','FontSize',12)
axis square

% Show full and reduced conditional estimates (for optimum DCM)
%--------------------------------------------------------------------------
subplot(2,2,3)
spm_plot_ci(qE(i),qC(i,i))
title('MAP connections (full)','FontSize',16)
axis square
a   = axis;

subplot(2,2,4)
spm_plot_ci(Ep(i),Cp(i,i))
title('MAP connections (optimum)','FontSize',16)
axis square
axis(a)

% Show structural and functional graphs
%--------------------------------------------------------------------------
spm_figure('Getwin','Graph'); clf

spm_dcm_graph(DCM.xY,DCM.Ep.A)

%-Save optimum and full DCM
%==========================================================================
DCM.P  = P;
DCM.PF = G;
DCM.PP = p;

% Reduced model (optimum)
%--------------------------------------------------------------------------
pth      = fileparts(P{1});
filename = fullfile(pth,'DCM_optimum.mat');
save(filename,'DCM','F','Ep','Cp', spm_get_defaults('mat.format'));

% Full model
%--------------------------------------------------------------------------
DCM      = FUL;
Ep       = FUL.Ep;
Cp       = FUL.Cp;
F        = FUL.F;
filename = fullfile(pth,'DCM_full');
save(filename,'DCM','F','Ep','Cp', spm_get_defaults('mat.format'));

author: spm


text:
function [z phi]=sim_rbf(Xc,X,W,k_i,basisfunction)
%simulates a radial basis function
%Xc is a N_r by N_dim matrix of rbf centres
%X is a N_p by N_dim matrix of the points to simulate
%W is the weight vector
%basisfunction may be 'gaussian' or 'polyharmonicspline'
%k_i is a prescaler for 'gaussian' rbf and function order for
%'polyharmonicspline'. See k_i(i)=0 for constant bias

%    Copyright Travis Wiens 2008
%
%    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/>.
%
%    If you would like to request that this software be licensed under a less
%    restrictive license (i.e. for commercial closed-source use) please
%    contact Travis at travis.mlfx@nutaksas.com

if nargin<4
    k_i=1;
end

if nargin<5
    basisfunction='gaussian';
end

N_r=size(Xc,1);%number of rbf centres
N_p=size(X,1);%number of points

if numel(k_i)==1
    k_i=k_i*ones(N_r);
end

phi=zeros(N_p,N_r);%rbf outputs

for i=1:N_r
    if k_i(i)==0
        phi(:,i)=1;
    else
        r=sqrt(sum((repmat(Xc(i,:),N_p,1)-X(:,:)).^2,2));%distance from point Xc to X

        switch basisfunction
            case {'gaussian','Gaussian'}
                phi(:,i)=exp(-k_i(i).*r.^2);
            case {'phs','polyharmonicspline'}
                if r==0
                    phi(:,i)=0;
                else
                    if round(k_i(i)/2)==k_i(i)/2%even
                        phi(:,i)=r.^k_i(i).*log(r);
                        phi(r==0,i)=0;%avoid log(0)
                    else
                        phi(:,i)=r.^k_i(i);
                    end
                end
            otherwise
                error('unknown basis function')
        end
    end
end

z=phi*W;


text:
function x_stdft = stdft(x, N, K, N_fft)

frames      = 1:K:(length(x)-N);
x_stdft     = zeros(length(frames), N_fft);

w           = ml_hanning(N);
x           = x(:);

for i = 1:length(frames)
    ii              = frames(i):(frames(i)+N-1);
	x_stdft(i, :) 	= fft(x(ii).*w, N_fft);
end

author: mpariente

repo: pystoi

save_path: github-repos/MATLAB/mpariente-pystoi


text:
function [y1, y2] = vl_nnpdist(x, x0, p, varargin)
%VL_NNPDIST CNN p-distance from target.
%   VL_NNPDIST(X, X0, P) computes the P distance raised of each feature
%   vector in X to the corresponding feature vector in X0:
%
%     Y(i,j,1) = (SUM_d (X(i,j,d) - X0(i,j,d))^P)^(1/P)
%
%   X0 should have the same size as X; the outoput Y has the same
%   height and width as X, but depth equal to 1. Optionally, X0 can
%   be a 1 x 1 x D x N array, in which case the same target feature
%   vector in X0 is compared to all feature vectors in X. In that case,
%   however, the DZDX0 are of size of X.
%
%   Setting the `noRoot` option to `true` does not take the 1/P power
%   in the formula, computing instead
%
%     Y(i,j,1) = SUM_d (X(i,j,d) - X0(i,j,d))^P
%
%   For example, `vl_nnpdist(x, x0, 2, 'noRoot', true)` computes the
%   squared L2 distance.
%
%   [DZDX, DZDX0] = VL_NNPDISTP(X, X0, P, DZDY) computes the derivative
%   of the block inputs projected onto DZDY. DZDX, DZDX0 and DZDY have the
%   same dimensions as X and Y, respectively.
%
%   VL_NNPDIST(___, 'OPT', VAL, ...) accepts the following options:
%
%   `NoRoot`:: `false`
%      If set to true, compute the P-distance to the P-th power.
%
%   `Epsilon`:: 1e-6
%      When computing derivatives, quantities that are divided in are
%      lower boudned by this value. For example, the L2 distance is
%      not smooth at the origin; this option prevents the
%      derivative from diverging.
%
%   `Aggregate`:: false
%      Instead of returning one scalar for each spatial location in
%      the inputs, sum all of them into a single scalar.
%
%   `InstanceWeights``:: `[]`
%      Optionally weight individual instances. This parameter can be
%      eigther a scalar or a weight mask, one for each pixel in the
%      input tensor.

% Copyright (C) 2015  Karel Lenc and Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

% -------------------------------------------------------------------------
%                                                             Parse options
% -------------------------------------------------------------------------

opts.noRoot = false ;
opts.epsilon = 1e-6 ;
opts.aggregate = false ;
opts.instanceWeights = [] ;
backMode = numel(varargin) > 0 && ~ischar(varargin{1}) ;
if backMode
  dzdy = varargin{1} ;
  opts = vl_argparse(opts, varargin(2:end), 'nonrecursive') ;
else
  dzdy = [] ;
  opts = vl_argparse(opts, varargin, 'nonrecursive') ;
end

% -------------------------------------------------------------------------
%                                                             Parse options
% -------------------------------------------------------------------------

d = bsxfun(@minus, x, x0) ;

if ~isempty(dzdy) && ~isempty(opts.instanceWeights)
  dzdy = bsxfun(@times, opts.instanceWeights, dzdy) ;
end

if ~opts.noRoot
  if isempty(dzdy)
    if p == 1
      y1 = sum(abs(d),3) ;
    elseif p == 2
      y1 = sqrt(sum(d.*d,3)) ;
    else
      y1 = sum(abs(d).^p,3).^(1/p) ;
    end
  else
    if p == 1
      y1 = bsxfun(@times, dzdy, sign(d)) ;
    elseif p == 2
      y1 = max(sum(d.*d,3), opts.epsilon).^(-0.5) ;
      y1 = bsxfun(@times, bsxfun(@times, dzdy, y1),  d) ;
    elseif p < 1
      y1 = sum(abs(d).^p,3).^((1-p)/p) ;
      y1 = bsxfun(@times, bsxfun(@times, dzdy, y1), max(abs(d), opts.epsilon).^(p-1) .* sign(d)) ;
    else
      y1 = max(sum(abs(d).^p,3), opts.epsilon).^((1-p)/p) ;
      y1 = bsxfun(@times, bsxfun(@times, dzdy, y1), abs(d).^(p-1) .* sign(d)) ;
    end
  end
else
  if isempty(dzdy)
    if p == 1
      y1 = sum(abs(d),3) ;
    elseif p == 2
      y1 = sum(d.*d,3) ;
    else
      y1 = sum(abs(d).^p,3) ;
    end
  else
    if p == 1
      y1 = bsxfun(@times, dzdy, sign(d)) ;
    elseif p == 2
      y1 = bsxfun(@times, 2 * dzdy, d) ;
    elseif p < 1
      y1 = bsxfun(@times, p * dzdy, max(abs(d), opts.epsilon).^(p-1) .* sign(d)) ;
    else
      y1 = bsxfun(@times, p * dzdy, abs(d).^(p-1) .* sign(d)) ;
    end
  end
end

if isempty(dzdy)
  if ~isempty(opts.instanceWeights)
    y1 = bsxfun(@times, opts.instanceWeights, y1) ;
  end
  if opts.aggregate
    y1 = sum(sum(y1)) ;
  end
end
if ~isempty(dzdy), y2 = -y1; end

repo: CREST-Release


text:
function gX = biasKernGradX(kern, X, X2)

% BIASKERNGRADX Gradient of BIAS kernel with respect to input locations.
% FORMAT
% DESC computes the gradident of the bias
% kernel with respect to the input positions where both the row
% positions and column positions are provided separately.
% ARG kern : kernel structure for which gradients are being
% computed.
% ARG x1 : row locations against which gradients are being computed.
% ARG x2 : column locations against which gradients are being computed.
% RETURN g : the returned gradients. The gradients are returned in
% a matrix which is numData2 x numInputs x numData1. Where numData1 is
% the number of data points in X1, numData2 is the number of data
% points in X2 and numInputs is the number of input
% dimensions in X.
%
% SEEALSO biasKernParamInit, kernGradX, biasKernDiagGradX
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006

% KERN

gX = zeros(size(X2, 1), size(X2, 2), size(X, 1));


text:
% von Bertalanffy 1.1 : a COBRA toolbox extension to thermodynamically constrain a metabolic model
% Ronan M. T. Fleming
% Science Institute and Center for Systems Biology, University of Iceland, Reykjavik.
% Ines Thiele
% Center for Systems Biology and Faculty of Industrial Engineering, Mechanical Engineering and
% Computer Science, University of Iceland, Reykjavik
%
%INSTALLATION
initVonBertalanffy
%TESTING
testSetupThermoModel
%PDF DOCUMENTATION
%
%Theory:
%R.M.T. Fleming†, I. Thiele, H.P. Nasheuer. Quantitative assignment of reaction directionality in
%constraint-based models of metabolism: Application to Escherichia coli.
%Biophysical Chemistry, 145, 47-56, 2009.
%
%Code:
%von Bertalanffy 1.0 : a COBRA toolbox extension to thermodynamically constrain a metabolic model
%(submitted)
%
%HTML DOCUMENTATION
%open YOURPATH/vonBertalanffy/doc/index.html in browser.
%
%html documentation can be renerated using m2html by running
%initVonBertalanffy.m if no doc folder is present
%http://www.artefact.tk/software/matlab/m2html/

author: opencobra

repo: cobratoolbox


text:
function scl = scale( disc, uFun )
% SCALE   Estimate the vertical scale from a vector of values.
%
%   SCL = SCALE( DISC, UFUN ) returns an estimate for the vertical scale
%   of UFUN. VALSDISCRETIZATION expects UFUN to be a cell array of sample data.

scl = max( cellfun(@(v) max(abs(v)), uFun) );

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function out=ctranspose(x)

out=conj(transpose(x));

%%%x=conj(x);
%%%precision=x(1).precision;
%%%out_rval=cell(fliplr(size(x)));
%%%out_ival=cell(fliplr(size(x)));
%%%
%%%xs=size(x)
%%%foo=1;
%%%for jj=1:xs(2)
%%% for ii=1:xs(1)
%%%  [xrval,xrexp,xival,xiexp]=getVals(x,foo);
%%%  out_rval{jj,ii}=xrval;
%%%  out_rexp{jj,ii}=xrexp;
%%%  out_ival{jj,ii}=xival;
%%%  out_iexp{jj,ii}=xiexp;
%%%  foo=foo+1;
%%% end
%%%end
%%%out=class(struct('rval',out_rval,...
%%%                 'rexp',out_rexp,...
%%%                 'ival',out_ival,...
%%%                 'iexp',out_iexp,...
%%%                 'precision',precision),'mp');


text:
function pols = llegepols1 ( degree, x )

%*****************************************************************************80
%
%% LLEGEPOLS1 evaluates orthogonal polynomials on the symmetric interval.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 July 2014
%
%  Author:
%
%    Original FORTRAN77 version by Hong Xiao, Zydrunas Gimbutas.
%    This MATLAB version by John Burkardt.
%
%  Reference:
%
%    Hong Xiao, Zydrunas Gimbutas,
%    A numerical algorithm for the construction of efficient quadrature
%    rules in two and higher dimensions,
%    Computers and Mathematics with Applications,
%    Volume 59, 2010, pages 663-676.
%
%  Parameters:
%
%    Input, integer DEGREE, the maximum degree.
%
%    Input, real X, the evaluation point.
%
%    Output, real POLS(DEGREE+1)), the orthogonal
%    polynomials evaluated at X.
%
  pols = zeros ( degree + 1, 1 );

  pkp1 = 1.0;
  pols(1) = pkp1;

  if ( degree == 0 )
    return
  end

  pk = pkp1;
  pkp1 = x;
  pols(2) = pkp1;

  if ( degree == 1 )
    return
  end

  for k = 1 : degree - 1

    pkm1 = pk;
    pk = pkp1;
    pkp1 = ( ( 2 * k + 1 ) * x * pk ...
           - (     k     ) * pkm1 ) ...
           / (     k + 1 );

    pols(k+2) = pkp1;

  end

  return
end


text:
function [cfg] = ft_neighbourplot(cfg, data)

% FT_NEIGHBOURPLOT visualizes neighbouring channels in a particular channel
% configuration. The positions of the channel are specified in a
% gradiometer or electrode configuration or from a layout.
%
% Use as
%   ft_neighbourplot(cfg)
% or as
%   ft_neighbourplot(cfg, data)
%
% Where the configuration can contain
%   cfg.verbose       = string, 'yes' or 'no', whether the function will print feedback text in the command window
%   cfg.neighbours    = neighbourhood structure, see FT_PREPARE_NEIGHBOURS (optional)
%   cfg.visible       = string, 'on' or 'off', whether figure will be visible (default = 'on')
%   cfg.enableedit    = string, 'yes' or 'no', allows you to interactively add or remove edges between vertices (default = 'no')
%
% and either one of the following options
%   cfg.layout        = filename of the layout, see FT_PREPARE_LAYOUT
%   cfg.elec          = structure with electrode positions or filename, see FT_READ_SENS
%   cfg.grad          = structure with gradiometer definition or filename, see FT_READ_SENS
%   cfg.opto          = structure with gradiometer definition or filename, see FT_READ_SENS
%
% If cfg.neighbours is not defined, this function will call
% FT_PREPARE_NEIGHBOURS to determine the channel neighbours. The
% following data fields may also be used by FT_PREPARE_NEIGHBOURS
%   data.elec         = structure with electrode positions
%   data.grad         = structure with gradiometer definition
%   data.opto         = structure with optode definition
%
% If cfg.neighbours is empty, no neighbouring sensors are assumed.
%
% Use cfg.enableedit to interactively add or remove edges in your own neighbour structure.
%
% See also FT_PREPARE_NEIGHBOURS, FT_PREPARE_LAYOUT

% Copyright (C) 2011, J?rn M. Horschig, 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/>.

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar data
ft_preamble provenance data
ft_preamble trackconfig

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% the data can be passed as input arguments or can be read from disk
hasdata = exist('data', 'var');

if hasdata
  % check if the input data is valid for this function
  data = ft_checkdata(data);
end

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'renamed', {'elecfile', 'elec'});
cfg = ft_checkconfig(cfg, 'renamed', {'gradfile', 'grad'});
cfg = ft_checkconfig(cfg, 'renamed', {'optofile', 'opto'});

% set the defaults
cfg.enableedit = ft_getopt(cfg, 'enableedit', 'no');
cfg.visible    = ft_getopt(cfg, 'visible', 'on');
cfg.renderer   = ft_getopt(cfg, 'renderer'); % let MATLAB decide on the default

if isfield(cfg, 'neighbours')
  cfg.neighbours = cfg.neighbours;
elseif hasdata
  cfg.neighbours = ft_prepare_neighbours(cfg, data);
else
  cfg.neighbours = ft_prepare_neighbours(cfg);
end

if ~isfield(cfg, 'verbose')
  cfg.verbose = 'no';
elseif strcmp(cfg.verbose, 'yes')
  cfg.verbose = true;
end

% get the the grad or elec
if hasdata
  sens = ft_fetch_sens(cfg, data);
else
  sens = ft_fetch_sens(cfg);
end

% insert sensors that are not in neighbourhood structure
if isempty(cfg.neighbours)
  nsel = 1:numel(sens.label);
else
  nsel = find(~ismember(sens.label, {cfg.neighbours.label}));
end

for i=1:numel(nsel)
  cfg.neighbours(end+1).label = sens.label{nsel(i)};
  cfg.neighbours(end).neighblabel = {};
end

[tmp, sel] = match_str(sens.label, {cfg.neighbours.label});
cfg.neighbours = cfg.neighbours(sel);

% give some graphical feedback
if all(sens.chanpos(:,3)==0)
  % the sensor positions are already projected on a 2D plane
  proj = sens.chanpos(:,1:2);
else
  % use 3-dimensional data for plotting
  proj = sens.chanpos;
end
hf = figure('visible', cfg.visible);
axis equal
axis vis3d
axis off
hold on;
hl = [];
for i=1:length(cfg.neighbours)
  this = cfg.neighbours(i);

  sel1 = match_str(sens.label, this.label);
  sel2 = match_str(sens.label, this.neighblabel);
  % account for missing sensors
  this.neighblabel = sens.label(sel2);
  for j=1:length(this.neighblabel)
    x1 = proj(sel1,1);
    y1 = proj(sel1,2);
    x2 = proj(sel2(j),1);
    y2 = proj(sel2(j),2);
    X = [x1 x2];
    Y = [y1 y2];
    if size(proj, 2) == 2
      hl(sel1, sel2(j)) = line(X, Y, 'color', 'r');
    elseif size(proj, 2) == 3
      z1 = proj(sel1,3);
      z2 = proj(sel2(j),3);
      Z = [z1 z2];
      hl(sel1, sel2(j)) = line(X, Y, Z, 'color', 'r');
    end
  end
end

% this is for putting the channels on top of the connections
hs = [];
for i=1:length(cfg.neighbours)
  this = cfg.neighbours(i);
  sel1 = match_str(sens.label, this.label);
  sel2 = match_str(sens.label, this.neighblabel);
  % account for missing sensors
  this.neighblabel = sens.label(sel2);
  if isempty(sel1)
    continue;
  end
  if size(proj, 2) == 2
    hs(i) = line(proj(sel1, 1), proj(sel1, 2),                                            ...
      'MarkerEdgeColor',  'k',                                        ...
      'MarkerFaceColor',  'k',                                        ...
      'Marker',           'o',                                        ...
      'MarkerSize',       .125*(2+numel(cfg.neighbours(i).neighblabel))^2, ...
      'UserData',         i,                                          ...
      'ButtonDownFcn',    @showLabelInTitle);

  elseif size(proj, 2) == 3
    hs(i) = line(proj(sel1, 1), proj(sel1, 2), proj(sel1, 3),                                ...
      'MarkerEdgeColor',  'k',                                        ...
      'MarkerFaceColor',  'k',                                        ...
      'Marker',           'o',                                        ...
      'MarkerSize',       .125*(2+numel(cfg.neighbours(i).neighblabel))^2, ...
      'UserData',         i,                        ...
      'ButtonDownFcn',    @showLabelInTitle);
  else
    ft_error('Channel coordinates are too high dimensional');
  end
end

hold off
title('[Click on a sensor to see its label]');

% store what is needed in UserData of figure
userdata.lastSensId = [];
userdata.cfg = cfg;
userdata.sens = sens;
userdata.hs = hs;
userdata.hl = hl;
userdata.quit = false;
hf = getparent(hf);
set(hf, 'UserData', userdata);

if istrue(cfg.enableedit)
  set(hf, 'CloseRequestFcn', @cleanup_cb);
  while ~userdata.quit
    uiwait(hf);
    userdata = get(hf, 'UserData');
  end
  cfg = userdata.cfg;

  hf = getparent(hf);
  delete(hf);
end

% remove SCALE and COMNT
desired = ft_channelselection({'all', '-SCALE', '-COMNT'}, {cfg.neighbours.label});

neighb_idx = ismember({cfg.neighbours.label}, desired);
cfg.neighbours = cfg.neighbours(neighb_idx);

% this is needed for the figure title
if isfield(cfg, 'dataname') && ~isempty(cfg.dataname)
  dataname = cfg.dataname;
elseif isfield(cfg, 'inputfile') && ~isempty(cfg.inputfile)
  dataname = cfg.inputfile;
elseif nargin>1
  dataname = arrayfun(@inputname, 2:nargin, 'UniformOutput', false);
else
  dataname = {};
end

% set the figure window title
if ~isempty(dataname)
  set(gcf, 'Name', sprintf('%d: %s: %s', double(gcf), mfilename, join_str(', ', dataname)));
else
  set(gcf, 'Name', sprintf('%d: %s', double(gcf), mfilename));
end
set(gcf, 'NumberTitle', 'off');

% set renderer if specified
if ~isempty(cfg.renderer)
  set(gcf, 'renderer', cfg.renderer)
end

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble trackconfig
ft_postamble previous data
ft_postamble provenance
ft_postamble savefig

% add a menu to the figure, but only if the current figure does not have subplots
menu_fieldtrip(gcf, cfg, false);

if ~ft_nargout
  % don't return anything
  clear cfg
end

end % main function

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function showLabelInTitle(gcbo, EventData, handles)

userdata    = get(gcf, 'UserData');
lastSensId  = userdata.lastSensId;
cfg         = userdata.cfg;
hs          = userdata.hs;
curSensId   = get(gcbo, 'UserData');

if lastSensId == curSensId

  title('[Click on a sensor to see its label]');
  set(hs(curSensId), 'MarkerFaceColor', 'k');
  userdata.lastSensId = [];

elseif isempty(lastSensId) || ~istrue(cfg.enableedit)

  userdata.lastSensId = curSensId;
  if istrue(cfg.enableedit)
    title(['Selected channel: ' cfg.neighbours(curSensId).label ' click on another to (dis-)connect']);
  else
    title(['Selected channel: ' cfg.neighbours(curSensId).label]);
  end
  if cfg.verbose
    str = sprintf('%s, ', cfg.neighbours(curSensId).neighblabel{:});
    if length(str)>2
      % remove the last comma and space
      str = str(1:end-2);
    end
    fprintf('Selected channel %s, which has %d neighbours: %s\n', ...
      cfg.neighbours(curSensId).label, ...
      length(cfg.neighbours(curSensId).neighblabel), ...
      str);
  end
  set(hs(curSensId), 'MarkerFaceColor', 'g');
  set(hs(lastSensId), 'MarkerFaceColor', 'k');

elseif istrue(cfg.enableedit)
  hl    = userdata.hl;
  sens  = userdata.sens;
  if all(sens.chanpos(:,3)==0)
    % the sensor positions are already projected on a 2D plane
    proj = sens.chanpos(:,1:2);
  else
    % use 3-dimensional data for plotting
    proj = sens.chanpos;
  end

  % find out whether they are connected
  connected1 = ismember(cfg.neighbours(curSensId).neighblabel, cfg.neighbours(lastSensId).label);
  connected2 = ismember(cfg.neighbours(lastSensId).neighblabel, cfg.neighbours(curSensId).label);

  if any(connected1) % then disconnect
    cfg.neighbours(curSensId).neighblabel(connected1) = [];
    cfg.neighbours(lastSensId).neighblabel(connected2) = [];
    title(['Disconnected channels ' cfg.neighbours(curSensId).label ' and ' cfg.neighbours(lastSensId).label]);
    delete(hl(curSensId, lastSensId));
    hl(curSensId, lastSensId) = 0;
    delete(hl(lastSensId, curSensId));
    hl(lastSensId, curSensId) = 0;
  else % then connect
    cfg.neighbours(curSensId).neighblabel{end+1} = cfg.neighbours(lastSensId).label;
    cfg.neighbours(lastSensId).neighblabel{end+1} = cfg.neighbours(curSensId).label;
    title(['Connected channels ' cfg.neighbours(curSensId).label ' and ' cfg.neighbours(lastSensId).label]);

    % draw new edge
    x1 = proj(curSensId,1);
    y1 = proj(curSensId,2);
    x2 = proj(lastSensId,1);
    y2 = proj(lastSensId,2);
    X = [x1 x2];
    Y = [y1 y2];
    if size(proj, 2) == 2
      hl(curSensId, lastSensId) = line(X, Y, 'color', 'r');
      hl(lastSensId, curSensId) = line(X, Y, 'color', 'r');
    elseif size(proj, 2) == 3
      z1 = proj(curSensId,3);
      z2 = proj(lastSensId,3);
      Z =[z1 z2];
      hl(curSensId, lastSensId) = line(X, Y, Z, 'color', 'r');
      hl(lastSensId, curSensId) = line(X, Y, Z, 'color', 'r');
    end

  end
  % draw nodes on top again
  delete(hs(curSensId));
  delete(hs(lastSensId));
  if size(proj, 2) == 2
    hs(curSensId) = line(proj(curSensId, 1), proj(curSensId, 2),                                            ...
      'MarkerEdgeColor',  'k',                                        ...
      'MarkerFaceColor',  'k',                                        ...
      'Marker',           'o',                                        ...
      'MarkerSize',       .125*(2+numel(cfg.neighbours(curSensId).neighblabel))^2, ...
      'UserData',         curSensId,                                          ...
      'ButtonDownFcn',    @showLabelInTitle);
    hs(lastSensId) = line(proj(lastSensId, 1), proj(lastSensId, 2),                                            ...
      'MarkerEdgeColor',  'k',                                        ...
      'MarkerFaceColor',  'k',                                        ...
      'Marker',           'o',                                        ...
      'MarkerSize',       .125*(2+numel(cfg.neighbours(lastSensId).neighblabel))^2, ...
      'UserData',         lastSensId,                                          ...
      'ButtonDownFcn',    @showLabelInTitle);

  elseif size(proj, 2) == 3
    hs(curSensId) = line(proj(curSensId, 1), proj(curSensId, 2), proj(curSensId, 3),                                ...
      'MarkerEdgeColor',  'k',                                        ...
      'MarkerFaceColor',  'k',                                        ...
      'Marker',           'o',                                        ...
      'MarkerSize',       .125*(2+numel(cfg.neighbours(curSensId).neighblabel))^2, ...
      'UserData',         curSensId,                        ...
      'ButtonDownFcn',    @showLabelInTitle);
    hs(lastSensId) = line(proj(lastSensId, 1), proj(lastSensId, 2), proj(lastSensId, 3),                                ...
      'MarkerEdgeColor',  'k',                                        ...
      'MarkerFaceColor',  'k',                                        ...
      'Marker',           'o',                                        ...
      'MarkerSize',       .125*(2+numel(cfg.neighbours(lastSensId).neighblabel))^2, ...
      'UserData',         lastSensId,                        ...
      'ButtonDownFcn',    @showLabelInTitle);
  else
    ft_error('Channel coordinates are too high dimensional');
  end

  if cfg.verbose
    str = sprintf('%s, ', cfg.neighbours(curSensId).neighblabel{:});
    if length(str)>2
      % remove the last comma and space
      str = str(1:end-2);
    end
    fprintf('Selected channel %s, which has %d neighbours: %s\n', ...
      cfg.neighbours(curSensId).label, ...
      length(cfg.neighbours(curSensId).neighblabel), ...
      str);
  end
  set(hs(curSensId), 'MarkerFaceColor', 'g');
  set(hs(lastSensId), 'MarkerFaceColor', 'k');
  userdata.lastSensId = curSensId;
  userdata.hl = hl;
  userdata.hs = hs;
  userdata.cfg = cfg;
  set(gcf, 'UserData', userdata);
  return;
else
  % can never happen, so do nothing
end

set(gcf, 'UserData', userdata);
end % subfunction

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cleanup_cb(h, eventdata)
userdata = get(h, 'UserData');
h   = getparent(h);
userdata.quit = true;
set(h, 'UserData', userdata);
uiresume
end % subfunction

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function h = getparent(h)
p = h;
while p~=0
  h = p;
  p = get(h, 'parent');
end
end % subfunction


text:
%% The Index of EBSD data
%
%%
% In previous chapters we have discussed how to select EBSD data by
% properties. In this chapter we are intested in the order the EBSD are
% stored within MTEX. Lets start by importing some sample data

mtexdata twins

%%
% and restricting it to very small rectangular subset

poly = [44 0 4 2];
ebsd = ebsd(inpolygon(ebsd,poly));

plot(ebsd,ebsd.orientations,'micronbar','off','edgecolor','k')

%%
% In the above plot each square corresponds to one entry in the variable
% |ebsd|. Lets visualize the order

text(ebsd,1:length(ebsd))

%%
% We may easily select specific measurement pixels by specifying their
% indeces

hold on
plot(ebsd(16:18),'edgeColor','red','facecolor','none','linewidth',4)
legend off
hold off

%%
% Whether lines or columns run first is not related to MTEX but inherits
% from the ordering of the imported EBSD data. Since, we have restricted
% our large EBSD map to the small subset the indece of restricted data does
% not coincide with the indece of the imported data anymore. However, the
% original indeces are still stored in |ebsd.id|. Lets visualize them

plot(ebsd,ebsd.orientations,'micronbar','off','edgecolor','k')
text(ebsd,ebsd.id)

%%
% In order to select EBSD data according to their original id use the
% option |'id'|, i.e.,

hold on
plot(ebsd('id',316:318),'edgeColor','red','facecolor','none','linewidth',4)
legend off
hold off

%% Square Grids
%
% In the cases of gridded data it is often useful to convert them into a
% matrix form.

ebsd = ebsd.gridify;

plot(ebsd,ebsd.orientations,'micronbar','off')

[i,j] = ndgrid(1:size(ebsd,1),1:size(ebsd,2));
str = arrayfun(@(a,b) ['(' int2str(a) ',' int2str(b) ')'],i,j,'UniformOutput',false);
text(ebsd,str)

%%
% This allows to select EBSD data simply by their coordinates within the
% grid, e.g., by

hold on
plot(ebsd(2,2:4),'edgeColor','red','facecolor','none','linewidth',4)
legend off
hold off

%%
% Note that the <EBSD.gridify.html |gridify|> command changes the order of
% measurements. They are now sort such that rows runs first and columns
% second, as this is the default convention how Matlab indexes matrices.

plot(ebsd,ebsd.orientations,'micronbar','off')
text(ebsd,1:length(ebsd))

%% Hexagonal Grid
%
% The command <EBSD.gridify.html |gridify|> may also be applied to EBSD
% data measured on a hexagonal grid.

mtexdata titanium silent

ebsd = ebsd.gridify

%%
% This rearranges the measurements in a matrix form which can be indexed
% similarly as in the square case.

ebsd = ebsd(10:16,68:79);

%%
% Lets visualize the matrix coordinates for the hexagonal grid

plot(ebsd,ebsd.orientations,'edgeColor','k','micronbar','off')
axis off

[i,j] = ndgrid(1:size(ebsd,1),1:size(ebsd,2));
str = arrayfun(@(a,b) ['(' int2str(a) ',' int2str(b) ')'],i,j,'UniformOutput',false);
text(ebsd,str)

%% Cube Coordinates
% In hexognal grids it is sometimes advantageous to use three digit cube
% coordinates to index the cell. This can be done using the commands
% <EBSDhex.hex2cube.html |hex2cube|> and <EBSDhex.cube2hex.html
% |cube2hex|>. Much more details on indexing hex grids can be found at
% <https://www.redblobgames.com/grids/hexagons/ here>.

plot(ebsd,ebsd.orientations,'edgeColor','k','micronbar','off')
axis off

[i,j] = ndgrid(1:size(ebsd,1),1:size(ebsd,2));
[x,y,z] = ebsd.hex2cube(i,j);
str = arrayfun(@(a,b,c) ['(' int2str(a) ',' int2str(b) ',' int2str(c) ')'],x,y,z,'UniformOutput',false);
text(ebsd,str)

author: mtex-toolbox

repo: mtex

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


text:
function [X_norm, mu, sigma] = featureNormalize(X)
%FEATURENORMALIZE Normalizes the features in X
%   FEATURENORMALIZE(X) returns a normalized version of X where
%   the mean value of each feature is 0 and the standard deviation
%   is 1. This is often a good preprocessing step to do when
%   working with learning algorithms.

% You need to set these values correctly
X_norm = X;
numColumns = size(X, 2);

mu = mean(X);
sigma = std(X);
for i = 1:numColumns
    X_norm(:,i) = (X(:, i) - mu(i)) / sigma(i);
end;

% ====================== YOUR CODE HERE ======================
% Instructions: First, for each feature dimension, compute the mean
%               of the feature and subtract it from the dataset,
%               storing the mean value in mu. Next, compute the
%               standard deviation of each feature and divide
%               each feature by it's standard deviation, storing
%               the standard deviation in sigma.
%
%               Note that X is a matrix where each column is a
%               feature and each row is an example. You need
%               to perform the normalization separately for
%               each feature.
%
% Hint: You might find the 'mean' and 'std' functions useful.
%

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

end

author: zsiciarz

repo: ml-coursera


text:
function sendJointsPositionsf( t_Kuka ,jPos)
%% Syntax:
% sendJointsPositionsf( t_Kuka ,jPos)

%% About
% This function is used to send the target positions of the joints for the
% direct servo motion
% this function is for direct-servoing in joint space

%% Arreguments
% jPos: is 7 cells array of doubles
% t: is the TCP/IP connection object

% Copy right, Mohammad SAFEEA, 3rd of May 2017

theCommand='jf_';
for i=1:7
    x=jPos{i};
    theCommand=[theCommand,sprintf('%0.5f',x),'_'];
end
fprintf(t_Kuka, theCommand);
end

repo: KST-Kuka-Sunrise-Toolbox


text:
function pass = test_uminus(pref)
% Test chebfun3/uminus.

if ( nargin < 1 )
    pref = chebfunpref;
end
tol = 1e5 * pref.cheb3Prefs.chebfun3eps;

dom = [-1 1 -1 1 -1 1;
       -2 2 -2 2 -2 2;
       -1 pi 0 2*pi -pi pi];

for j = 1 : size(dom,1)
    f = chebfun3(@(x,y,z) cos(x.*y.*z), dom(j,:));

    uminusF = chebfun3(@(x,y,z) -cos(x.*y.*z), dom(j,:));

    tolk = norm(dom(j, :), inf) * tol;

    pass(j) = norm((-f) - uminusF) < tolk;
end

end


text:
function y = getNumericScalarFieldSafe(data,fld,mode)
% y = getNumericScalarFieldSafe(data,fld,mode)
% mode: optional, either 'safe' or 'safest' (default).
% y: col vector

if nargin < 3
    mode = 'safest';
end

if isempty(data)
    y = zeros(size(data));
else
    switch mode
        case 'safest'
            y = {data.(fld)}';
            cellfun(@(x)assert((isnumeric(x)||islogical(x))&&isscalar(x)),y);
            y = cell2mat(y);
            assert(isequal(size(y),[numel(data) 1]));
        case 'safe'
            y = cat(1,data.(fld));
            assert(isequal(size(y),[numel(data) 1]));
        otherwise
            assert(false);
    end
end

author: kristinbranson

repo: JAABA


text:
function gt=gabtight(varargin)
%GABTIGHT  Canonical tight window of Gabor frame
%   Usage:  gt=gabtight(a,M,L);
%           gt=gabtight(g,a,M);
%           gt=gabtight(g,a,M,L);
%           gd=gabtight(g,a,M,'lt',lt);
%
%   Input parameters:
%         g     : Gabor window.
%         a     : Length of time shift.
%         M     : Number of modulations.
%         L     : Length of window. (optional)
%         lt    : Lattice type (for non-separable lattices).
%   Output parameters:
%         gt    : Canonical tight window, column vector.
%
%   `gabtight(a,M,L)` computes a nice tight window of length *L* for a
%   lattice with parameters *a*, *M*. The window is not an FIR window,
%   meaning that it will only generate a tight system if the system
%   length is equal to *L*.
%
%   `gabtight(g,a,M)` computes the canonical tight window of the Gabor frame
%   with window *g* and parameters *a*, *M*.
%
%   The window *g* may be a vector of numerical values, a text string or a
%   cell array. See the help of |gabwin| for more details.
%
%   If the length of *g* is equal to *M*, then the input window is assumed to
%   be a FIR window. In this case, the canonical dual window also has
%   length of *M*. Otherwise the smallest possible transform length is
%   chosen as the window length.
%
%   `gabtight(g,a,M,L)` returns a window that is tight for a system of
%   length *L*. Unless the input window *g* is a FIR window, the returned
%   tight window will have length *L*.
%
%   `gabtight(g,a,M,'lt',lt)` does the same for a non-separable lattice
%   specified by *lt*. Please see the help of |matrix2latticetype| for a
%   precise description of the parameter *lt*.
%
%   If $a>M$ then an orthonormal window of the Gabor Riesz sequence with
%   window *g* and parameters *a* and *M* will be calculated.
%
%   Examples:
%   ---------
%
%   The following example shows the canonical tight window of the Gaussian
%   window. This is calculated by default by |gabtight| if no window is
%   specified:::
%
%     a=20;
%     M=30;
%     L=300;
%     gt=gabtight(a,M,L);
%
%     % Simple plot in the time-domain
%     figure(1);
%     plot(gt);
%
%     % Frequency domain
%     figure(2);
%     magresp(gt,'dynrange',100);
%
%   See also:  gabdual, gabwin, fir2long, dgt

%   AUTHOR : Peter L. Søndergaard.
%   TESTING: TEST_DGT
%   REFERENCE: OK

%% ------------ decode input parameters ------------

if nargin<3
    error('%s: Too few input parameters.',upper(mfilename));
end;

if numel(varargin{1})==1
  % First argument is a scalar.

  a=varargin{1};
  M=varargin{2};

  g='gauss';

  varargin=varargin(3:end);
else
  % First argument assumed to be a vector.

  g=varargin{1};
  a=varargin{2};
  M=varargin{3};

  varargin=varargin(4:end);
end;

definput.keyvals.L=[];
definput.keyvals.lt=[0 1];
definput.keyvals.nsalg=0;
[flags,kv,L]=ltfatarghelper({'L'},definput,varargin);

%% ------ step 2: Verify a, M and L
if isempty(L)
    if isnumeric(g)
        % Use the window length
        Ls=length(g);
    else
        % Use the smallest possible length
        Ls=1;
    end;

    % ----- step 2b : Verify a, M and get L from the window length ----------
    L=dgtlength(Ls,a,M,kv.lt);

else

    % ----- step 2a : Verify a, M and get L

    Luser=dgtlength(L,a,M,kv.lt);
    if Luser~=L
        error(['%s: Incorrect transform length L=%i specified. Next valid length ' ...
               'is L=%i. See the help of DGTLENGTH for the requirements.'],...
              upper(mfilename),L,Luser);
    end;

end;

%% ----- step 3 : Determine the window

[g,info]=gabwin(g,a,M,L,kv.lt,'callfun',upper(mfilename));

if L<info.gl
  error('%s: Window is too long.',upper(mfilename));
end;

R=size(g,2);

% -------- Are we in the Riesz sequence of in the frame case

scale=1;
if a>M*R
  % Handle the Riesz basis (dual lattice) case.
  % Swap a and M, and scale differently.
  scale=sqrt(a/M);
  tmp=a;
  a=M;
  M=tmp;
end;

% -------- Compute -------------

if kv.lt(2)==1
    % Rectangular case

    if (info.gl<=M) && (R==1)

        % Diagonal of the frame operator
        d = gabframediag(g,a,M,L);
        gt=g./sqrt(long2fir(d,info.gl));

    else

        % Long window case

        % Just in case, otherwise the call is harmless.
        g=fir2long(g,L);

        gt=comp_gabtight_long(g,a,M)*scale;

    end;

else

    % Just in case, otherwise the call is harmless.
    g=fir2long(g,L);

    if (kv.nsalg==1) || (kv.nsalg==0 && kv.lt(2)<=2)

        mwin=comp_nonsepwin2multi(g,a,M,kv.lt,L);

        gtfull=comp_gabtight_long(mwin,a*kv.lt(2),M)*scale;

        % We need just the first vector
        gt=gtfull(:,1);

    else

        [s0,s1,br] = shearfind(L,a,M,kv.lt);

        if s1 ~= 0
            p1 = comp_pchirp(L,s1);
            g = p1.*g;
        end

        b=L/M;
        Mr = L/br;
        ar = a*b/br;

        if s0 == 0
            gt=comp_gabtight_long(g,ar,Mr);
        else
            p0=comp_pchirp(L,-s0);
            g = p0.*fft(g);
            gt=comp_gabtight_long(g,L/Mr,L/ar)*sqrt(L);
            gt = ifft(conj(p0).*gt);
        end

        if s1 ~= 0
            gt = conj(p1).*gt;
        end

    end;

    if (info.gl<=M) && (R==1)
        gt=long2fir(gt,M);
    end;

end;

% --------- post process result -------

if isreal(g) && (kv.lt(2)<=2)
  % If g is real and the lattice is either rectangular or quinqux, then
  % the output is known to be real.
  gt=real(gt);
end;

if info.wasrow
  gt=gt.';
end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [ps,ix] = dpsimplify(p,tol)

% Recursive Douglas-Peucker Polyline Simplification, Simplify
%
% [ps,ix] = dpsimplify(p,tol)
%
% dpsimplify uses the recursive Douglas-Peucker line simplification
% algorithm to reduce the number of vertices in a piecewise linear curve
% according to a specified tolerance. The algorithm is also know as
% Iterative Endpoint Fit. It works also for polylines and polygons
% in higher dimensions.
%
% In case of nans (missing vertex coordinates) dpsimplify assumes that
% nans separate polylines. As such, dpsimplify treats each line
% separately.
%
% For additional information on the algorithm follow this link
% http://en.wikipedia.org/wiki/Ramer-Douglas-Peucker_algorithm
%
% Input arguments
%
%     p     polyline n*d matrix with n vertices in d
%           dimensions.
%     tol   tolerance (maximal euclidean distance allowed
%           between the new line and a vertex)
%
% Output arguments
%
%     ps    simplified line
%     ix    linear index of the vertices retained in p (ps = p(ix))
%
% Examples
%
% 1. Simplify line
%
%     tol    = 1;
%     x      = 1:0.1:8*pi;
%     y      = sin(x) + randn(size(x))*0.1;
%     p      = [x' y'];
%     ps     = dpsimplify(p,tol);
%
%     plot(p(:,1),p(:,2),'k')
%     hold on
%     plot(ps(:,1),ps(:,2),'r','LineWidth',2);
%     legend('original polyline','simplified')
%
% 2. Reduce polyline so that only knickpoints remain by
%    choosing a very low tolerance
%
%     p = [(1:10)' [1 2 3 2 4 6 7 8 5 2]'];
%     p2 = dpsimplify(p,eps);
%     plot(p(:,1),p(:,2),'k+--')
%     hold on
%     plot(p2(:,1),p2(:,2),'ro','MarkerSize',10);
%     legend('original line','knickpoints')
%
% 3. Simplify a 3d-curve
%
%     x = sin(1:0.01:20)';
%     y = cos(1:0.01:20)';
%     z = x.*y.*(1:0.01:20)';
%     ps = dpsimplify([x y z],0.1);
%     plot3(x,y,z);
%     hold on
%     plot3(ps(:,1),ps(:,2),ps(:,3),'k*-');
%
%
%
% Author: Wolfgang Schwanghart, 13. July, 2010.
% w.schwanghart[at]unibas.ch

if nargin == 0
    help dpsimplify
    return
end

error(nargchk(2, 2, nargin))

% error checking
if ~isscalar(tol) || tol<0;
    error('tol must be a positive scalar')
end

% nr of dimensions
nrvertices    = size(p,1);
dims    = size(p,2);

% anonymous function for starting point and end point comparision
% using a relative tolerance test
compare = @(a,b) abs(a-b)/max(abs(a),abs(b)) <= eps;

% what happens, when there are NaNs?
% NaNs divide polylines.
Inan      = any(isnan(p),2);
% any NaN at all?
Inanp     = any(Inan);

% if there is only one vertex
if nrvertices == 1 || isempty(p);
    ps = p;
    ix = 1;

% if there are two
elseif nrvertices == 2 && ~Inanp;
    % when the line has no vertices (except end and start point of the
    % line) check if the distance between both is less than the tolerance.
    % If so, return the center.
    if dims == 2;
        d    = hypot(p(1,1)-p(2,1),p(1,2)-p(2,2));
    else
        d    = sqrt(sum((p(1,:)-p(2,:)).^2));
    end

    if d <= tol;
        ps = sum(p,1)/2;
        ix = 1;
    else
        ps = p;
        ix = [1;2];
    end

elseif Inanp;

    % case: there are nans in the p array
    % --> find start and end indices of contiguous non-nan data
    Inan = ~Inan;
    sIX = strfind(Inan',[0 1])' + 1;
    eIX = strfind(Inan',[1 0])';

    if Inan(end)==true;
        eIX = [eIX;nrvertices];
    end

    if Inan(1);
        sIX = [1;sIX];
    end

    % calculate length of non-nan components
    lIX = eIX-sIX+1;
    % put each component into a single cell
    c   = mat2cell(p(Inan,:),lIX,dims);

    % now call dpsimplify again inside cellfun.
    if nargout == 2;
        [ps,ix]   = cellfun(@(x) dpsimplify(x,tol),c,'uniformoutput',false);
        ix        = cellfun(@(x,six) x+six-1,ix,num2cell(sIX),'uniformoutput',false);
    else
        ps   = cellfun(@(x) dpsimplify(x,tol),c,'uniformoutput',false);
    end

    % write the data from a cell array back to a matrix
    ps = cellfun(@(x) [x;nan(1,dims)],ps,'uniformoutput',false);
    ps = cell2mat(ps);
    ps(end,:) = [];

    % ix wanted? write ix to a matrix, too.
    if nargout == 2;
        ix = cell2mat(ix);
    end

else

% if there are no nans than start the recursive algorithm
ixe     = size(p,1);
ixs     = 1;

% logical vector for the vertices to be retained
I   = true(ixe,1);

% call recursive function
p   = simplifyrec(p,tol,ixs,ixe);
ps  = p(I,:);

% if desired return the index of retained vertices
if nargout == 2;
    ix  = find(I);
end

end

% _________________________________________________________
function p  = simplifyrec(p,tol,ixs,ixe)

    % check if startpoint and endpoint are the same
    % better comparison needed which included a tolerance eps

    c1 = num2cell(p(ixs,:));
    c2 = num2cell(p(ixe,:));

    % same start and endpoint with tolerance
    sameSE = all(cell2mat(cellfun(compare,c1(:),c2(:),'UniformOutput',false)));

    if sameSE;
        % calculate the shortest distance of all vertices between ixs and
        % ixe to ixs only
        if dims == 2;
            d    = hypot(p(ixs,1)-p(ixs+1:ixe-1,1),p(ixs,2)-p(ixs+1:ixe-1,2));
        else
            d    = sqrt(sum(bsxfun(@minus,p(ixs,:),p(ixs+1:ixe-1,:)).^2,2));
        end
    else
        % calculate shortest distance of all points to the line from ixs to ixe
        % subtract starting point from other locations
        pt = bsxfun(@minus,p(ixs+1:ixe,:),p(ixs,:));

        % end point
        a = pt(end,:)';

        beta = (a' * pt')./(a'*a);
        b    = pt-bsxfun(@times,beta,a)';
        if dims == 2;
            % if line in 2D use the numerical more robust hypot function
            d    = hypot(b(:,1),b(:,2));
        else
            d    = sqrt(sum(b.^2,2));
        end
    end

    % identify maximum distance and get the linear index of its location
    [dmax,ixc] = max(d);
    ixc  = ixs + ixc;

    % if the maximum distance is smaller than the tolerance remove vertices
    % between ixs and ixe
    if dmax <= tol;
        if ixs ~= ixe-1;
            I(ixs+1:ixe-1) = false;
        end
    % if not, call simplifyrec for the segments between ixs and ixc (ixc
    % and ixe)
    else
        p   = simplifyrec(p,tol,ixs,ixc);
        p   = simplifyrec(p,tol,ixc,ixe);

    end

end
end


text:
function tapas_BLT_runBuffer(params)
% tapas_BLT_runBuffer runs an initial fixation cross

Screen('DrawTexture', params.screen.window, params.text.imgiti, [], params.screen.fit, 0);
Screen('Flip', params.screen.window);
wait2_escapeOption(params.dur.introiti , params);

author: translationalneuromodeling

repo: tapas


text:
function [x, infos] = palm_sparse_smooth_nmf(V, rank, in_options)
% PALM framework with smoothness and sparsity constraints for non-negative
% matrix factorization (PALM-Sparse-Smooth-NMF)
%
% The problem of interest is defined as
%
%           min || WH - V ||_F^2 + lambda * || W ||_1 + eta || HT ||_F^2
%                                   + betaW ||W||_F^2 + betaH ||H||_F^2,
%           where
%           {V, W, H} >= 0,
%
%           Algorithm for NMF with eucidian norm as objective function and
%           L1 constraint on W for sparse paterns and Tikhonov regularization
%           for smooth activation coefficients.
%
% Inputs:
%       matrix      V
%       rank        rank
%       options     options
%           lambda      weight for the L1 sparsity penalty (default: 0)
%           eta         weight for the smoothness constraint.
%           gamma1: 	constant > 1 for the gradient descend step of W.
%           gamma2:     constant > 1 for the gradient descend step of W.
%           betaH:      constant. L-2 constraint for H.
%           betaW:      constant. L-2 constraint for W.
% Output:
%       w           solution of w
%       infos       information
%
% References:
%
%
% This file is part of NMFLibrary.
%
% This file has been ported from
%       palm_nmf.m at https://github.com/raimon-fa/palm-nmf
%
%{
The MIT License
Copyright (c) 2017 Raimon Fabregat
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.
%}
%
% Ported by M.Horie and H.Kasai on June 21, 2022 for NMFLibrary
%
% Change log:
%

    % set dimensions and samples
    [m, n] = size(V);

    % set local options
    local_options = [];
    local_options.lambda = 0;   % sparsity = lambda
    local_options.eta   = 0;    % smoothness = eta
    local_options.betaW = 0.1;
    local_options.betaH = 0.1;
    local_options.gamma1 = 1.001;
    local_options.gamma2 = 1.001;
    local_options.sub_mode = 'std';

    % check input options
    if ~exist('in_options', 'var') || isempty(in_options)
        in_options = struct();
    end
    % merge options
    options = mergeOptions(get_nmf_default_options(), local_options);
    options = mergeOptions(options, in_options);

    % initialize factors
    init_options = options;
    [init_factors, ~] = generate_init_factors(V, rank, init_options);
    W = init_factors.W;
    H = init_factors.H;

    % initialize for PAML-NMF
    TTp_norm = 0;
    TTp = zeros(n);
    if options.lambda == 0 && options.eta == 0
        % NMF
        % In the case without constraints it can be shown that
        % the gammas can be divided by 2 (Bolte 2014)
        options.gamma1 = options.gamma1 / 2;
        options.gamma2 = options.gamma2 / 2;
        options.betaW = 0;
        options.betaH = 0;
    elseif options.lambda > 0 && options.eta == 0
        % sparse NMF
        options.betaW = 0;
        options.sub_mode = 'sparse';
    elseif options.lambda == 0 && options.eta > 0
        % smooth NMF
        % Tikhonov regularization matrix
        T = eye(n) - diag(ones(n-1, 1),-1);
        T = T(:, 1:end-1);
        TTp = T * T';
        TTp_norm = norm(TTp, 'fro');
        options.betaH = 0;
        options.sub_mode = 'smooth';
    elseif options.lambda > 0 && options.eta > 0
        % smooth and sparse NMF
        % Tikhonov regularization matrix
        T = eye(n) - diag(ones(n-1, 1),-1);
        T = T(:, 1:end-1);
        TTp = T * T';
        TTp_norm = norm(TTp, 'fro');
        options.sub_mode = 'smooth & sparse';
    else
        error('Give positive values to the parameters')
    end

    % initialize
    method_name = sprintf('PALM-Sparse-Smooth (%s)', options.sub_mode);
    epoch = 0;
    grad_calc_count = 0;

    if options.verbose > 0
        fprintf('%s: started ...\n', method_name);
    end

    % store initial info
    clear infos;
    [infos, f_val, optgap] = store_nmf_info(V, W, H, [], options, [], epoch, grad_calc_count, 0);

    if options.verbose > 1
        fprintf('%s: Epoch = 0000, cost = %.16e, optgap = %.4e\n', method_name, options.sub_mode, f_val, optgap);
    end

    % set start time
    start_time = tic();

    % main loop
    while true

        % check stop condition
        [stop_flag, reason, max_reached_flag] = check_stop_condition(epoch, infos, options);
        if stop_flag
            display_stop_reason(epoch, infos, options, method_name, reason, max_reached_flag);
            break;
        end

        % update W
        c = options.gamma1 * 2 * (norm(H * H', 'fro') + options.betaW);
        z1 = W - (1 / c) * 2 * ((W * H - V) * H' + options.betaW * W);
        W = max(z1 - 2 * options.lambda / c, 0);

        % update H
        d = options.gamma2 * 2 * (norm(W*W', 'fro') + options.eta * TTp_norm + options.betaH);
        z2 = H - (1 / d) * 2 * (W' * (W * H - V) + options.eta * (H * TTp) + options.betaH * H);
        H = max(z2, 0);

        % measure elapsed time
        elapsed_time = toc(start_time);

        % measure gradient calc count
        grad_calc_count = grad_calc_count + m*n;

        % update epoch
        epoch = epoch + 1;

        % store info
        infos = store_nmf_info(V, W, H, [], options, infos, epoch, grad_calc_count, elapsed_time);

        % display info
        display_info(method_name, epoch, infos, options);
    end

    x.W = W;
    x.H = H;

end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function [ a, u, v, seed ] = r8sm_random ( m, n, seed )

%*****************************************************************************80
%
%% R8SM_RANDOM randomizes a R8SM matrix.
%
%  Discussion:
%
%    The R8SM storage format is used for an M by N Sherman Morrison matrix B,
%    which is defined by an M by N matrix A, an M vector U, and
%    an N vector V, by B = A - U * V'
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns of the matrix.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real A(M,N), the R8SM matrix.
%
%    Output, real U(M), V(N), the R8SM vectors.
%
%    Output, integer SEED, an updated seed for the random number generator.
%
  for j = 1 : n
    for i = 1 : m
      [ a(i,j), seed ] = r8_uniform_01 ( seed );
    end
  end

  for i = 1 : m
    [ u(i), seed ] = r8_uniform_01 ( seed );
  end

  for j = 1 : n
    [ v(j), seed ] = r8_uniform_01 ( seed );
  end

  return
end


text:
function mb = trimeshMeanBreadth(vertices, faces)
%TRIMESHMEANBREADTH Mean breadth of a triangular mesh.
%
%   MB = trimeshMeanBreadth(VERTICES, FACES)
%   Computes the mean breadth (proporitonal to the integral of mean
%   curvature) of a triangular mesh.
%
%   Example
%     [V, F] = createCube;
%     F2 = triangulateFaces(F);
%     MB = trimeshMeanBreadth(V, F2)
%     MB =
%         1.5000
%
%   See also
%   meshes3d, trimeshSurfaceArea, trimeshEdgeFaces, polyhedronMeanBreadth
%
%   References
%   Stoyan D., Kendall W.S., Mecke J. (1995) "Stochastic Geometry and its
%       Applications", John Wiley and Sons, p. 26
%   Ohser, J., Muescklich, F. (2000) "Statistical Analysis of
%       Microstructures in Materials Sciences", John Wiley and Sons, p.352

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2015-08-19,    using Matlab 8.5.0.197613 (R2015a)
% Copyright 2015 INRA - Cepia Software Platform.

%% Check input validity

if size(faces, 2) ~= 3
    error('meshes3d:trimeshMeanBreadth:NonTriangularMesh', ...
        'Requires a triangular mesh as input');
end

%% Compute edge and edgeFaces arrays
% Uses the same code as in trimeshEdgeFaces

% compute vertex indices of each edge (in increasing index order)
edges = sort([faces(:,[1 2]) ; faces(:,[2 3]) ; faces(:,[3 1])], 2);

% create an array to keep indices of faces "creating" each edge
nFaces = size(faces, 1);
edgeFaceInds = repmat( (1:nFaces)', 3, 1);

% sort edges, keeping indices
[edges, ia, ib] = unique(edges, 'rows'); %#ok<ASGLU>
nEdges = size(edges, 1);

% allocate memory for result
edgeFaces = zeros(nEdges, 2);

% iterate over edges, to identify incident faces
for iEdge = 1:nEdges
    inds = find(ib == iEdge);
    edgeFaces(iEdge, 1:length(inds)) = edgeFaceInds(inds);
end

%% Compute dihedral angle for each edge

% compute normal of each face
normals = meshFaceNormals(vertices, faces);

% allocate memory for resulting angles
alpha = zeros(nEdges, 1);

% iterate over edges
for iEdge = 1:nEdges
    % indices of adjacent faces
    indFace1 = edgeFaces(iEdge, 1);
    indFace2 = edgeFaces(iEdge, 2);

    % normal vector of adjacent faces
    normal1 = normals(indFace1, :);
    normal2 = normals(indFace2, :);

    % compute dihedral angle of two vectors
    alpha(iEdge) = vectorAngle3d(normal1, normal2);
end

%% Compute mean breadth
% integrate the dihedral angles weighted by the length of each edge

% compute length of each edge
lengths = meshEdgeLength(vertices, edges);

% compute product of length by angles
mb = sum(alpha .* lengths) / (4*pi);


text:
function v=shear_bulk2poisson(mu,k)

v=((3.*k)-(2.*mu))./((6.*k)+(2.*mu));

end

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


text:
function copyFolderContents(pnamesrc, pnamedst)

copyfile([pnamesrc, '/*'], pnamedst);
copyfile([pnamesrc, '/*.*'], pnamedst);
copyfile([pnamesrc, '/.*'], pnamedst);

author: BUNPC


text:
classdef ThirdBodyGravityValidator < AbstractLaunchVehicleDataValidator
    %ThirdBodyGravityValidator Summary of this class goes here
    %   Detailed explanation goes here

    properties
        lvdData LvdData
    end

    methods
        function obj = ThirdBodyGravityValidator(lvdData)
            obj.lvdData = lvdData;
        end

        function [errors, warnings] = validate(obj)
            errors = LaunchVehicleDataValidationError.empty(0,1);
            warnings = LaunchVehicleDataValidationWarning.empty(0,1);

            warnEvtNoBodiesNums = [];
            warnEvtHasBodiesNums = [];
            evts = obj.lvdData.script.evts;
            for(i=1:length(evts))
                evt = evts(i);
                if(evt.propagatorObj == evt.forceModelPropagator)
                    forceModels = evt.forceModelPropagator.forceModels;
                    stateLogEntry = obj.lvdData.stateLog.getFirstStateLogForEvent(evt);

                    if(ismember(ForceModelsEnum.Gravity3rdBody, forceModels) && ...
                       isempty(stateLogEntry.thirdBodyGravity.bodies))
                        warnEvtNoBodiesNums(end+1) =  evt.getEventNum(); %#ok<AGROW>
                    end

                    if(not(ismember(ForceModelsEnum.Gravity3rdBody, forceModels)) && ...
                       not(isempty(stateLogEntry.thirdBodyGravity.bodies)))
                        warnEvtHasBodiesNums(end+1) =  evt.getEventNum(); %#ok<AGROW>
                    end
                end
            end

            if(not(isempty(warnEvtNoBodiesNums)))
                eventStr = makeEventsStr(unique(warnEvtNoBodiesNums));
                str = sprintf('Some events have the Third Body Gravity force model active but there are no third body gravity sources active in the event''s initial state. (Events: %s)', eventStr);
                warnings(end+1) = LaunchVehicleDataValidationWarning(str);
            end

            if(not(isempty(warnEvtHasBodiesNums)))
                eventStr = makeEventsStr(unique(warnEvtHasBodiesNums));
                str = sprintf('Some events have third body gravity sources active in the event''s initial state but the Third Body Gravity force model is inactive. (Events: %s)', eventStr);
                warnings(end+1) = LaunchVehicleDataValidationWarning(str);
            end
        end
    end
end


text:
function [leftpt,leftcurve,rightpt,rightcurve]=slicesurf3(node,elem,p1,p2,p3,step,minangle)
%
% [leftpt,leftcurve,rightpt,rightcurve]=slicesurf3(node,elem,p1,p2,p3,step,minangle)
%
% Slice a closed surface by a plane and extract the landmark nodes along
% the intersection between p1 and p3, then output into 2 segments: between
% p2 to p1 (left half), and p2 to p3 (right half)
%
% author: Qianqian Fang (q.fang at neu.edu)
%
% input:
%    node: an N x 3 array defining the 3-D positions of the mesh
%    elem: an N x 3 interger array specifying the surface triangle indices;
%    p1: 3D position of the start on the curve-of-interest
%    p2: 3D position of the middle on the curve-of-interest
%    p3: 3D position of the end on the curve-of-interest
%    step: (optional) a percentage (0-100) specifying the spacing of the
%        output landmark nodes; step=20 means the landmarks on the left
%        curve are spaced as 20% of the total lengths of the left-half, and
%        those on the right-curve are spaced at 20% of the right-half,
%        starting from p2.
%    minangle: (optional) a positive minangle will ask this function to
%        call polylinesimplify to remove sharp turns on the curve.
%
% output:
%    leftpt: the equal-spaced landmark nodes on the left-half (p2-p1)
%            intersection curve; spacing between these nodes are
%            (step% * length of the curve between p2-p1)
%    leftcurve: all nodes on the left-half (p2-p1) intersection curve
%    rightpt: the equal-spaced landmark nodes on the right-half (p2-p3)
%            intersection curve; spacing between these nodes are
%            (step% * length of the curve between p2-p3)
%    rightcurve: all nodes on the left-half (p2-p1) intersection curve
%
% -- this function is part of brain2mesh toolbox (http://mcx.space/brain2mesh)
%    License: GPL v3 or later, see LICENSE.txt for details
%

fullcurve=slicesurf(node, elem, [p1;p2;p3]);
if(nargin>=7 && minangle>0)
    fullcurve=polylinesimplify(fullcurve,minangle);
end

[fulllen, fullcurve]=polylinelen(fullcurve, p1,p3,p2);

[leftlen,  leftcurve]=polylinelen(fullcurve, p2, p1);
if(nargin>=6)
    [idx, weight, leftpt]=polylineinterp(leftlen, sum(leftlen)*(step:step:(100-step*0.5))*0.01, leftcurve);
else
    leftpt=leftcurve;
end

if(nargout>2)
    [rightlen, rightcurve]=polylinelen(fullcurve, p2, p3);
    if(nargin>=6)
        [idx, weight, rightpt]=polylineinterp(rightlen, sum(rightlen)*(step:step:(100-step*0.5))*0.01, rightcurve);
    else
        rightpt=rightcurve;
    end
end

author: fangq


text:
function [layers, net_output_info, dag_net]=gen_network_featnet_imgraw(featnet_config)

dag_net=dagnn.DagNN();
start_var_name='data_input';

layer_gen_info=[];
layer_gen_info.one_output_dim=3;
layer_gen_info.one_outputs={start_var_name};

layer_gen_info=do_add_img_input_scratch_path(dag_net, featnet_config, layer_gen_info);

net_output_info=[];
net_output_info.var_names=layer_gen_info.one_outputs;
net_output_info.var_dims=layer_gen_info.one_output_dim;

My_net_util.fix_padding_resnet(dag_net);

% input_vars=net.getInputs;
% dag_net.print({input_vars{1}, [400 400 3]}, 'Format', 'dot');

layers=cell(0);

one_layer=[];
one_layer.type='my_custom';
one_layer.custom_type='dagnn_wrapper';
one_layer.forward_fn=@cnn_layer_dagnn_wrapper_forward;
one_layer.backward_fn=@cnn_layer_dagnn_wrapper_backward;
one_layer.layer_update_fn=[];

one_layer.input_var_names=dag_net.getInputs();
one_layer.output_var_names=dag_net.getOutputs();

one_layer.use_single_output=true;

assert(length(one_layer.input_var_names)==1);
assert(length(one_layer.output_var_names)==1);
assert(all(ismember(net_output_info.var_names, one_layer.output_var_names)));

layers{end+1}=one_layer;

end

function layer_gen_info=do_add_img_input_scratch_path(dag_net, featnet_config, layer_gen_info)

layer_gen_info=add_scratch_img_input_block(dag_net, featnet_config, layer_gen_info);

conv_num=featnet_config.conv_num_one_block;
conv_block_num=featnet_config.conv_block_num;
for b_idx=1:conv_block_num

    one_output_dim=featnet_config.filter_num_blocks(b_idx);
    layer_name_prefix=sprintf('inputblock%d', b_idx);

    if b_idx>=2
        layer_gen_info=do_add_res_conv_block_downsample(dag_net, layer_gen_info, one_output_dim, conv_num, layer_name_prefix);
    else
        layer_gen_info=My_net_util.add_res_conv_block(dag_net, layer_gen_info, one_output_dim, conv_num, layer_name_prefix);
    end

end

output_var_name=layer_gen_info.one_outputs{1};
output_var_idx=dag_net.getVarIndex(output_var_name);
dag_net.vars(output_var_idx).fanout=0;

end

function layer_gen_info=add_scratch_img_input_block(dag_net, featnet_config, layer_gen_info)

input_dim=layer_gen_info.one_output_dim;
one_output_dim=featnet_config.filter_num_blocks(1);

filter_size=[3 3 input_dim one_output_dim];
one_inputs=layer_gen_info.one_outputs;
filter_name='imginput_c1';
one_outputs={[filter_name '_outvar']};
My_net_util.add_conv_layer_dagnn(dag_net, one_inputs, one_outputs, false, filter_size, filter_name, true);

tmp_l_idx=dag_net.getLayerIndex(filter_name);
dag_net.layers(tmp_l_idx).block.stride=2;

one_outputs=My_net_util.add_relu_dagnn(dag_net, filter_name, one_outputs);

layer_gen_info.one_output_dim=one_output_dim;
layer_gen_info.one_outputs=one_outputs;

filter_size=[3 3 one_output_dim one_output_dim];
one_inputs=layer_gen_info.one_outputs;
filter_name='imginput_c2';
one_outputs={[filter_name '_outvar']};
My_net_util.add_conv_layer_dagnn(dag_net, one_inputs, one_outputs, false, filter_size, filter_name, true);

layer_gen_info.one_outputs=one_outputs;

end

function layer_gen_info=do_add_res_conv_block_downsample(dag_net, layer_gen_info, output_dim, conv_num, layer_name_prefix)

one_output_dim=layer_gen_info.one_output_dim;
feat_dim_before=one_output_dim;
feat_dim_after=output_dim;
one_outputs_dimred=layer_gen_info.one_outputs;

[one_outputs_dimred, layer_name_dimred]=My_net_util.add_dim_reduce_layer(dag_net, one_outputs_dimred{1}, feat_dim_before, feat_dim_after);

layer_gen_info.one_output_dim=feat_dim_after;

tmp_l_idx=dag_net.getLayerIndex(layer_name_dimred);
dag_net.layers(tmp_l_idx).block.stride=2;

layer_gen_info.one_outputs=one_outputs_dimred;

layer_gen_info=My_net_util.add_res_conv_block(dag_net, layer_gen_info, output_dim, conv_num, layer_name_prefix);

end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function f = foldHOG(w)
% Condense HOG features into one orientation histogram.
%   f = foldHOG(w)
%
%   Used for displaying features and filters

% AUTORIGHTS
% -------------------------------------------------------
% Copyright (C) 2011-2012 Ross Girshick
% Copyright (C) 2008, 2009, 2010 Pedro Felzenszwalb, Ross Girshick
% Copyright (C) 2007 Pedro Felzenszwalb, Deva Ramanan
%
% This file is part of the voc-releaseX code
% (http://people.cs.uchicago.edu/~rbg/latent/)
% and is available under the terms of an MIT-like license
% provided in COPYING. Please retain this notice and
% COPYING if you use this file (or a portion of it) in
% your project.
% -------------------------------------------------------

% Return the contrast insensitive orientations
f = w(:,:,19:27);

% f=max(w(:,:,1:9),0)+max(w(:,:,10:18),0)+max(w(:,:,19:27),0);

author: rbgirshick

repo: voc-dpm

save_path: github-repos/MATLAB/rbgirshick-voc-dpm