text:
function d=distitpf(pf1,pf2,mode)
%DISTITPF calculates the Itakura spectral distance between power spectra D=(PF1,PF2,MODE)
%
% Inputs: PF1,PF2     Power spectra to be compared. Each row represents a power spectrum: the first
%                     and last columns represent the DC and Nyquist terms respectively.
%                     PF1 and PF2 must have the same number of columns.
%
%         MODE        Character string selecting the following options:
%                         'x'  Calculate the full distance matrix from every row of PF1 to every row of PF2
%                         'd'  Calculate only the distance between corresponding rows of PF1 and PF2
%                              The default is 'd' if PF1 and PF2 have the same number of rows otherwise 'x'.
%
% Output: D           If MODE='d' then D is a column vector with the same number of rows as the shorter of PF1 and PF2.
%                     If MODE='x' then D is a matrix with the same number of rows as PF1 and the same number of columns as PF2'.
%
% If ave() denotes the average over +ve and -ve frequency, the Itakura spectral distance is
%
%                               log(ave(pf1/pf2)) - ave(log(pf1/pf2))
%
% The Itakura distance is gain-independent, i.e. distitpf(g*pf1,pf2) is independent of g.

% The Itakura distance can also be calculated directly from AR coefficients; providing np is large
% enough, the values of d0 and d1 in the following will be very similar:
%
%         np=255; d0=distitar(ar1,ar2); d1=distitpf(lpcar2pf(ar1,np),lpcar2pf(ar2,np))
%

% Ref: A.H.Gray Jr and J.D.Markel, "Distance measures for speech processing", IEEE ASSP-24(5): 380-391, Oct 1976
%      L. Rabiner abd B-H Juang, "Fundamentals of Speech Recognition", Section 4.5, Prentice-Hall 1993, ISBN 0-13-015157-2
%      F. Itakura, "Minimum prediction residual principle applied to speech recognition", IEEE ASSP-23: 62-72, 1975

%      Copyright (C) Mike Brookes 1997
%      Version: $Id: distitpf.m 713 2011-10-16 14:45:43Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[nf1,p2]=size(pf1);
p1=p2-1;
nf2=size(pf2,1);
if nargin<3 | isempty(mode) mode='0'; end
if any(mode=='d') | (mode~='x' & nf1==nf2)
   nx=min(nf1,nf2);
   r=pf1(1:nx,:)./pf2(1:nx,:);
   q=log(r);
   d=log((sum(r(:,2:p1),2)+0.5*(r(:,1)+r(:,p2)))/p1)-(sum(q(:,2:p1),2)+0.5*(q(:,1)+q(:,p2)))/p1;
else
   r=permute(pf1(:,:,ones(1,nf2)),[1 3 2])./permute(pf2(:,:,ones(1,nf1)),[3 1 2]);
   q=log(r);
   d=log((sum(r(:,:,2:p1),3)+0.5*(r(:,:,1)+r(:,:,p2)))/p1)-(sum(q(:,:,2:p1),3)+0.5*(q(:,:,1)+q(:,:,p2)))/p1;
end

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% Tutorial for FAIR: elastic regularizer, matrix based and matrix free
%
% illustrates the tools for L2-norm based regularization
%
%  S(y) = alpha/2 * norm(B*(y-yRef)^2,
%
% where
%  alpha regularization parameter, weights regularization versus
%        distance in the joint objective function, alpha = 1 here
%  yRef  is a reference configuration, e.g. yRef = x or a
%        pre-registration result
%  B     a discretized partial differential operator either in explicit
%        matrix form or as a structure containing the necessary
%        parameters to compute B*y, here B = elastic
%==============================================================================

clear, close all, help(mfilename);

%% 2D example
omega  = [0,1,0,1]; % physical domin
m      = [16,12];   % number of discretization points
hd     = prod((omega(2:2:end)-omega(1:2:end))./m);

alpha  = 1;         % regularization parameter, irrelevant in this tutorial
mu     = 1;         % Lame constants, control elasticity properties
lambda = 0;         % Youngs modulus and Poisson ratio
para   = {'alpha',alpha,'mu',mu,'lambda',lambda};

yRef   = getStaggeredGrid(omega,m); % reference for regularization
yc     = rand(size(yRef));          % random transformation

%% 1. build elasticity operator on a staggered grid
%  - explicitly
%  - via regularizer.m
%  - via A(:,j) = B(e_j), B is the matrix free version

B = getElasticMatrixStg(omega,m,mu,lambda);
Amb = hd*alpha*(B'*B);

regularizer('reset','regularizer','mbElastic',para{:});
% regularizer('disp');
[Smb,dSmb,d2Smb] = regularizer(yc-yRef,omega,m);

regularizer('reset','regularizer','mfElastic',para{:});
% regularizer('disp');
[Scmf, dSmf, d2Smf] = regularizer(yc,omega,m);

% note: B*y can also be computed without storing B. The following loop
% computes A(:,j) = mfBy(e(j)), where mfBy is the matrix free operation B*y
% and e(j) is the j-th unit vector; hence A == B,

e = @(j) ((1:size(B,2))' == j); % j-th unit vector
A = sparse(size(B,1),size(B,2));
for j=1:size(B,2), A(:,j) = d2Smf.By(e(j),omega,m); end;

% now the transpose
e = @(j) ((1:size(B,1))' == j); % j-th unit vector
C = sparse(size(B,2),size(B,1));
for j=1:size(B,1), C(:,j) = d2Smf.BTy(e(j),omega,m); end;

% visualize matrices and transpose
FAIRfigure(1); clf;
subplot(2,3,1); spy(B);    title('B =elastic operator on staggered grid')
subplot(2,3,2); spy(A);    title('B from mfBy');
subplot(2,3,3); spy(B-A);  title('difference |B-mfBy()|');
subplot(2,3,4); spy(B');   title('B''')
subplot(2,3,5); spy(C);    title('B'' from mfBy');
subplot(2,3,6); spy(B'-C); title('difference |B''-mfBy(...,''BTy'')|');

%% 3D example
omega  = [0,4,0,2,0,1]; % physical domin
m      = [16,12,8];     % number of discretization points
hd     = prod((omega(2:2:end)-omega(1:2:end))./m);

yRef   = getStaggeredGrid(omega,m); % reference for regularization
yc     = rand(size(yRef));          % random transformation

% build elasticity operator on a staggered grid
B = getElasticMatrixStg(omega,m,mu,lambda);
Amb = hd*alpha*(B'*B);

regularizer('reset','regularizer','mbElastic',para{:});
% regularizer('disp');
[Smb,dSmb,d2Smb] = regularizer(yc-yRef,omega,m);

regularizer('reset','regularizer','mfElastic',para{:});
% regularizer('disp');
[Smf, dSmf, d2Smf] = regularizer(yc-yRef,omega,m);

zc     = randn(size(B,1),1);        % random for the adjoint

FAIRfigure(2); clf;
subplot(2,2,1); spy(B);          title('B elastic on staggered grid')
subplot(2,2,3); spy(Amb);        title('mb: hd\alpha B''*B')
subplot(2,2,4); spy(d2Smb);      title('mb: d2S')
subplot(2,2,2); spy(Amb-d2Smb);  title('difference')

Splain = 0.5*(yc-yRef)'*Amb*(yc-yRef);

fprintf('Splain = %-10s\n',num2str(Splain))
fprintf('Smb    = %-10s, Splain-Smb =%-10s\n',num2str(Smb),num2str(Splain-Smb))
fprintf('Smf    = %-10s, Splain-Smf =%-10s\n',num2str(Smf),num2str(Splain-Smf))

compare = @(s,l,r) fprintf('%-25s = %s\n',s,num2str(norm(l-r)));
compare('||B*yc-mfBy(y)||',B*yc,d2Smf.By(yc,omega,m));
compare('||B''*z-mfBy(z,''BTy'')||',B'*zc,d2Smf.BTy(zc,omega,m));

%% check elasticity

%% 3D example
omega  = [0,4,0,2,0,1]; % physical domin
m      = [32,16,8];     % number of discretization points
hd     = prod((omega(2:2:end)-omega(1:2:end))./m);

yRef   = getStaggeredGrid(omega,m); % reference for regularization
yc     = rand(size(yRef));          % random transformation

% build elasticity operator on a staggered grid
B = getElasticMatrixStg(omega,m,mu,lambda);

regularizer('reset','regularizer','mfElastic',para{:});
% regularizer('disp');
[Smf, dSmf, d2Smf] = regularizer(yc-yRef,omega,m);

H.d2S = d2Smf;
H.omega  = omega;
H.m      = m;
H.alpha  = regularizer('get','alpha');
H.mu     = regularizer('get','mu');
H.lambda = regularizer('get','lambda');
H.regularizer = regularizer;

By = B*yc;
testMF  = norm(By-d2Smf.By(yc,omega,m));
testMFT = norm(B'*By-d2Smf.BTy(By,omega,m));
fprintf('|By-mfBy|     = %s\n',num2str(testMF));
fprintf('|B''*By-mfBy''| = %s\n',num2str(testMFT));
% ---------------------------------------------------------------------
% using multigrid to solve A*uc = fc, where A = I + alpha*hd*B'*B

xc  = getStaggeredGrid(omega,m);
M   = speye(length(xc),length(xc)); % idenity matrix
A   = M + hd*alpha*B'*B;
fc  = A*xc;
u0  = zeros(size(xc));

% prepare for multigrid
H.MGlevel      = log2(m(1))+1;
H.MGcycle      = 4;
H.MGomega      = 2/3;   %% !!! 0.5 should be better
H.MGsmoother   = 'mfJacobi';
H.MGpresmooth  = 10;
H.MGpostsmooth = 10;
H.d2D.M        = full(diag(M));

uMG = mfvcycle(H,u0,fc,1e-12,H.MGlevel,5);

testMG = norm( A\fc - uMG);
fprintf('|A\\fc-uMG|=%s\n',num2str(testMG));

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

%==============================================================================
FAIRmessage(sprintf('<%s> done',mfilename));
%==============================================================================


text:
function [K, Kc] = sdrbfKernCompute(sdrbfKern, t1, t2)

% SDRBFKERNCOMPUTE Compute the SDRBF kernel given the parameters and t1.
% FORMAT
% DESC computes the kernel matrices for the switching dynamical radial
% basis function kernel. It is used together with the switching dynamical
% latent force model. The kernel matrices are stored in a cell structure
% and each component of the cell corresponds to the kernel matrix of a
% particular sequence of latent forces.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG t1 : the input points associated with the rows of the kernel.
% ARG t2 : the input points associated with the columns of the kernel.
% RETURN K : the kernel matrices computed at the given points.
% RETURN Kc : the kernel matrices grouped in cells. Suitable form for the
% sparse approximations.
%
% FORMAT
% DESC computes the kernel matrices for the switching dynamical radial
% basis function kernel. It is used together with the switching dynamical
% latent force model. The kernel matrices are stored in a cell structure
% and each component of the cell corresponds to the kernel matrix of a
% particular sequence of latent forces.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG t1 : input data matrix in the form of a design matrix.
% RETURN K : the kernel matrix computed at the given points.
% RETURN Kc : the kernel matrices grouped in cells. Suitable form for the
% sparse approximations.
%
% SEEALSO : rbfKernParamInit, kernCompute, kernCreate, rbfKernDiagCompute
%
% COPYRIGHT : Mauricio A. Alvarez, 2010

% KERN

if nargin < 3
    t2 = t1;
end

if size(t1, 2) > 1 || size(t2, 2) > 1
    error('Input can only have one column');
end

%Form the basic kernels
rbfKern = struct();

% Create structures that will make easy the computation of the kernels
% spVector = [sdrbfKern.switchingTimes t1(end)+0.1]; % For the last interval include the last point

spVector = [cumsum(sdrbfKern.switchingTimes) t1(end)+50];

dim1 = zeros(1, sdrbfKern.nIntervals);
dim2 = zeros(1, sdrbfKern.nIntervals);

for i =1:sdrbfKern.nIntervals
    for j=1:sdrbfKern.nlfPerInt
        rbfKern(i,j).variance = sdrbfKern.variance;
        rbfKern(i,j).inverseWidth = sdrbfKern.inverseWidth(j,i);
        rbfKern(i,j).isNormalised = sdrbfKern.isNormalised;
    end
    newt1 = t1(t1> spVector(i) & t1<spVector(i+1));
    newt2 = t2(t2> spVector(i) & t2<spVector(i+1));
    dim1(i) = length(newt1);
    dim2(i) = length(newt2);
end

if sum(dim1)~=length(t1) || sum(dim2)~=length(t2)
    error('A problem with the dimensions of the switching intervals occured')
end

% Compute some necessary constants

K = cell(sdrbfKern.nlfPerInt,1);
Kc = cell(sdrbfKern.nlfPerInt,1);

for q=1:sdrbfKern.nlfPerInt
    startValOne = 1;
    endValOne   = 0;
    startValThree = 1;
    endValThree = 0;
    tK = zeros(sum(dim1), sum(dim2));
    for i=1:sdrbfKern.nIntervals
        endValOne = endValOne + dim1(i);
        endValThree = endValThree + dim2(i);
        tK(startValOne:endValOne, startValThree:endValThree) = rbfKernCompute(rbfKern(i,q),...
            t1(startValOne:endValOne) - spVector(i), t2(startValThree:endValThree) - spVector(i));
        Kc{q}{i} = tK(startValOne:endValOne, startValThree:endValThree);
        startValThree = endValThree + 1;
        startValOne = endValOne + 1;
    end
    tK = real(tK);
    K{q} = tK;
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [axis,angle]=rotation2axisangle(R)
cosq = 0.5*(R(1,1)+R(2,2)+R(3,3)-1);
if(cosq>1)
    cosq=1;
end
if(cosq<-1)
    cosq=-1;
end
angle = acos(cosq);
axis = zeros(3,1);
if(abs(angle)<1e-4)
elseif(abs(angle-pi)<1e-4)
    for i=1:3
        a = 0.5*(R(i,i)+1);
        if abs(a)<1e-4
            axis(i) = 0;
        else
            axis(i) = a^0.5;
        end
    end

    maxind = 1;
    if axis(2)>axis(1)
        maxind = 2;
    end
    if axis(3)>axis(2)
        maxind = 3;
    end
    if(maxind == 1)
        a1a2 = R(1,2) + R(2,1);
        a1a3 = R(1,3) + R(3,1);
        if(a1a2<0)
            axis(2) = -axis(2);
        end
        if(a1a3<0)
            axis(3) = -axis(3);
        end
    end
    if(maxind == 2)
        a1a2 = R(1,2) + R(2,1);
        a2a3 = R(2,3) + R(3,2);
        if(a1a2<0)
            axis(1) = -axis(1);
        end
        if(a2a3<0)
            axis(3) = -axis(3);
        end
    end
    if(maxind == 3)
        a1a3 = R(1,3) + R(3,1);
        a2a3 = R(2,3) + R(3,2);
        if(a1a3<0)
            axis(1) = -axis(1);
        end
        if(a2a3<0)
            axis(2) = -axis(2);
        end
    end
    axis = axis/norm(axis);
else
    axis(1) = R(3,2) - R(2,3);
    axis(2) = R(1,3) - R(3,1);
    axis(3) = R(2,1) - R(1,2);
    axis = axis/norm(axis);
end

end

author: xuhuairuogu

repo: V-REP-Simulation-Projects

path: github-repos/MATLAB/xuhuairuogu-V-REP-Simulation-Projects/V-REP-Simulation-Projects-841b944af4ea3a8fb250578d36434515f577f411/admittance control(adapted from an example in the book--A Systematic Approach to Learning Robot Programming with ROS)/iiwa14_kinematics/rotation2axisangle.m


text:
function element_num = sphere_grid_q16_element_num ( nelemx, nelemy )

%*****************************************************************************80
%
%% SPHERE_GRID_Q16_ELEMENT_NUM counts the elements in a Q16 sphere grid.
%
%  Example:
%
%    Input:
%
%      NELEMX = 3, NELEMY = 2
%
%    Output:
%
%      ELEMENT_NUM = NELEMX * NELEMY = 6
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NELEMX, NELEMY, the number of elements along the
%    X and Y directions.
%
%    Output, integer ELEMENT_NUM, the number of elements in the grid.
%
  element_num = nelemx * nelemy;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function i4_factorial_test ( )

%*****************************************************************************80
%
%% I4_FACTORIAL_TEST tests I4_FACTORIAL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 October 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4_FACTORIAL_TEST:\n' );
  fprintf ( 1, '  I4_FACTORIAL evaluates the factorial function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     X       Exact F       I4_FACTORIAL(X)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, n, fn ] = i4_factorial_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fn2 = i4_factorial ( n );

    fprintf ( 1, '  %2d  %10d  %10d\n', n, fn, fn2 );

  end

  return
end


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   PARAMETERS Returns a data structure containing the parameters of the
%   example 4 DOF planar robot.
%
%   Author: Arturo Gil. Universidad Miguel Hernandez de Elche.
%   email: arturo.gil@umh.es date:   05/03/2012
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.
function robot = parameters()

robot.name='Example 4DOF planar arm';

%kinematic data DH parameters
robot.DH.theta='[q(1) q(2) q(3) q(4)]';
robot.DH.d='[0  0  0  0]';
robot.DH.a='[1  1  1  1]';
robot.DH.alpha='[0  0  0  0]';

%number of degrees of freedom
robot.DOF = 4;

%Jacobian matrix
% if not defined, compute it with manipulator_jacobian
robot.J=[];
%the proper rows of J to compute manipulability
robot.selJ =[1,2,6];
robot.kind=['R' 'R' 'R' 'R'];

%Function name to compute inverse kinematic
robot.inversekinematic_fn = 'inverse_kinematics_4dofplanar(robot, T, q)';
robot.directkinematic_fn = 'directkinematic(robot, q)';

%minimum and maximum rotation angle in rad
robot.maxangle =[deg2rad(-180) deg2rad(180); %Axis 1, minimum, maximum
                deg2rad(-180) deg2rad(180);
                deg2rad(-180) deg2rad(180)]; %Axis 2, minimum, maximum

            %minimum and maximum rotation angle in rad
robot.maxangle =[deg2rad(-90) deg2rad(90); %Axis 1, minimum, maximum
                deg2rad(-90) deg2rad(90);
                deg2rad(-180) deg2rad(180)]; %Axis 2, minimum, maximum

%maximum absolute speed of each joint rad/s or m/s
robot.velmax = []; %empty, not available

robot.accelmax=robot.velmax/0.1; % 0.1 is here an acceleration time
% end effectors maximum velocity
robot.linear_velmax = 0; %m/s, example, not available

%base reference system
robot.T0 = eye(4);

%INITIALIZATION OF VARIABLES REQUIRED FOR THE SIMULATION
%position, velocity and acceleration
robot=init_sim_variables(robot);
robot.path = pwd;

%this is an ad hoc transformation between the piece and the robot's end
%effector
% robot.Tcoupling=[-1 0 0  0;
%      0 1 0  0;
%      0 0 -1 0;
%      0 0 0 1];
% Tc1 = [cos(-pi/2) 0 sin(-pi/2) 0;
%        0 1 0 0;
%        -sin(-pi/2) 0 cos(-pi/2) 0;
%        0  0  0  1       ];
% Tc2 = [1 0 0 0;
%        0 cos(pi) -sin(pi) 0;
%        0  sin(pi) cos(pi) 0;
%        0 0 0 1];
% robot.Tcoupling =Tc1*Tc2;
robot.Tcoupling = eye(4);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%GRAPHICS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%read graphics files
robot.graphical.has_graphics=1;
robot.graphical.color = [255 20 40]./255;
%for transparency
robot.graphical.draw_transparent=0;
%draw DH systems
robot.graphical.draw_axes=1;
%DH system length and Font size, standard is 1/10. Select 2/20, 3/30 for
%bigger robots
robot.graphical.axes_scale=1;
%adjust for a default view of the robot
robot.axis=[-3.5 3.5 -3.5 3.5 0 1]
robot = read_graphics(robot);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%DYNAMIC PARAMETERS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
robot.has_dynamics=1;

%consider friction in the computations
robot.dynamics.friction=0;

%link masses (kg)
robot.dynamics.masses=[1 1 1] ;

%COM of each link with respect to own reference system
robot.dynamics.r_com=[-0.5      0         0; %(rx, ry, rz) link 1
                      -0.5      0         0;
                      -0.5      0         0];%(rx, ry, rz) link 2

%link masses
m1 = robot.dynamics.masses(1);
m2 = robot.dynamics.masses(2);
m3 = robot.dynamics.masses(3);

a=eval(robot.DH.a);
L1 = a(1);
L2 = a(2);
L3 = a(3);

%Momentos de inercia de cada eslabon.
% Ixx	Iyy	Izz	Ixy	Iyz	Ixz, por cada fila
robot.dynamics.Inertia=[0   m1*L1^2/3   m1*L1^2/3    0	0	0;
                        0   m2*L2^2/3   m2*L2^2/3    0	0	0;
                        0   m3*L3^2/3   m3*L3^2/3    0	0	0];

%Actuator rotor inertia
robot.motors.Inertia=[0 0 0];
%Reduction ratio motor/joint speed
robot.motors.G=[1 1 1];
%consider friction
robot.friction=0;
%Viscous friction factor, motor referred
%robot.B = [1e-3  1e-3 1e-3];
robot.motors.Viscous = [10  10 10];
%Coulomb friction, motor referred
robot.motors.Coulomb = [0 0;
            0 0;
            0 0];

%SPECIAL PARAMETERS TO SOLVE THE INVERSE KINEMATICS
robot.parameters.step_time=0.1;
%Error in XYZ to stop inverse kinematics
robot.parameters.epsilonXYZ=0.01;
%Error in Quaternion to stop inverse kinematics.
robot.parameters.epsilonQ=0.01;
robot.parameters.stop_iterations=500;

repo: ARTE


text:
% test_remove_vertex_from_tri.m

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2013 University of Oxford
% Version: 0.1.0
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% 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. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see
% <http://www.gnu.org/licenses/>.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Remove opposite vertices from a cube
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% create a surface mesh that is a cube
x = [
    0 0 0
    0 0 1
    0 1 0
    0 1 1
    1 0 0
    1 0 1
    1 1 0
    1 1 1
    ];

tri = [
    1     3     5
    2     1     5
    2     3     1
    4     3     2
    6     2     5
    4     2     6
    7     6     5
    3     7     5
    4     7     3
    8     6     7
    4     8     7
    4     6     8
    ];

% plot mesh
hold off
trisurf(tri, x(:,1), x(:,2), x(:,3), ones(1, size(x, 1)))
axis equal

% vertices to remove
idx = [2 7];

% plot vertices that are going to be removed
hold on
plot3(x(idx, 1), x(idx, 2), x(idx, 3), 'ro')

% remove vertices
[tri2, x2] = remove_vertex_from_tri(tri, x, idx);

% plot result
hold off
trisurf(tri2, x2(:,1), x2(:,2), x2(:,3), ones(1, size(x2, 1)))
hold on
plot3(x(idx, 1), x(idx, 2), x(idx, 3), 'ro')
axis equal

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Remove neighbour vertices from a cube
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% create a surface mesh that is a cube
x = [
    0 0 0
    0 0 1
    0 1 0
    0 1 1
    1 0 0
    1 0 1
    1 1 0
    1 1 1
    ];

tri = [
    1     3     5
    2     1     5
    2     3     1
    4     3     2
    6     2     5
    4     2     6
    7     6     5
    3     7     5
    4     7     3
    8     6     7
    4     8     7
    4     6     8
    ];

% plot mesh
hold off
trisurf(tri, x(:,1), x(:,2), x(:,3), ones(1, size(x, 1)))
axis equal

% vertices to remove
idx = [1 2];

% plot vertices that are going to be removed
hold on
plot3(x(idx, 1), x(idx, 2), x(idx, 3), 'ro')

% remove vertices
[tri2, x2] = remove_vertex_from_tri(tri, x, idx);

% plot result
hold off
trisurf(tri2, x2(:,1), x2(:,2), x2(:,3), ones(1, size(x2, 1)))
hold on
plot3(x(idx, 1), x(idx, 2), x(idx, 3), 'ro')
axis equal

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function draw_loudspeakers(x0,dimensions,conf)
%DRAW_LOUDSPEAKERS draws loudspeaker symbols or "o" at the given positions
%
%   Usage: draw_loudspeakers(x0,[dimensions],conf)
%
%   Input parameters:
%       x0          - positions and directions of the loudspeakers / m
%       dimensions  - dimension defining the plane in which the loudspeaker
%                     symbol should be plotted. For example [1 1 0] corresponds
%                     to the xy-plane. Default: [1 1 0]
%       conf        - configuration struct (see SFS_config)
%
%   DRAW_LOUDSPEAKERS(x0,dimensions,conf) draws loudspeaker symbols or filled
%   points at the given secondary source positions. This can be controlled by
%   the conf.plot.realloudspeakers setting. The loudspeaker symbols are pointing
%   in their given direction and are colored accordingly to their weights,
%   whereby the weights are scaled to an absolute maximum of 1 before.
%
%   See also: plot_sound_field

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% 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.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input parameter =====================================
nargmin = 2;
nargmax = 3;
narginchk(nargmin,nargmax);
isargsecondarysource(x0)
nls = size(x0,1);
if nargin<nargmax
    conf = dimensions;
    dimensions = [1 1 0];
end

%% ===== Configuration ===================================================
p.realloudspeakers = conf.plot.realloudspeakers;
p.lssize = conf.plot.lssize;

%% ===== Plotting ========================================================
% Check in which plane we should plot the secondary sources and shift them
% around accordingly
if ~any(dimensions)
    return;
elseif ~dimensions(1)
    x0(:,1:2) = x0(:,2:3);
elseif ~dimensions(2)
    x0(:,2) = x0(:,3);
end

% Plot only "o" at the loudspeaker positions, use this as default for all cases
% that are not the x-y-plane
if ~p.realloudspeakers || ~(dimensions(1)&&dimensions(2))

    if p.realloudspeakers && ~(dimensions(1)&&dimensions(2))
        warning('%s: Real loudspeaker can only be drawn in the x-y-plane', ...
            upper(mfilename));
    end
    % Plot all secondary sources with the same color + symbol
    plot(x0(:,1),x0(:,2),'o', ...
        'MarkerFaceColor','k', ...
        'MarkerEdgeColor','k', ...
        'MarkerSize',4);

else

    % Weightings of the single loudspeakers (scale maximum to 1)
    weights = x0(:,7) / max(abs(x0(:,7)));
    % Set fill color for active loudspeakers
    % Scale the color: sc = 1 => black, sc = 0.5 => gray.
    sc = 0.6;
    fc = [(1-sc.*weights), ...
          (1-sc.*weights), ...
          (1-sc.*weights)];

    % Plot a real speaker symbol at the desired position
    % vertex coordinates with height and width after lssize
    w = p.lssize;
    h = p.lssize;
    v1 = [-h -h -h/2 -h/2 -h ; -w/2 w/2 w/2 -w/2 -w/2 ; 0 0 0 0 0];
    v2 = [-h/2 0 0 -h/2 ; -w/6 -w/2 w/2 w/6 ; 0 0 0 0 ];

    hold on;

    % Draw loudspeakers
    for n=1:nls

        % Get the azimuth direction of the secondary sources
        [phi,~] = cart2pol(x0(n,4),x0(n,5));

        % Rotation matrix (orientation of the speakers)
        % R = [cos(phi(n)) -sin(phi(n));sin(phi(n)) cos(phi(n))];
        R = rotation_matrix(phi);
        for k=1:length(v1)
            vr1(:,k) = R * v1(:,k);
        end
        for k=1:length(v2)
            vr2(:,k) = R * v2(:,k);
        end

        % Shift
        v01(1,:) = vr1(1,:) + x0(n,1);
        v01(2,:) = vr1(2,:) + x0(n,2);
        v02(1,:) = vr2(1,:) + x0(n,1);
        v02(2,:) = vr2(2,:) + x0(n,2);

        % Draw loudspeakers
        fill(v01(1,:),v01(2,:),fc(n,:));
        fill(v02(1,:),v02(2,:),fc(n,:));

    end

    hold off;

end

% Set equal axis ratio
axis image

author: sfstoolbox

repo: sfs-matlab


text:
clear all
% explicit euler
% m = 1; %adjust value of m
% y(1) = 1;%input your initial condition
% dt = 0.2; %adjust your step size
% T = 0:dt:15; %set up your time domain, here I have [0,5]
% for i = 2:length(T) %construct a 'for' loop
% y(i) = y(i-1) + m*dt*(-3*(i-2)/(i-1))*y(i-1)+m*dt*(2*(1+(i-2)^3))*exp(-(i-2));
% end
% plot(T,y)
%
% Implicit euler
% n = 1; %adjust value of m
% w(1) = 1;%input your initial condition
% dtt = 0.2; %adjust your step size
% T1 = 0:dtt:15; %set up your time domain, here I have [0,5]
% for j = 2:length(T1) %construct a 'for' loop
%      w(j) = (w(j-1)+ dtt*(2*(j^3))*exp(-j+1))/(1+3*dtt*(j-1)/j) ;
% end
% plot(T1,w)
%
%
% Crank-Nicolson
m = 1; %adjust value of m
y(1) = 1;%input your initial condition
dt = 0.2; %adjust your step size
T = 0:dt:15; %set up your time domain, here I have [0,5]
for i = 2:length(T) %construct a 'for' loop
    y(i)=(dt*i^3*exp(-i+1)+dt*(i-1)^3*exp(-i+2)+...
           y(i-1)*(1-0.5*dt*(-3*i-6/(i-1))))/(1+0.5*dt*(3*i+3/i));
end
plot(T,y)

repo: Aero-matlab


text:
close all;
clearvars;
clc;
rng default;
N = 32;
K = 4;
gen = spx.data.synthetic.SparseSignalGenerator(N, K);
rep =  gen.biUniform();
figure;
stem(rep, '.');
export_fig images/demo_sparse_biuniform_1.png -r120 -nocrop;

rep =  gen.biUniform(2, 4);
figure;
stem(rep, '.');
export_fig images/demo_sparse_biuniform_2.png -r120 -nocrop;

author: indigits

repo: sparse-plex


text:
clear all; close all; clc

load catData_w.mat; load dogData_w.mat; CD=[dog_wave cat_wave];
train=[dog_wave(:,1:60) cat_wave(:,1:60)];
test=[dog_wave(:,61:80) cat_wave(:,61:80)];
label=[ones(60,1); -1*ones(60,1)].';

A=label*pinv(train); test_labels=sign(A*test);
subplot(4,1,1), bar(test_labels,'FaceColor',[.6 .6 .6],'EdgeColor','k'), axis off
subplot(4,1,2), bar(A,'FaceColor',[.6 .6 .6],'EdgeColor','k'), axis([0 1024 -0.002 0.002]), axis off
figure(2), subplot(2,2,1)
A2=flipud(reshape(A,32,32)); pcolor(A2), colormap(gray), axis off

figure(1), subplot(4,1,3)
A=lasso(train.',label.','Lambda',0.1).';
test_labels=sign(A*test);
bar(test_labels,'FaceColor',[.6 .6 .6],'EdgeColor','k'), axis off
subplot(4,1,4)
bar(A,'FaceColor',[.6 .6 .6],'EdgeColor','k'), axis([0 1024 -0.008 0.008]), axis off
figure(2), subplot(2,2,2)
A2=flipud(reshape(A,32,32)); pcolor(A2), colormap(gray), axis off

author: dynamicslab


text:
%% Example of de-meaning
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

%% Generate random dataset
ds=cosmo_synthetic_dataset('nchunks',4,'ntargets',3);

% add some constant to all data
ds.samples=ds.samples+2;

% show dataset
subplot(2,2,1);
imagesc(ds.samples,[-4 4])
title('before demeaning');
subplot(2,2,2);
hist(ds.samples(:),10)
xlim([-6 6]);

%% Split the dataset by chunks
% >@@>
splits=cosmo_split(ds,{'chunks'},1);
% <@@<
nsplits=numel(splits);

% allocate space for output
outputs=cell(nsplits,1);

% treat each element in splits seperately, and subtract the mean for each
% feature seperately
for k=1:nsplits
    d=splits{k};
    % >@@>

    % mean over samples, for each feature
    mu=mean(d.samples,1);

    % subtract the mean.
    % equivalent, but less efficient, is:
    %     nsamples=size(d.samples,1);
    %     d.samples=d.samples-repmat(mu,nsamples,1);
    %
    d.samples=bsxfun(@minus,d.samples,mu);
    % <@@<

    % store output
    outputs{k}=d;
end

ds_demeaned=cosmo_stack(outputs);

% show dataset
subplot(2,2,3);
imagesc(ds_demeaned.samples,[-4 4])
title('after demeaning');
subplot(2,2,4);
hist(ds_demeaned.samples(:),10);
xlim([-6 6]);

%% Alternative approach to demeaning

% note: the samples in the output are in a different order than the input,
% but otherwise the same
demeaner=@(x)bsxfun(@minus,x,mean(x,1)); % function handle as helper
ds_demeaned_alt=cosmo_fx(ds,demeaner,'chunks');

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
% Script: call_bootfitF.m
% -----------------------
% Script to input parameters for a single Omori law fit
% Calculation is done in function calc_bootfitF
% The function works on the catalog newt2!!!
% J.Woessner
% last update: 10.07.03

report_this_filefun(mfilename('fullpath'));
% Get input parameters
prompt  = {,'Enter length of learning period (days)','Enter forecast period (days):','Enter number of bootstraps:'};
title   = 'Parameters ';
lines= 1;
def     = {'50','50','100'};
answer  = inputdlg(prompt,title,lines,def);
time = str2double(answer{1});
timef = str2double(answer{2});
bootloops = str2double(answer{3});

% maepi is the mainshock of sequence
calc_bootfitF(newt2,time,timef,bootloops,maepi)

author: CelsoReyes


text:
  function fs = de_ftab_xform(ftab, fm)
%|function fs = de_ftab_xform(ftab, fm)
%|
%| convert nonlinear BH function "f" to "f*" using jacobian
%| in
%|	ftab	struct	from de_ftab_build()
%|	fm	[n? M]	f_m samples
%| out
%|	fs	[n? L]	f^*_l samples
%|
%| Copyright 2002-2-15, Jeff Fessler, University of Michigan

if nargin < 2, ir_usage, end

[LL MM] = size(ftab.T);
dim = size(fm);
if MM == 1
	dim = [dim 1];
end
if dim(end) ~= MM, keyboard, error 'bad input dim', end
fm = reshape(fm, prod(dim(1:end-1)), MM); % [*n M]
fs = fm * ftab.T'; % [*n L]
fs = reshape(fs, [dim(1:end-1) LL]); % [n? L]

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function showsteerable(coeff)
% Show the subbands of steerable pyramid

figure();
ht=length(coeff);

for i=2:length(coeff)-1
    for j=1:4
        subplot(ht-2,4,(i-2)*4+j);
        imshow(coeff{i}{j},[]);
    end
end

author: Sable

repo: mcbench-benchmarks


text:
function xyr = Xr(s)
x = 1+2*s-2*s^2 ;
y = s ;

xyr = [x ; y] ;


text:
function [sig,a] = spm_shoot_blur(t,prm,its,sig)
% A function for blurring ("smoothing") tissue probability maps
% FORMAT [sig,a_new] = spm_shoot_blur(t,prm,its,sig)
%     t   - sufficient statistics
%     prm - regularisation parameters (1,1,1, 0.01,0.02,1)
%     its - max no. iterations (12)
%     sig - optional starting estimates
%
%     sig - "smoothed" average
%     a   - parameters
%
% The core of this procedure is described in:
%     John Ashburner & Karl J. Friston.
%     "Computing Average Shaped Tissue Probability Templates"
%     NeuroImage, In Press, Accepted Manuscript, Available online 24 December 2008
%
% However, there is an additional modification such that the the null space
% of the parameters is rotated out.
%________________________________________________________
% (c) Wellcome Trust Centre for NeuroImaging (2009)

% John Ashburner
% $Id: spm_shoot_blur.m 7387 2018-08-03 15:13:57Z john $

d   = [size(t),1,1,1];
if nargin<3, its = 16;                         end % Maximum no. iterations
if nargin<2, prm = [1,1,1, 0.01 0.02 1];       end % Default regularisation
rits = [1 1]; % No. cycles and no. relaxation iterations

W    = zeros([d(1:3) round(((d(4)-1)*d(4))/2)],'single'); % 2nd derivatives
gr   = zeros([d(1:3),d(4)-1],'single');                   % 1st derivatives

% Re-organise sufficient statistics to a form that is easier to work with
t    = max(t,eps('single')*1000);
s    = sum(t,4);
for k=1:d(4)
    t(:,:,:,k) = t(:,:,:,k)./s;
end
maxs   = max(s(:)); % Used for scaling the regularisation
prm(4) = prm(4)+maxs*d(4)*1e-6;

% Only d(4)-1 fields need to be estimated because sum(a,4) = 0.  This matrix
% is used to rotate out the null space
R     = null(ones(1,d(4)));

% Initial starting estimates (if sig is passed)
a     = zeros([d(1:3),d(4)-1],'single');
if nargin>=4
    for z=1:d(3) % Loop over planes
        sz = sig(:,:,z,:);
        sz = min(max(sz,0),1);
        sz(~isfinite(sz)) = 1/d(4);
        sz = squeeze(log(double(sz*(1-d(4)*1e-3)+1e-3)));
        for j1=1:(d(4)-1)
            az = zeros(d(1:2));
            for j2=1:d(4)
                az = az + R(j2,j1)*sz(:,:,j2); % Note the rotation
            end
            a(:,:,z,j1) = az;
        end
        clear sz az
    end
end

for i=1:its

    ll  = 0;
    for z=1:d(3) % Loop over planes

        % Compute softmax for this plane
        sig = double(reshape(sftmax(a(:,:,z,:),R),[d(1:2),d(4)]));

        % -ve log likelihood of the likelihood
        ll  = ll - sum(sum(sum(log(sig).*reshape(t(:,:,z,:),[d(1:2),d(4)]),3).*s(:,:,z)));

        % Compute first derivatives (d(4)-1) x 1
        grz = sig - double(reshape(t(:,:,z,:),[d(1:2),d(4)]));
        for j1=1:(d(4)-1)
            gr(:,:,z,j1) = 0;
            for j2=1:d(4)
                gr(:,:,z,j1) = gr(:,:,z,j1) + R(j2,j1)*grz(:,:,j2); % Note the rotation
            end
            gr(:,:,z,j1) = gr(:,:,z,j1).*s(:,:,z);
        end

        % Compute d(4) x d(4) matrix of second derivatives at each voxel.
        % These should be positive definate, but rounding errors may prevent this.
        % Regularisation is included to enforce +ve definateness.
        wz = zeros([d(1:2),d(4),d(4)]);
        for j1=1:d(4)
            wz(:,:,j1,j1) =   (1-sig(:,:,j1)).*sig(:,:,j1).*s(:,:,z);
            for j2=1:(j1-1)
                wz(:,:,j1,j2) = -sig(:,:,j1) .*sig(:,:,j2).*s(:,:,z);
                wz(:,:,j2,j1) = wz(:,:,j1,j2);
            end
        end

        % First step of rotating 2nd derivatives to (d(4)-1) x (d(4)-1)
        % by R'*W*R
        wz1 = zeros([d(1:2),d(4),d(4)-1]);
        for j1=1:d(4)
            for j2=1:(d(4)-1)
                tmp = zeros(d(1:2));
                for j3=1:d(4)
                    tmp = tmp + wz(:,:,j1,j3)*R(j3,j2);
                end
                wz1(:,:,j1,j2) = tmp;
            end
        end

        % Second step of rotating 2nd derivatives to (d(4)-1) x (d(4)-1)
        % by R'*W*R
        wz = zeros([d(1:2),d(4)-1,d(4)-1]);
        for j1=1:(d(4)-1)
            for j2=1:(d(4)-1)
                tmp = zeros(d(1:2));
                for j3=1:d(4)
                    tmp = tmp + R(j3,j1)*wz1(:,:,j3,j2);
                end
                wz(:,:,j1,j2) = tmp;
            end
        end

        % First pull out the diagonal of the 2nd derivs
        for j1=1:d(4)-1
            W(:,:,z,j1) = wz(:,:,j1,j1);% + maxs*sqrt(eps('single'))*d(4)^2;
        end

        % Then pull out the off diagonal parts (note that matrices are symmetric)
        jj = d(4);
        for j1=1:d(4)-1
           for j2=(j1+1):(d(4)-1)
               W(:,:,z,jj) = wz(:,:,j2,j1);
               jj = jj+1;
           end
        end
    end

    % ss1 and ss2 are for examining how close the 1st derivatives are to zero.
    % At convergence, the derivatives from the likelihood term should match those
    % from the prior (regularisation) term.
    ss1 = sum(sum(sum(sum(gr.^2))));
    gr1 = spm_field('vel2mom',a,prm);        % 1st derivative of the prior term
    ll1 = 0.5*sum(sum(sum(sum(gr1.*a)))); % -ve log probability of the prior term
    gr  = gr + gr1;                       % Combine the derivatives of the two terms
    ss2 = sum(sum(sum(sum(gr.^2))));      % This should approach zero at convergence
    mx  = max(max(max(sum(gr.^2,4))));

    fprintf('%2d %8.4f %8.4f %8.4f %g\n', i, ll/prod(d(1:3)),ll1/prod(d(1:3)), (ll+ll1)/prod(d(1:3)), (ss2)/prod(d(1:3)));

    reg = double(0.01*sqrt(mx)*d(4));
   %reg = double(0.1*sqrt(ss2/prod(d(1:3))));
    a   = a - spm_field(W,gr,[prm(1:3) prm(4)+reg prm(5:6) rits]); % Gauss-Newton update

    if ss2/ss1<1e-4, break; end        % Converged?
end

sig = sftmax(a,R);
%________________________________________________________

%________________________________________________________
function sig = sftmax(a,R)
% Softmax function

d     = [size(a) 1 1 1];
sig   = zeros([d(1:3),d(4)+1],'single');
trunc = log(realmax('single')*(1-eps('single'))/(d(4)+1));

for j=1:size(a,3) % Loop over planes

    % Rotate the null-space back in to the data
    aj  = double(reshape(a(:,:,j,:),[d(1:2),d(4)]));
    sj  = zeros([d(1:2),d(4)+1]);
    for j1=1:d(4)+1
        sj(:,:,j1) = 0;
        for j2=1:d(4)
            sj(:,:,j1) = sj(:,:,j1) + R(j1,j2)*aj(:,:,j2);
        end
    end

    % Compute softmax
    sj = min(max(sj,-trunc),trunc);
    sj = exp(sj)+eps('single')*(d(4)+1);
    s  = sum(sj,3);
    for i=1:d(4)+1
        sig(:,:,j,i) = single(sj(:,:,i)./s);
    end

end
%________________________________________________________

author: spm


text:
function [Y,W,SetupStruc] = Process_DSB(s,Transfer,SetupStruc)
K = SetupStruc.DSB.K;
hop = SetupStruc.DSB.hop;
win = hanning(K,'periodic');
win = win/sqrt(sum(win(1:hop:K).^2));
SetupStruc.DSB.win = win;  % Preserve 'win' in 'SetupStruc'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = size(s,2);
for i = 1:N
    X(:,:,i) = fft(enframe(s(:,i),win,hop)');
end
frame_N = size(X,2);
K_m = K/2+1;
Num = size(Transfer,3);
Y = zeros((frame_N-1)*hop+K,Num);
Y_f = zeros(size(X,1),size(X,2),Num);
%%%%%%%%%%%%%%%%%%%%%%%%%% Obtain processing matrix 'W'
W = conj(permute(Transfer,[3 2 1]))/N;
for i = 2:K_m
    X_f = permute(X(i,:,:),[3 2 1]);
    W_f = permute(W(:,:,i),[1 2 3]);
    Y_ = W_f*X_f;
    Y_f(i,:,:) = Y_.';
    if(i~=K_m)
        Y_f(K+2-i,:,:) = Y_';
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Recover signals
if(K/hop==2)
    win = ones(K,1);
end
for i = 1:Num
    Y(:,i) = overlapadd(real(ifft(Y_f(:,:,i)))',win,hop);
end

return;

repo: Voice-Separation-and-Enhancement


text:
function [ m, d ] = mothers_day ( y )

%*****************************************************************************80
%
%% MOTHERS_DAY computes the date of Mother's Day (US) for a Common year.
%
%  Discussion:
%
%    Mother's Day occurs on the second Sunday in May.
%
%  Example:
%
%    Input:
%
%      Y = 2003
%
%    Output:
%
%      M = 5
%      D = 11
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, the year.
%
%    Output, integer M, D, the month and day of Mother's Day.
%

%
%  Determine the day of the week for 8 May, the earliest day
%  that Mother's day can occur.
%
  m = 5;
  d = 8;
  f = 0.0;

  w = ymdf_to_weekday_common ( y, m, d, f );
%
%  W = 1 means this day is Sunday, and day D is Mother's day.
%  Otherwise, figure out how to increment W to 8 (Sunday again);
%  The same increment makes D the correct day number.
%
  if ( w ~= 1 )
    d = d + 8 - w;
  end

  return
end


text:
function [modelCoupled] = coupleRxnList2Rxn(model, rxnList, rxnC, c, u)
% This function adds coupling constraints to the fluxes `vi` of a given list of reactions
% (`RxnList`).The constraints are proportional to the flux `v` of a specified
% reaction `rxnC`, so that for all reactions in `RxnList` `vi ~ vrxnC`.
% For all reactions, a threshold `u` on flux is set (default value: 0.01).
%
% To add a coupling constraint to a reaction, a coupling vector `c` is
% determined (default value 1000). `c` is multiplied by `vrxnC`, so that for
% all irreversible reactions in `RxnList vi - c * vrxnC <= u`.
%
% For all reversible reactions, the following equation holds true for the
% reverse direction: `vi + c * vrxnC >= u`.
%
% The output is a coupled model (`modelCoupled`), in which for every new
% entry in `modelCoupled.b` a "slack" variable has been added
% to `modelCoupled.mets`.
%
% USAGE:
%
%    [modelCoupled] = coupleRxnList2Rxn(model, rxnList, rxnC, c, u)
%
% INPUTS:
%    model:           model structure
%    rxnList:         array of reaction names
%    rxnC:            reaction that should be coupled with each reaction in the
%                     reaction list
%    c:               vector of coupling factors for each rxn in rxnList (default c = 1000)
%    u:               vector of lower bounds one reaction couples (default u = 0.01)
%
% OUTPUT:
%    modelCoupled:    coupled model
%
% .. Authors:
%    - Sept 2011 AH/IT
%    - May 2012: Added a warning if the coupled reaction is not in model. AH

if ~exist('rxnList', 'var') || isempty(rxnList)
    rxnList = model.rxns;
end

if ischar(rxnC)
    rxnC = {rxnC};
end

if ~exist('c', 'var') || isempty(c)
    c = 1000;
end

if ~exist('u', 'var') || isempty(u)
    u = 0.01;
end

nRxnList = numel(rxnList);
% create the constraint IDs.
ctrs = [strcat('slack_',rxnList)';strcat('slack_',rxnList,'_R')'];
ctrs = ctrs(:);
% get those reactions which are not reversible.
[pres,pos] = ismember(rxnList,model.rxns);
revs = model.lb(pos(pres)) < 0;
toRemove = [false(1,nRxnList);~revs'];
toRemove = toRemove(:);

% if rxnC is not part of the rxnList, we add it to the end for constraint
% addition.
if isempty(intersect(rxnList,rxnC))
    rxnList(end+1) = rxnC;
end
% get the rxnC Position;
rxnCID = find(ismember(rxnList,rxnC));

plusminus = [ones(1,nRxnList);-ones(1,nRxnList)];
plusminus = plusminus(:);

% generate the coefficient matrix
coefs = sparse(2*nRxnList,nRxnList + numel(setdiff(rxnC,rxnList)));
% determine the indices for the coefficients.
rxnInd =  [1:nRxnList;1:nRxnList];
rxnInd = rxnInd(:);
constInd = 1:2*nRxnList;
coefs(sub2ind(size(coefs),constInd',rxnInd)) = 1;
% set the coupling coefficients
cs = - plusminus * c;
coefs(:,rxnCID ) = cs;
% determine the senses for the indices
dsenses = [repmat('L',1,nRxnList);repmat('G',1,nRxnList)];
dsenses = dsenses(:);
ds = plusminus * u;

% remove non reversible reactions
ds = ds(~toRemove);
ctrs = ctrs(~toRemove);
coefs = coefs(~toRemove,:);
dsenses = dsenses(~toRemove);

modelCoupled = addCOBRAConstraints(model,rxnList,ds,'c', coefs,'dsense',dsenses, 'ConstraintID',ctrs);

author: opencobra

repo: cobratoolbox


text:
%%%  e8.4_FLY
clear
%The maximum control speed and maximum navigation angle can be modified as appropriate
path( './icon/',path);
load vehicle_local_position.mat;

%Constant value
RAD2DEG = 57.2957795;
DEG2RAD = 0.0174533;
%throttle when UAV is hovering
THR_HOVER = 0.609;

%% control parameter
%Attitude PID parameters
Kp_PITCH_ANGLE =6.5;
Kp_PITCH_AngleRate = 0.1;
Ki_PITCH_AngleRate = 0.02;
Kd_PITCH_AngleRate = 0.001;
Kp_ROLL_ANGLE =6.5;
Kp_ROLL_AngleRate = 0.1;
Ki_ROLL_AngleRate = 0.02;
Kd_ROLL_AngleRate = 0.001;

Kp_YAW_ANGLE = 3;
Kp_YAW_AngleRate = 0.5;
Ki_YAW_AngleRate = 0.01;
Kd_YAW_AngleRate = 0.00;

%Position PID parameters
Kpxp = 1.0;
Kpyp = 1.0;
Kpzp = 4.0;
Kvxp = 3.0; Kvxi = 0.1; Kvxd = 0.01;
Kvyp = 3.0; Kvyi = 0.1; Kvyd = 0.01;
Kvzp = 0.45; Kvzi = 0.01; Kvzd = 0.005;

Saturation_I_RP_Max = 0.3;
Saturation_I_RP_Min = -0.3;
Saturation_I_Y_Max = 0.2;
Saturation_I_Y_Min = -0.2;
Saturation_I_ah = 3.43;
Saturation_I_az = 5;

%max control angle,default 35deg
MAX_CONTROL_ANGLE_ROLL = 35;
MAX_CONTROL_ANGLE_PITCH  = 35;
%Maximum navigation angle, used in position control to avoid large attitude angles
MAX_CONTROL_ANGLE_NAV_ROLL = 15;
MAX_CONTROL_ANGLE_NAV_PITCH  = 15;
%max control angle rate,rad/s
MAX_CONTROL_ANGLE_RATE_PITCH = 220;
MAX_CONTROL_ANGLE_RATE_ROLL = 220;
MAX_CONTROL_ANGLE_RATE_Y = 200;
%Maximum control speed, m/s
MAX_CONTROL_VELOCITY_XY = 0.5;
MAX_CONTROL_VELOCITY_Z = 1;
%Throttle amplitude
MAX_MAN_THR = 0.9;
MIN_MAN_THR = 0.05;
%% run simulink model
e8_4_FLY

author: RflySim

repo: RflyExpCode

save_path: github-repos/MATLAB/RflySim-RflyExpCode


text:
function [dateV] = doy2date(doyV,yearV)
% DOY2DATE.m will convert a vector of day of year numbers and years
% and convert them to MATLAB date format.
%
% Sample Call:
%  doyV = [54;200.4315];
%  yearV = [2009;2009];
%  [dateV] = doy2date(doyV,yearV);
%
% Inputs:
%  doyV -> vector of day of year numbers (n x 1)
%  yearV -> vector of years (n x 1)
%
% Outputs:
%  dateV -> vector of MATLAB dates (n x 1)
%
% AUTHOR    : A. Booth (ashley [at] boothfamily [dot] com)
% DATE      : 22-May-2009 09:34:53
% Revision  : 1.00
% DEVELOPED : 7.4.0.287 (R2007a) Windows XP
% FILENAME  : doy2date.m

%Check size of incoming vectors
if (size(doyV,1)== 1) | (size(doyV,2) == 1) %Make sure only a nx1 vector
    if size(doyV,1)<size(doyV,2)
        doyV = doyV';
        colflip = 1; %take note if the rows were flipped to col
    else
        colflip = 0;
    end
else %check to see that vectors are columns:
    error('DOY vector can not be a matrix')
end

%year vector
if (size(yearV,1)== 1) | (size(yearV,2) == 1) %Make sure only a nx1 vector
    if size(yearV,1)<size(yearV,2)
        yearV = yearV';
%         colflip = 1; %take note if the rows were flipped to col
    else
%         colflip = 0;
    end
else %check to see that vectors are columns:
    error('Year vector can not be a matrix')
end

%Check to make sure sizes of the vectors are the same
if ~min(size(doyV) == size(yearV))
    error('Day of year vector and year vector must be the same size')
end

%Make year into date vector
z = zeros(length(yearV),5);
dv = horzcat(yearV,z);

%Calc matlab date
dateV = doyV + datenum(dv);

% flip output if input was flipped
if colflip
    dateV = dateV';
end

% disp('Completed doy2date.m')
% ===== EOF [doy2date.m] ======

author: goGPS-Project


text:
function [W,Idx] = barycentricWeights(X,xBnd,N)
% [W,Idx] = barycentricWeights(X,xBnd,N)
%
% This function computes the weights and indicies such that
%
% INPUTS:
%   X = [d , n] = n points of interest, each is d-dimensional
%   xBnd = [d , 2] = [min,max] along each dimension
%   N = [1, d] = number of grid points along each dimension
%
% OUTPUTS:
%   W = [d+1, n] weight to apply to each index
%   Idx = [d+1, n] linear index corresponding to each weight
%
% NOTES:
%   --> y = f(x), where x is d-dimensional, and y is scalar
%   --> Y = Y(x1,x2,...,xd)
%   --> n = size(Y);
%   --> y(X(:,i)) = dot(W(:,i)*Y(Idx(:,i));
%
% ASSUME:
%   --> all data is uniformly spaced:
%       {xi} = linspace(xBnd(i,1), xBnd(i,2), n(i));
%
% REFERENCE:
%
% Based on paper: "Multidimensional Triangulation and Interpolation for
% Reinforcement Learning" by Scott Davies, NIPS 1996
%
% See Also: barycentricInterpolate

[d, n] = size(X);
W = zeros(d+1,n);
Idx = zeros(d+1,n);
for i=1:n
    [W(:,i),Idx(:,i)] = barycentricWeightsCore(X(:,i),xBnd,N);
end

end

function [w,idx] = barycentricWeightsCore(x,xBnd,N)
% [w,idx] = barycentricWeights(x,xBnd,n)
%
% This function computes the weights and indicies such that y = Y(idx)*w,
% where Y is some data, stored on a length(N) dimensional grid.
%
% INPUTS:
%   x = [d , 1] point of interest
%   xBnd = [d , 2] = [min,max] along each dimension
%   N = [1, d] = number of grid points along each dimension
%
% OUTPUTS:
%   w = [d+1,1] weight to apply to each index
%   idx = [d+1,1] linear index corresponding to each weight
%
% NOTES:
%   --> y = f(x), where x is d-dimensional, and y is scalar
%   --> Y = Y(x1,x2,...,xd)
%   --> n = size(Y);
%   --> y = dot(Y(idx),w);
%
% ASSUME:
%   --> all data is uniformly spaced:
%       {xi} = linspace(xBnd(i,1), xBnd(i,2), n(i));
%
% REFERENCE:
%
% Based on paper: "Multidimensional Triangulation and Interpolation for
% Reinforcement Learning" by Scott Davies, NIPS 1996

% Coerce any data that is out of range:
checkLow = x<xBnd(:,1); x(checkLow) = xBnd(checkLow,1);
checkUpp = x>xBnd(:,2); x(checkUpp) = xBnd(checkUpp,2);

% Translate and scale such that each grid cell is a unit cube:
x = x-xBnd(:,1);
x = x./(xBnd(:,2)- xBnd(:,1));
x = x.*(N-1);

% Seperate out into bin number and fractional length:
xBaseIdx = floor(x);  %Figure out the lower index (bin) number
checkCeil = xBaseIdx==(N-1); newCeil = N-2;
xBaseIdx(checkCeil) = newCeil(checkCeil);  %Fix special case for points on upper grid
x = x-xBaseIdx; %Keep the fraction along each dimension:

% Figure out which dimensions to build simplex along:
[~,I] = sort(x); %x must be a column vector!

% Build simplex by walking along each successive dimension
% The matrix X stores the verticies of the simplex, where each vertex is on
% the unit hypercube. The first vertex is always zeros(n,1) and the last is
% always ones(n,1). The intermediate verticies are selected based on the
% sort step. Since X will be triangular, and all entries are unity, it is
% trivial to solve the system: x = X*w for the unknown weights w by row
% reduction.
n = length(N);
X = zeros(n,n+1);
w = zeros(n+1,1);
wSum = 0;
for i=1:n
    X(I(i), (n+2-i):(n+1)) = 1;
    w(n+2-i) = x(I(i)) - wSum;
    wSum = wSum + w(n+2-i);
end
w(1) = 1-wSum;  %Final convex combination

% Compute the linear indices:
idx = zeros(n+1,1);
for i=1:(n+1)
    idx(i) = sub2idx(N, X(:,i)+xBaseIdx+1, n);
end

end

%%%% Custom version of sub2ind, that takes a vector input for the desired
%%%% input dimension, rather than a comma seperated list.
function idx = sub2idx(n,idxVec,d)
%
% Based on Matlab's sub2ind command
% Modified by Matthew Kelly, April 8, 2015
%
% CHANGES:
%   - must explicitly enter all dimensions: length(n) == length(idxVec)
%   - pass through for scalar: length(n)==1 is valid input
%   - must pass desired subscripts as a vector (not a comma-seperated-list)
%   - disabled most error checking
%
% INPUTS:
%   n = size(dataMatrix)
%   idxVec = list of desired indices. For example A(2,1) -> [2,1]
%
% OUTPUTS:
%   idx = linear index, such that:
%
%       A(idx) = A(idxVec(1), ... , idxVec(length(n));
%

if d == 1
    idx = idxVec;
elseif d ==2
    idx = idxVec(1) + (idxVec(2) - 1).*n(1);
else
    %Compute linear indices
    k = [1 cumprod(n(1:end-1))'];
    idx = 1;
    for i = 1:d
        v = idxVec(i);
        idx = idx + (v-1)*k(i);
    end
end

end

author: MatthewPeterKelly

repo: dscTutorials


text:
%bilinearFit
% Fits data to bilinear function
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Affiliation: Covert Lab, Department of Bioengineering, Stanford University
% Last Updated: 1/11/2012
function [a, b, cutX, gof] = bilinearFit(x, y)

tf = ~isnan(x) & ~isnan(y);
x = x(tf);
y = y(tf);

[~, order] = sort(x);
x = x(order);
y = y(order);

a = NaN(numel(x), 2);
b = NaN(numel(x), 2);
sse = NaN(numel(x), 2);
rsquare = NaN(numel(x), 2);
adjrsquare = NaN(numel(x), 2);
for i = 2:numel(x)-2
    [fitobject1, gof1] = fit(x(1:i), y(1:i), 'poly1');
    [fitobject2, gof2] = fit(x(i+1:end), y(i+1:end), 'poly1');
    a(i, :) = [fitobject1.p2 fitobject2.p2];
    b(i, :) = [fitobject1.p1 fitobject2.p1];
    sse(i, :) = [gof1.sse gof2.sse];
    rsquare(i, :) = [gof1.rsquare gof2.rsquare];
    adjrsquare(i, :) = [gof1.adjrsquare gof2.adjrsquare];
end

[~, cutIdx] = nanmin(sum(sse, 2));
cutX = mean(x(cutIdx:cutIdx+1));

a = a(cutIdx, :);
b = b(cutIdx, :);
sse = sse(cutIdx, :);
rsquare = rsquare(cutIdx, :);
adjrsquare = adjrsquare(cutIdx, :);

sse0 = sum(sse);
rsquare0 = 1 - sse0 / sum((y-mean(y)).^2);
n = numel(x);
p = 3;
adjrsquare0 = 1 - (1 - rsquare0) * (n-1) / (n-p-1);

gof = struct(...
    'sse', [sse0 sse], ...
    'rsquare', [rsquare0 rsquare], ...
    'adjrsquare', [adjrsquare0 adjrsquare]);

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function [Z,order,Md] = som_linkage(sM,varargin)

%SOM_LINKAGE Make a hierarchical linkage of the SOM map units.
%
% [Z,order,Dist] = som_linkage(sM, [[argID,] value, ...])
%
%  Z = som_linkage(sM);
%  Z = som_linkage(D,'complete');
%  Z = som_linkage(sM,'single','ignore',find(~som_hits(sM,D)));
%  Z = som_linkage(sM,pdist(sM.codebook,'mahal'));
%  som_dendrogram(Z);
%
%  Input and output arguments ([]'s are optional):
%   sM       (struct) map or data struct to be clustered
%            (matrix) size dlen x dim, a data set: the matrix must not
%                     contain any NaN's!
%   [argID,  (string) See below. The values which are unambiguous can
%    value]  (varies) be given without the preceeding argID.
%
%   Z        (matrix) size dlen-1 x 3, the linkage info
%                     Z(i,1) and Z(i,2) hold the indeces of clusters
%                     combined on level i (starting from bottom). The new
%                     cluster has index dlen+i. The initial cluster
%                     index of each unit is its linear index in the
%                     original data matrix. Z(i,3) is the distance
%                     between the combined clusters. See LINKAGE
%                     function in the Statistics Toolbox.
%                     The ignored samples are listed at the
%                     end of the Z-matrix and have Z(*,3) == Inf
%   Dist     (matrix) size dlen x dlen, pairwise distance matrix
%
% Here are the valid argument IDs and corresponding values. The values
% which are unambiguous (marked with '*') can be given without the
% preceeding argID.
%   'linkage' *(string) the linkage criteria to use: 'single' (the
%                       default), 'average' or 'complete'
%   'topol'   *(struct) topology struct
%   'connect' *(string) 'neighbors' or 'any' (default), whether the
%                       connections should be allowed only between
%                       neighbors or between any vectors
%              (matrix) size dlen x dlen indicating the connections
%                       between vectors
%              (scalar) the N-neighborhood upto which the connections
%                       should be formed (implies 'neighbors')
%   'ignore'   (vector) the units/vectors which should be ignored
%   'dist'     (matrix) size dlen x dlen, pairwise distance matrix to
%                       be used instead of euclidian distances
%              (vector) as the output of PDIST function
%              (scalar) distance norm to use (euclidian = 2)
%   'mask'     (vector) size dim x 1, the search mask used to
%                       weight distance calculation. By default
%                       sM.mask or a vector of ones is used.
%
% Note that together 'connect'='neighbors' and 'ignore' may form
% areas on the map which will never be connected: connections
% across the ignored map units simply do not exist.
%
% See also KMEANS_CLUSTERS, LINKAGE, PDIST, DENDROGRAM.

% Copyright (c) 2000 by Juha Vesanto
% Contributed to SOM Toolbox on June 16th, 2000 by Juha Vesanto
% http://www.cis.hut.fi/projects/somtoolbox/

% Version 2.0beta juuso 160600

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% input arguments

% the data
if isstruct(sM),
  if isfield(sM,'data'), D = sM.data; sT = []; mask = [];
  else D = sM.codebook; sT = sM.topol; mask = sM.mask;
  end
else
  D = sM; sT = []; mask = [];
end
[dlen dim] = size(D);
if isempty(mask), mask = ones(dim,1); end
if any(isnan(D(:))), error('Data matrix must not have any NaNs.'); end

% varargin
Md = 2;
linkage = 'single';
ignore_units = [];
constrained = 0;
i=1;
while i<=length(varargin),
  argok = 1;
  if ischar(varargin{i}),
    switch varargin{i},
      % argument IDs
     case {'topol','som_topol','sTopol'}, i=i+1; sT = varargin{i};
     case 'connect', i=i+1;
       if ischar(varargin{i}), constrained = ~strcmp(varargin{i},'any');
       else constrained = varargin{i}; end
     case 'ignore',  i=i+1; ignore_units = varargin{i};
     case 'dist',    i=i+1; Md = varargin{i};
     case 'linkage', i=i+1; linkage = varargin{i};
     case 'mask',    i=i+1; mask = varargin{i};
     case 'tracking',i=i+1; tracking = varargin{i};
      % unambiguous values
     case 'neighbors', constrained = 1;
     case 'any',       constrained = 0;
     case {'single','average','complete'}, linkage = varargin{i};
     otherwise argok=0;
    end
  elseif isstruct(varargin{i}) & isfield(varargin{i},'type'),
    switch varargin{i}(1).type,
     case 'som_topol', sTopol = varargin{i};
     otherwise argok=0;
    end
  else
    argok = 0;
  end
  if ~argok,
    disp(['(som_linkage) Ignoring invalid argument #' num2str(i+1)]);
  end
  i = i+1;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% distance matrix

% given distance matrix % jh corrected this place totally 27.3. 03
if (prod(size(Md))==1), % no explicit distance matrix, set flag
  q=2; % 17.2.03 kr added some brackets
else
  if (prod(size(Md))<dlen^2), % check pdist form
    Md = squareform(Md);   % transform to ordinary square diastance matrix
  end
  % jh: 27.3. 03 "calculate pairwise dist. matrix" see approx. 20 lines below
  % sets self-distance to Inf! This must be set here also,
  % otherwise clustering fails for user defined dist. matrix!
  Md(eye(dlen)==1)=Inf;
end

% neighborhood constraint
if length(constrained)==1 & constrained>0,
  Ne1 = som_unit_neighs(sT);
  Conn = som_neighborhood(Ne1,constrained);
  Conn(~isfinite(Conn(:))) = 0;
else Conn = constrained; end
if ~isempty(Conn), for i=1:dlen, C(i,i) = 1; end, end

% pairwise distance matrix across connected units
n = size(D,1);
if prod(size(Md))>1,
  % remove distances between non-neighbors
  if constrained, for i = 1:n, Md(i,find(Conn(i,:)==0)) = Inf; end, end
else
  % calculate pairwise distance matrix
  q = Md;
  Md = zeros(n,n)+Inf;
  if ~constrained & q==2, % fast for the default case
    for i = 1:n-1,
      x = D(i,:);
      inds = [(i+1):n];
      Diff = D(inds,:) - x(ones(n-i,1),:);
      Md(inds,i) = sqrt((Diff.^2)*mask);
      Md(i,inds) = Md(inds,i)';
    end
  else
    for i = 1:n-1,
      inds = find(Conn(i,:)==1);
      inds = inds(find(inds>i));
      Diff = abs(D(inds,:) - D(i*ones(length(inds),1),:));
      switch q,
       case 1,    dist = Diff*mask;
       case 2,    dist = sqrt((Diff.^2)*mask);
       case Inf,  dist = max(Diff,[],2);
       otherwise, dist = ((Diff.^q)*mask).^(1/q);
      end
      Md(inds,i) = dist;
      Md(i,inds) = dist';
    end
  end
end

% set distances to ignored units to Inf
if ~isempty(ignore_units),
  Md(ignore_units,:) = Inf;
  Md(:,ignore_units) = Inf;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% construct dendrogram

Z = zeros(n-1,3)+NaN;  % merged clusters and distance for each step
clusters = 1:dlen;     % each vector is at first in its own cluster
Cd = Md;               % distances between clusters

h = waitbar(0,'Constructing hierarchical clustering');

for i=1:n-1,

  % tracking
  waitbar(i/(n-1),h);

  %% combine two closest clusters
  % find the clusters which are closest to each other (c1 and c2)
  [d,ind] = min(min(Cd));
  if ~isfinite(d), break; end   % no more connected clusters
  [d,c1] = min(Cd(:,ind));      % cluster1
  c2 = clusters(ind);           % cluster2
  % combine the two clusters
  c1_inds = find(clusters==c1); % vectors belonging to c1
  c2_inds = find(clusters==c2); % vectors belonging to c2
  c_inds = [c1_inds, c2_inds];  % members of the new cluster
  % new cluster index = bigger cluster
  if length(c2_inds)>length(c1_inds), c=c2; k=c1; else c=c1; k=c2; end
  clusters(c_inds) = c;         % update cluster info
  Z(i,:) = [c, k, d];           % save info into Z

  %% update cluster distances
  % remove the subclusters from the Cd table
  Cd(c_inds,c_inds) = Inf;      % distance of the cluster to its members = Inf
  k_inds = c_inds(c_inds ~= c); % vectors of the smaller cluster
  Cd(k_inds,:) = Inf;           % distance of the subclusters to
  Cd(:,k_inds) = Inf;           % other clusters = Inf
  % update the distance of this cluster to the other clusters
  cl = unique(clusters(clusters ~= c)); % indeces of all other clusters
  if ~isempty(cl), % added since v6.5 works differently than 6.1
    for l=cl,
      o_inds = find(clusters==l); % indeces belonging to cluster k
      vd = Md(c_inds,o_inds);     % distances between vectors in c and k
      vd = vd(isfinite(vd(:)));   % remove infinite distances (no connection)
      len = length(vd);
      if ~len, % if the two clusters are not connected, their distance in Inf
	sd = Inf;
      else   % otherwise, calculate the distance between them
	switch linkage,
	 case 'single',   sd = min(vd);
	 case 'average',  sd = sum(vd)/len;
	 case 'complete', sd = max(vd);
	 otherwise, error(['Unknown set distance: ' linkage]);
	end
      end
      Cd(c,l) = sd; Cd(l,c) = sd;
    end
  end
end
close(h);

last = Z(i,1);
if isnan(last),
  last = Z(i-1,1);
  rest = setdiff(unique(clusters),last);
  Z(i:n-1,1) = rest';
  Z(i:n-1,2) = last;
  Z(i:n-1,3) = Inf;
  i = i - 1;
else
  rest = [];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% return values

% calculate the order of samples
order = last;
% go through the combination matrix from top to down
for k=i:-1:1,
  c = Z(k,1); k = Z(k,2);                     % what (k) change to what (c)
  j = find(order==c);                         % the occurance of c in order
  if j == length(order), order = [order k];   % put k behind c
  else order = [order(1:j) k order(j+1:end)];
  end
end
order = [rest, order];

% to maintain compatibility with Statistics Toolbox, the values in
% Z must be yet transformed so that they are similar to the output
% of LINKAGE function

Zs = Z;
current_cluster = 1:dlen;
for i=1:size(Z,1),
  Zs(i,1) = current_cluster(Z(i,1));
  Zs(i,2) = current_cluster(Z(i,2));
  current_cluster(Z(i,[1 2])) = dlen+i;
end

Z = Zs;

return;

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

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function cvt_3d_sampling ( g_num, it_num, s_num )

%*****************************************************************************80
%
%% CVT_3D_SAMPLING carries out the Lloyd algorithm in a 3D unit box.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 June 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer G_NUM, the number of generators.
%    A value of 50 is reasonable.
%
%    Input, integer IT_NUM, the number of CVT iterations.
%    A value of 20 or 50 might be reasonable.
%
%    Input, integer S_NUM, the number of sample points to use
%    when estimating the Voronoi regions.
%    A value of 1,000 is too low.  A value of 1,000,000 is somewhat high.
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CVT_3D_SAMPLING\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Use sampling to approximate Lloyd''s algorithm\n' );
  fprintf ( 1, '  in the 3D unit cube.\n' );

  if ( nargin < 1 )
    g_num = input ( '  Enter number of generators, G_NUM: ' );
  elseif ( ischar ( g_num ) )
    g_num = str2num ( g_num );
  end

  if ( nargin < 2 )
    it_num = input ( '  Enter number of iterations, IT_NUM: ' );
  elseif ( ischar ( it_num ) )
    it_num = str2num ( it_num );
  end

  if ( nargin < 3 )
    s_num = input ( '  Enter number of sample points, S_NUM: ' );
  elseif ( ischar ( s_num ) )
    s_num = str2num ( s_num );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of generators is %d\n', g_num );
  fprintf ( 1, '  Number of iterations is %d\n', it_num );
  fprintf ( 1, '  Number of samples is %d\n', s_num );
%
%  Initialize the generators.
%  Choice 1: RAND
%
  g = rand ( g_num, 3 );
%
%  Carry out the iteration.
%
  step = 1 : it_num;
  e = nan ( it_num, 1 );
  gm = nan ( it_num, 1 );

  for it = 1 : it_num

    figure ( 1 )
    scatter3 ( g(:,1), g(:,2), g(:,3), 'filled' )
    xlabel ( '<---X--->' )
    ylabel ( '<---Y--->' )
    zlabel ( '<---Z--->' )
    title ( sprintf ( 'Generator locations before iteration %d', it ) );
%
%  Compute the Delaunay triangle information T for the current nodes.
%
    t = DelaunayTri ( g );
%
%  Generate sample points.
%
    s = rand ( s_num, 3 );
%
%  For each sample point, find K, the index of the nearest generator.
%  We do this efficiently by using the Delaunay information with
%  Matlab's DSEARCH command, rather than a brute force nearest neighbor
%  computation.
%
    k = nearestNeighbor ( t, s );

    m(:,1) = accumarray ( k, ones(s_num,1) );

    g_new(:,1) = accumarray ( k, s(:,1) ) ./ m(:,1);
    g_new(:,2) = accumarray ( k, s(:,2) ) ./ m(:,1);
    g_new(:,3) = accumarray ( k, s(:,3) ) ./ m(:,1);
%
%  Compute the average energy.
%
    e(it,1) = sum ( ( s(:,1) - g(k(:,1),1) ).^2 ...
                  + ( s(:,2) - g(k(:,1),2) ).^2 ...
                  + ( s(:,3) - g(k(:,1),3) ).^2 ) / s_num;
%
%  Display the energy.
%
    figure ( 3 )
    plot ( step, log ( e ), 'm-*' )
    title ( 'Log (Energy)' )
    xlabel ( 'Step' )
    ylabel ( 'Energy' )
    grid
%
%  Compute the generator motion.
%
    gm(it,1) = sum ( ( g_new(:,1) - g(:,1) ).^2 ...
                   + ( g_new(:,2) - g(:,2) ).^2 ...
                   + ( g_new(:,3) - g(:,3) ).^2 ) / g_num;
%
%  Display the generator motion.
%
    figure ( 4 )
    plot ( step, log ( gm ), 'm-*' )
    title ( 'Log (Average generator motion)' )
    xlabel ( 'Step' )
    ylabel ( 'Motion' )
    grid
%
%  Continue?
%
    s = input ( 'RETURN, or Q to quit: ', 's' );

    if ( s == 'q' | s == 'Q' )
      break
    end
%
%  Update the generators.
%
    g = g_new;

  end
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CVT_3D_SAMPLING\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function timestamp ( )

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

  return
end


text:
function a = c8vec_nint ( n, a )

%*****************************************************************************80
%
%% C8VEC_NINT rounds the entries of a C8VEC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vector.
%
%    Input, complex A(N), the vector to be rounded.
%
%    Output, complex A(N), the rounded vector.
%
  for i = 1 : n
    a(i) = c8_nint ( a(i) );
  end

  return
end


text:
function templates_test ( )

%*****************************************************************************80
%
%% TEMPLATES_TEST tests the TEMPLATES library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 March 2006
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEMPLATES_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the TEMPLATES library.\n' );

  bicg_test ( );
  bicgstab_test ( );
  cg_test ( );
  cgs_test ( );
  cheby_test ( );
  gmres_test ( );
  gs_test ( );
  jacobi_test ( );
  lehmer_test ( );
  mm_to_msm_test ( );
  poisson_test ( );
  qmr_test;
  sor_test ( );
  wathen_test ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEMPLATES_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function F=dwt_fuse(I1,I2,zt)
%    DWT
%    Input:
%    I1 - input image A
%    I2 - input image B
%    zt - maximum decomposition level
%    Output:
%    F  - fused image
%
%    The code is edited by Yu Liu, 01-09-2014.

%-------------------------------------------------------------------------%
%                                DWT
%-------------------------------------------------------------------------%
I1=double(I1);
I2=double(I2);
tempA=I1;
tempB=I2;

X=cell(zt,4);
Y=cell(zt,4);
Z=cell(zt,4);
for i=1:zt
    [X{i,1},X{i,2},X{i,3},X{i,4}]=dwt2(tempA,'db1','mode','per');
    tempA=X{i,1};
    [Y{i,1},Y{i,2},Y{i,3},Y{i,4}]=dwt2(tempB,'db1','mode','per');
    tempB=Y{i,1};
end

%-------------------------------------------------------------------------%
%                               low-pass fusion
%-------------------------------------------------------------------------%
Z{zt,1}=(X{zt,1}+Y{zt,1})/2;

%-------------------------------------------------------------------------%
%                               high-pass fusion
%-------------------------------------------------------------------------%
for i=zt:-1:1
    for j=2:4
        Z{i,j}=selc(X{i,j},Y{i,j},3);
    end
end

%-------------------------------------------------------------------------%
%                               IDWT
%-------------------------------------------------------------------------%
for i=zt:-1:1
    if i>1
        Z{i-1,1}=idwt2(Z{i,1},Z{i,2},Z{i,3},Z{i,4},'db1','mode','per');
    else
        F=idwt2(Z{i,1},Z{i,2},Z{i,3},Z{i,4},'db1','mode','per');
    end
end

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function s_out = structmerge(s1, s2, varargin)
% STRUCTMERGE - Merge struct variables into a common struct
%
%  S_OUT = STRUCTMERGE(S1, S2, ...)
%
%  Merges the structures S1 and S2 into a common structure S_OUT
%  such that S_OUT has all of the fields of S1 and S2. When
%  S1 and S2 share the same fieldname, the value of S2 is taken.
%
%  The behavior of the function can be altered by passing additional
%  arguments as name/value pairs.
%
%  Parameter (default)     | Description
%  ------------------------------------------------------------
%  ErrorIfNewField (0)     | (0/1) Is it an error if S2 contains a
%                          |  field that is not present in S1?
%
%  See also: STRUCT

ErrorIfNewField = 0;

assign(varargin{:});

f1 = fieldnames(s1);
f2 = fieldnames(s2);

if ErrorIfNewField,
	[c,f1i] = setdiff(f2,f1);
	if ~isempty(c),
		error(['Some fields of the second structure are not in the first: ' cell2str(c) '.']);
	end;
end;

s_out = s1;

for i=1:length(f2),
	s_out = setfield(s_out,f2{i},getfield(s2,f2{i}));
end;

author: flatironinstitute

repo: CaImAn-MATLAB


text:
% simple file for demonstrating the Feldkamp - Devis - Kress algorithm
% for tomo reconstruction according to cone beam geometry

clear all,
close all

global size_proj % size of the projection
global Sinog % make available projections to read_1proj function

%
NI=128; % NI is a scalar that specifies the number of rows and col. in the phantom.
phant = phantom(NI); % generating phantom image

% Generate set of projections around 360 degrees with a fan beam geometry
Radius=150; % distance focal spot - object ratation axis
[F,Floc,Fangles] = fanbeam(phant,150,'FanRotationIncrement',1, 'FanSensorGeometry' ,'line');
Sinog=rot90(F,-1); % sinogram
%figure, imshow(Floc,Fangles,Sinog,[ ],'n'), axis normal
figure, imshow(Sinog,[ ]), axis normal
ylabel('Rotation Angles (degrees)')
xlabel('Sensor Positions (degrees)')
title('sinogram');
colormap(hot), colorbar

% Recontruction from sinthetic projections
tic
Ifb = ifanbeam(F,Radius,'FanSensorGeometry' ,'line');
toc
figure, imshow(phant,[]); title('phantom');
figure, imshow(Ifb); title(' ifanbeam');

% Using FDK cone beam code to reconstruct the central slice

% filter parameters
df=1;
filter='ram-lak' % high pass for sintetic images
% interpolation parameter : notice that we use bilinear in cone beam
%and linear in fan beam
interp='nearest neighbor',  % interpolation method
interp='bilinear',  % interpolation method
interp='n_n_demo',  % interpolation method

% not needed for this demo
dir='dummy' % radiographic directoryname
fl_prefix='pcon_'; % prefix tomographic data i
%

range_slices=[ ] ; % only central slice is reconst
step=1; % all slices
Proj_Crop=[1 size(Sinog,2) ]; % i.e. no crop because the sinogram in centered

% call iconebeam function for filterd backprojection
NI=0;
tic

[Vol] = ...
    iconebeamdemo2007(dir,fl_prefix,-Fangles,interp,filter,df,NI,step,Radius,range_slices,Proj_Crop);

toc

[svx,svy,svz]=size(Vol);
for is=1:svz
    vct=Vol(:,:,is); vct=vct(:);
%b=max(vct);
%a=min(vct);
a=prctile(Vol(:),3);
b=prctile(Vol(:),97);
end

if svz >1
% Rendering !!
figure,
hs=slice(Vol,[round(svx/2)],[round(svy/2)],[round(svz/3), round(svz*2/3)]);%,'linear') ; %SLICE(V,Sx,Sy,Sz)
view(3); %grid off;
set(hs,'FaceColor','interp','Edgecolor','none')
colormap(gray(256));
%draws slices along the x,y,z directions at the points in the vectors Sx,Sy,Sz.
else
   Img=Vol(:,:,1)'; Img=fliplr(Img);
   figure, imshow(Img,[a  b]) ; title(' iconebeam');

end

%save filename Vol % to save vol on disk !
% author: schena@units.it


text:
function y=remove_CP(x,Ncp,Noff)
% Remove cyclic prefix
if nargin<3
    Noff=0;
end
y=x(:,Ncp+1-Noff:end-Noff);

author: LyricYang


text:
function res = sensors(this, type, newsens)
% Sets and gets sensor fields for EEG and MEG
% returns empty matrix if no sensors are defined.
% FORMAT res = sensors(this, type, newsens)
%   type - 'EEG' or 'MEG'
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: sensors.m 4790 2012-07-17 16:00:20Z vladimir $

if nargin<2
    error('Sensor type (EEG or MEG) must be specified');
end

switch lower(type)
    case 'eeg'
        if nargin < 3
            if isfield(this.sensors, 'eeg')
                res = this.sensors.eeg;
            else
                res = [];
            end
        else
            this.sensors(1).eeg = ft_datatype_sens(newsens);
            res = this;
        end
    case 'meg'
        if nargin < 3
            if isfield(this.sensors, 'meg')
                res = this.sensors.meg;
            else
                res = [];
            end
        else
            this.sensors(1).meg = ft_datatype_sens(newsens);
            res = this;
        end
    otherwise
        error('Unsupported sensor type');
end


text:
function b = r8vec_sht ( n, a  )

%*****************************************************************************80
%
%% R8VEC_SHT computes a "slow" Hartley transform of an R8VEC.
%
%  Discussion:
%
%    The discrete Hartley transform B of a set of data A is
%
%      B(I) = 1/sqrt(N) * Sum ( 0 <= J <= N-1 ) A(J) * CAS(2*PI*I*J/N)
%
%    Here, the data and coefficients are indexed from 0 to N-1.
%
%    With the above normalization factor of 1/sqrt(N), the Hartley
%    transform is its own inverse.
%
%    This routine is provided for illustration and testing.  It is inefficient
%    relative to optimized routines.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Ralph Hartley,
%    A More Symmetrical Fourier Analysis Applied to Transmission Problems,
%    Proceedings of the Institute of Radio Engineers,
%    Volume 30, pages 144-150, 1942.
%
%  Parameters:
%
%    Input, integer ( kind = 4 ) N, the number of data values.
%
%    Input, real A(N), the data to be transformed.
%
%    Output, real B(N), the transformed data.
%
  b(1:n) = 0.0;

  for i = 1 : n
    for j = 1 : n
      theta = 2.0 * pi * mod ( ( i - 1 ) * ( j - 1 ), n ) / n;
      b(i) = b(i) + a(j) * ( cos ( theta ) + sin ( theta ) );
    end
  end

  b(1:n) = b(1:n) / sqrt ( n );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% [HEIGHT] = spyrHt(INDICES)
%
% Compute height of steerable pyramid with given index matrix.

% Eero Simoncelli, 6/96.

function [ht] =  spyrHt(pind)

nbands = spyrNumBands(pind);

% Don't count lowpass, or highpass residual bands
if (size(pind,1) > 2)
  ht = (size(pind,1)-2)/nbands;
else
  ht = 0;
end

repo: SelfExSR


text:
function [h_1d, h_2d] = create_blur_kernel(Kx,opt)

FIX_HSIZE = 1; % hsize = 11
HSIZE = 15;
eps = 0.01;

if FIX_HSIZE
    V = [Kx(1:floor(HSIZE/2)+1); Kx(opt.N-floor(HSIZE/2)+1:opt.N)];
    V = circshift(V, floor(HSIZE/2));
else
    for i = 1:length(Kx)
        if i>2 && Kx(i) < eps
            halflen = i-2;
            break
        end
    end

    len = 2*halflen + 1;
    V = zeros(1, len);
    for i = 1:halflen
        V(i) = Kx(opt.N - halflen + i);
    end
    for i = halflen+1:len
        V(i) = Kx(i-halflen);
    end
end

V = V / sum(V);

h_1d = V';
h_2d = h_1d'*h_1d;

author: thomas-koehler

repo: SupER


text:
function output_info=cnn_layer_dagnn_wrapper_forward(input_info, layer, work_info_batch)

run_trn=work_info_batch.ref.run_trn;
group_info=get_current_work_group_idx(work_info_batch);

net=group_info.dag_net;

if input_info.is_group_data
    input_var_names=layer.input_var_names;
    input_num=length(input_var_names);
    input_data=cell(0);
    for v_idx=1:input_num
        one_x=input_info.data_child_groups{v_idx}.x;
        assert(~isempty(one_x));
        input_data=cat(2, input_data, {input_var_names{v_idx}, one_x});
        assert(~isempty(one_x));
    end
else
    assert(length(layer.input_var_names)==1);
    input_var_name=layer.input_var_names{1};
    input_data={input_var_name, input_info.x};
    assert(~isempty(input_info.x));
end

do_bp_run=group_info.net_info.ref.do_bp && run_trn;

if do_bp_run

    net.mode = 'normal' ;
    net.do_forward_trn(input_data);

else

    net.mode = 'test' ;
    net.eval(input_data) ;
end

output_var_names=layer.output_var_names;
output_num=length(output_var_names);
output_var_idxes=zeros(output_num,1);
for o_idx=1:output_num
    output_var_idxes(o_idx)=net.getVarIndex(output_var_names{o_idx});
end
assert(all(output_var_idxes>0));

if layer.use_single_output
    output_x = net.vars(output_var_idxes(end)).value ;

    output_info=[];
    output_info.is_group_data=false;
    output_info.x=output_x;

    assert(~isempty(output_info.x));
else

    output_num=length(output_var_idxes);
    data_child_groups=cell(output_num, 1);
    for o_idx=1:output_num
        data_child_groups{o_idx}.is_group_data=false;
        data_child_groups{o_idx}.x = net.vars(output_var_idxes(o_idx)).value;

        assert(~isempty(data_child_groups{o_idx}.x));
    end

    output_info=[];
    output_info.is_group_data=true;
    output_info.data_child_groups=data_child_groups;
end

output_info=my_init_input_info(output_info);

% do the print here, generate the network graph :

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

% for cascaded refinenets
% input_vars=net.getInputs;
% net.print({input_vars{1}, [13 13 2048], input_vars{2}, [25 25 1024], input_vars{3}, [50 50 512], input_vars{4}, [100 100 256]}, 'Format', 'dot');

end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function mkdirKPM(dname)
% function mkdirKPM(dname)
% If directory does not exist, make it
%
% mkdirKPM('foo\bar\baz') makes subdirectories bar and baz
% mkdirKPM('foo\bar\baz.txt') makes subdirectories bar

% convert foo\bar\baz to pathstr=foo\bar, name=baz
% convert foo\bar\baz.txt to pathstr=foo\bar, name=baz
[pathstr, name, ext, versn] = fileparts(dname);
if ~isempty(ext) % we stripped off something after final period
  % convert foo\bar to pathstr=foo, name=bar
  % convert foo\bar.bad to pathstr=foo, name=bar.bad
  [pathstr, name, ext, versn] = fileparts(pathstr);
  name = sprintf('%s%s', name, ext); % in case there is a period in the directory name
end

dname = fullfile(pathstr, name);
if ~exist(dname, 'dir')
  %fprintf('mkdirKPM: making %s, %s \n', pathstr, name);
  mkdir(pathstr, name)
else
  %fprintf('mkdirKPM: %s already exists\n', dname)
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
classdef tt_function
    %TT_FUNCTION Summary of this class goes here
    %   Detailed explanation goes here

    properties
        fun;
        coeff;
    end

    methods
        %-----------------------------------------------------------------%
        function obj = tt_function(fun, coeff)
            obj.fun = fun;
            obj.coeff = coeff;
        end

        %-----------------------------------------------------------------%
        function tt = to_tt_tensor(obj, x)
        %TO_TT_TENSOR Converts tt_function to tt_tensor
        %   x - d x 1 cell array of coordinates for evaluation

        d = obj.coeff.d;
        r = obj.coeff.r;
        n = obj.coeff.n;
        ps = obj.coeff.ps;
        cr = obj.coeff.core;
        n_tt = cellfun(@numel, x);
        ps_tt = cumsum([1 ; n_tt.*r(1:d).*r(2:d+1)]);
        cr_tt = zeros(ps_tt(d+1)-1, 1);

        % compute basis functions
        basis = tt_function.get_basis(x, obj.fun, n);

        % evaluate cores
        for i = 1: d
            c = reshape(cr(ps(i):ps(i+1)-1), [r(i) n(i) r(i+1)]);
            t = ten_conv(c, 2, basis{i});
            cr_tt(ps_tt(i):ps_tt(i+1)-1) = reshape(t, [r(i)*n_tt(i)*r(i+1) 1]);
        end

        tt = tt_tensor;
        tt.d = d;
        tt.r = r;
        tt.n = n_tt;
        tt.ps = ps_tt;
        tt.core = cr_tt;

        end

        %-----------------------------------------------------------------%
        function val = eval(obj, x)
        %EVAL Evaluates tt_function at a given d-dimensional point
        %   Inefficient code, simple solution

        d = numel(x);
        t = cell(d,1);
        for i = 1: d
            t{i} = x(i);
        end
        tt = obj.to_tt_tensor(t);
        val = tt(ones(1, d));

        end

    end

    methods (Static)
        %-----------------------------------------------------------------%
        function basis = get_basis(x, fun, n)
        %GET_BASIS Constructs a cell array of values of basis functions.

        d = numel(n);
        basis = cell(d, 1);

        for i = 1: d
            basis{i} = zeros(n(i), numel(x{i}));
            for j = 1: n(i)
                for k = 1: numel(x{i})
                    basis{i}(j,k) = fun(i, j, x{i}(k));
                end
            end
        end

        end

    end

end

author: oseledets

repo: TT-Toolbox


text:
classdef SynchronizationFile < handle
    properties (Access = public)
        synchronizationPointsMap;
        fileName;
        count;%dependent property
    end

    methods
        function c = get.count(obj)
            c = obj.synchronizationPointsMap.Count;
        end
    end

    methods (Access = public)
        function obj = SynchronizationFile(samples,frames)
            obj.synchronizationPointsMap = containers.Map('KeyType','uint64','ValueType','uint32');
            if nargin > 0
                obj.synchronizationPointsMap.keys = samples;
                obj.synchronizationPointsMap.values = frames;
            end
        end

        function isValid = isValidSynchronizationFile(obj)
            isValid = (length(obj.synchronizationPointsMap.keys) >= 2);
        end

        function setSynchronizationPoint(obj,sample,frame)
            obj.synchronizationPointsMap(uint64(sample)) = uint32(frame);
        end

        function removeSynchronizationPointWithSample(obj,sample)
            obj.synchronizationPointsMap.remove(sample);
        end

        function videoFrame = sampleToVideoFrame(obj, sample)
            if ~obj.isValidSynchronizationFile()
                videoFrame = [];
            else
                [sample1, frame1] = obj.getFirstSynchronizationPoint();
                [sample2, frame2] = obj.getLastSynchronizationPoint();

                a = double(frame2 - frame1) / double(sample2 - sample1);
                videoFrame = a * (sample - double(sample1)) + double(frame1);

                if videoFrame < 1
                    videoFrame = 1;
                end

                videoFrame = uint32(videoFrame);
            end
        end

        function sample = videoFrameToSample(obj, videoFrame)
            if ~obj.isValidSynchronizationFile()
                sample = [];
            else
                [sample1, frame1] = obj.getFirstSynchronizationPoint();
                [sample2, frame2] = obj.getLastSynchronizationPoint();

                a = (double(sample2) - double(sample1)) / (double(frame2) - double(frame1));
                sample = a * (double(videoFrame) - double(frame1)) + double(sample1);

                if sample < 1
                    sample = 1;
                end

                sample = uint32(sample);
            end
        end
    end

    methods (Access = private)

        function [sample, frame] = getFirstSynchronizationPoint(obj)
            samples = obj.synchronizationPointsMap.keys;
            frames = obj.synchronizationPointsMap.values;
            if isempty(samples)
                sample = [];
                frame = [];
            else
                [~, minIdx] = min(cell2mat(samples));
                sample = samples{minIdx};
                frame = frames{minIdx};
            end
        end

        function [sample, frame] = getLastSynchronizationPoint(obj)
            samples = obj.synchronizationPointsMap.keys;
            frames = obj.synchronizationPointsMap.values;
            if isempty(samples)
                sample = [];
                frame = [];
            else
                [~, maxIdx] = max(cell2mat(samples));
                sample = samples{maxIdx};
                frame = frames{maxIdx};
            end
        end
    end
end

author: avenix

repo: WDK

save_path: github-repos/MATLAB/avenix-WDK


text:
function patches = img2patch_index(img, patchSize, stride)

m = size(img, 1); n = size(img, 2);
a = patchSize(1); b = patchSize(2); c = stride;

if stride > 1
	aCut = mod(m, a); bCut = mod(n, b);
	img = img(1:end-aCut, 1:end-bCut);
end

% Extended image from the original image
imgExt = [ img img(:, 1:b-c); img(1:a-c, :) img(1:a-c, 1:b-c) ];

if stride == 1
	patches = im2col(imgExt, patchSize, 'sliding');
return
end

iMat = zeros(size(img), class(img));
iMat( 1:c:end, 1:c:end ) = 1; % Take patches in distances of 'stride'
iPatch = find(iMat);
[i, j] = ind2sub(size(img), iPatch);
patches = zeros(prod(patchSize), length(iPatch));
for k = 1:length(iPatch)
	patch = imgExt(i(k):i(k)+a-1, j(k):j(k)+b-1);
	patches(:, k) = patch(:);
end

end % img2patch_index

author: JeffFessler

repo: mirt


text:
function [alpha,beta]=moments23TermRecur(mu)
%%MOMENTS23TERMRECUR Given 2*N+1 moments of a scalar weighting function
%          from 0 to 2*N, the kth moment is
%          mu(k+1)=integral_a^b w(x)*x^k dx where w(x) is a weighting
%          function and a and b need not be finite, obtain the three term
%          recusion for orthogonal polynomials that share the same moments.
%          The three term recursion coefficients cna be used with the
%          orthoPolyZerosFromRecur to obtain quadrature points and weights
%          when using the same weighting function. The three-term recursion
%          has the form:
%          x*p_{k-1}(x)=beta_{k-1}*p_{k-2}(x)+alpha_k*p_{k-1}+beta_kp_k(x)
%          for k=1,...,N starting with p_{-1}(x)=0. This can be useful for
%          obtaining cubature points for integration over continuous
%          probability distributions with easy-to-express noncentral
%          moments (e.g. from the moment generting function of the
%          distribution.
%
%INPUTS: mu A length (2*N+1) array holding the moments from orders 0 to 2*N
%           of the desired weighting function (which could be, but needn't
%           be, a continuous probability distribution).
%
%OUTPUTS: alpha A length N vector of coefficients.
%          beta A length N-1 vector of coefficients.S
%
%The algorithm is taken from Section 4 of [1]. To use these with the
%orthoPolyZerosFromRecur function, one must pass alpha, beta, and mu(1).
%See the example below.
%
%EXAMPLE:
%In this example, we get the values for a three-term recursion of the
%standard normal distribution. We then plus the values into
%orthoPolyZerosFromRecur and get cubature points for integration voer the
%distribution. We then show that these points are identical tot he ones
%obtained from the quadraturePoints1D function (differences are within
%finite precision limits).
% meanVal=0;
% variance=1;
% N=5;
% mu=zeros(2*N+1,1);
% for k=0:2*N
%     mu(k+1)=GaussianD.momentGenFun(meanVal,variance,k);
% end
% [alpha,beta]=moments23TermRecur(mu);
% [xi,w]=orthoPolyZerosFromRecur(N,alpha,beta,[],mu(1));
% [xiGauss,wGauss]=quadraturePoints1D(N,0);
% max(abs((xi-xiGauss)))
% max(abs(w-wGauss))
%
%REFERENCES:
%[1] G. H. Golub and J. H. Welsh, "Calculation of Gauss quadrature rules,"
%    Mathematics of Computation, vol. 23, pp. 221-230, 1969.
%
%January 2021 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

N=(length(mu)-1)/2;

if(fix(N)~=N)
    error('This function requires an odd number of moments.')
end

%The Gram Matrix in Section 4.
M=zeros(N+1,N+1);
for i1=1:(N+1)
    for i2=1:(N+1)
        M(i1,i2)=mu(i1+i2-1);
    end
end

%Using cholSemiDef instead of chol, it is more robust to semidefinite
%matrices.
R=cholSemiDef(M,'upper');

%Equation 4.3 in [1] for both alpha and beta.
alpha=zeros(N,1);
alpha(1)=R(1,2)/R(1,1);
for i=2:N
    alpha(i)=R(i,i+1)/R(i,i)-R(i-1,i)/R(i-1,i-1);
end

beta=zeros(N-1,1);
for i=1:(N-1)
    beta(i)=R(i+1,i+1)/R(i,i);
end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
% Artificial data motivated by mass-preserving imaging
%==============================================================================

checkDataFile
if expfileExists, return; end;

viewPara = {'viewImage','viewImage2D','colormap',gray(256)};
imgPara  = {'imgModel','splineInter','regularizer','moments','theta',1e-2};
traPara  = {'trafo','affine2D'};
disPara  = {'distance','SSD'};
regPara  = {'regularizer', 'mbHyperElastic','alpha',1,...
  'alphaLength',1000,'alphaArea',0,'alphaVolume',100};

expfile = jpgs2data('','Gauss-T.jpg','Gauss-R.jpg', ...
  'omega',[0,1,0,1],'m',[256,256],...
  'viewPara',viewPara,'imgPara',imgPara);

save(expfile,'-append','viewPara','imgPara','traPara','disPara','regPara');
checkDataFile
%==============================================================================