text:
% Using A Binary Sensor Mask Example
%
% This example demonstrates how to use a binary sensor mask for the
% detection of the pressure field generated by an initial pressure
% distribution within a two-dimensional homogeneous propagation medium. It
% builds on the Homogeneous Propagation Medium Example.
%
% author: Bradley Treeby
% date: 29th June 2009
% last update: 22nd September 2012
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox

% This file is part of k-Wave. k-Wave 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.
%
% k-Wave 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 Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.

clear all;

% =========================================================================
% SIMULATION
% =========================================================================

% create the computational grid
Nx = 128;           % number of grid points in the x (row) direction
Ny = 128;           % number of grid points in the y (column) direction
dx = 0.1e-3;        % grid point spacing in the x direction [m]
dy = 0.1e-3;        % grid point spacing in the y direction [m]
kgrid = makeGrid(Nx, dx, Ny, dy);

% define the properties of the propagation medium
medium.sound_speed = 1500; 	% [m/s]
medium.alpha_coeff = 0.75;  % [dB/(MHz^y cm)]
medium.alpha_power = 1.5;

% create initial pressure distribution using makeDisc
disc_magnitude = 5; % [Pa]
disc_x_pos = 50;    % [grid points]
disc_y_pos = 50;    % [grid points]
disc_radius = 8;    % [grid points]
disc_1 = disc_magnitude*makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius);

disc_magnitude = 3; % [Pa]
disc_x_pos = 80;    % [grid points]
disc_y_pos = 60;    % [grid points]
disc_radius = 5;    % [grid points]
disc_2 = disc_magnitude*makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius);

source.p0 = disc_1 + disc_2;

% define a binary sensor mask
sensor_x_pos = Nx/2;        % [grid points]
sensor_y_pos = Ny/2;        % [grid points]
sensor_radius = Nx/2 - 22;  % [grid points]
sensor_arc_angle = 3*pi/2;  % [radians]
sensor.mask = makeCircle(Nx, Ny, sensor_x_pos, sensor_y_pos, sensor_radius, sensor_arc_angle);

% run the simulation
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor);

% reorder the simulation data
sensor_data_reordered = reorderSensorData(kgrid, sensor, sensor_data);

% =========================================================================
% VISUALISATION
% =========================================================================

% plot the initial pressure and sensor distribution
figure;
imagesc(kgrid.y_vec*1e3, kgrid.x_vec*1e3, source.p0 + disc_magnitude*sensor.mask, [-1 1]);
colormap(getColorMap);
ylabel('x-position [mm]');
xlabel('y-position [mm]');
axis image;

% plot the simulated sensor data
figure;
imagesc(sensor_data, [-1, 1]);
colormap(getColorMap);
ylabel('Sensor Position');
xlabel('Time Step');
colorbar;

% plot the re-ordered sensor data
figure;
imagesc(sensor_data_reordered, [-1, 1]);
colormap(getColorMap);
ylabel('Sensor Position');
xlabel('Time Step');
colorbar;


text:
function point = geod2cart(src, curve, normal)
%GEOD2CART Convert geodesic coordinates to cartesian coord.
%
%   PT2 = geod2cart(PT1, CURVE, NORMAL)
%   CURVE and NORMAL are both [N*2] array with the same length, and
%   represent positions of the curve, and normal to each point.
%   PT1 is the point to transform, in geodesic  coordinate (first coord is
%   distance from the curve start, and second coord is distance between
%   point and curve).
%
%   The function return the coordinate of PT1 in the same coordinate system
%   than for the curve.
%
%   TODO : add processing of points not projected on the curve.
%   -> use the closest end
%
%   See also
%   polylines2d, cart2geod, curveLength

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2004-04-08
% Copyright 2004-2022 INRA - TPV URPOI - BIA IMASTE

t = parametrize(curve);
N = size(src, 1);
ind = zeros(N, 1);
for i=1:N
    indices = find(t>=src(i,1));
    ind(i) = indices(1);
end

theta = lineAngle([zeros(N,1) zeros(N,1) normal(ind,:)]);
d = src(:,2);
point = [curve(ind,1)+d.*cos(theta), curve(ind,2)+d.*sin(theta)];

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function result=RS_DF_Test(rdata1,rtrue1,grid1,hstep,el,bootMC)

% Construct PITs
pit = [];
T = size(rdata1,1);

for i = 1:T

    ygrid = grid1;
    rdata = rdata1(i,1:length(ygrid))';

    fitdens1 = fitdist(ygrid, 'normal', 'frequency',round(rdata));

    z1 = normcdf(rtrue1(i,:),fitdens1.mu, fitdens1.sigma);

    pit = [pit; z1];
end

% Construct Histogram of PITS
bin  = 10;
m    = length(pit);
rvec = (0:0.001:1);
hs   = histc(pit,0:1/bin:1);

result.rvec      = rvec;
result.m         = m;
result.bin       = bin;
result.histogram = hs;

% Test Statistic
for r = 1:size(rvec,2)
    ecdf(:,r) = mean(pit < rvec(:,r));
end
result.ecdf=ecdf;

% Critical Values
if hstep == 1
    results1CS = bear.RS_DF_rstest(pit);
    table1 = results1CS(:,1:2);
elseif hstep >= 1
%     table1 = RS_DF_CVfinalbootstrap(el,bootMC,pit,rvec); % 20.9.2018
    table1 = bear.RS_DF_CVfinalbootstrapInoue(el,bootMC,pit,rvec);
end

result.output=round(table1.*1000)./1000;

author: european-central-bank

repo: BEAR-toolbox


text:
function out = OrthNorm(n)
% =======================================================================
% Compute a [n x n] random orthonormal matrix
% =======================================================================
% out = OrthNorm(n)
% -----------------------------------------------------------------------
% INPUT
%   - n: size of the matrix
% -----------------------------------------------------------------------
% OUTPUT
%   - out: [n x n] random orthonormal matrix
% =======================================================================
% Ambrogio Cesa Bianchi, March 2015
% ambrogio.cesabianchi@gmail.com

% [n x n]  matrix of N(0,1) random variables
X = randn(n,n);

% QR decomposition of X
[Q, ~] = qr(X);

% Check precision
if sum(sum(Q*Q'))> n + 1.0e-5
    error('Q*transpose(Q) is not equal to identity')
end

% Random orthonormal matrix such that Q*Q'=I
out = Q;

author: ambropo

repo: VAR-Toolbox

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


text:
function [c,ceq] = nl_constraints(t, x, u, mpcModel, varargin) %add battery SOC constraints

    c(1) = x(2)-mpcModel.battery.range(2);
    c(2) = mpcModel.battery.range(1)-x(2);

    ceq =  [];

end

author: juchengquan


text:
function [identified]=irf_disp_IV(n,endo,IRFperiods,IRFt,irf_estimates,D_estimates,gamma_estimates,pref,signreslabels,strctident)

%determine number of identified shocks
if IRFt==4||IRFt==6
Index=find(~contains(signreslabels,'shock'));
identified=min(Index)-1;
identified=n-numel(Index);
elseif IRFt==5
identified=1;
signreslabels{1,1} = strcat('Shock identified by IV (',strctident.Instrument,')');
else
identified = n;
end

count=1;
namecount=1;
irf_estimates = irf_estimates';
hf = figure;
   if IRFt==1
   irfname='impulse response functions (no structural identifcation)';
   elseif IRFt==2
   irfname='impulse response functions (structural identification by Cholesky ordering)';
   elseif IRFt==3
   irfname='impulse response functions (structural identification by triangular factorisation)';
   elseif IRFt==4
   irfname=['impulse response functions (structural identification by ',strctident.hbartext_signres,strctident.hbartext_zerores,strctident.hbartext_magnres,strctident.hbartext_relmagnres,strctident.hbartext_FEVDres,strctident.hbartext_CorrelInstrumentShock,':::',' restrictions)'];
   irfname=erase(irfname,', :::'); % delete the last ,
   elseif IRFt==5
   irfname=['impulse response functions (structural identification by IV ',strctident.Instrument,')'];
   elseif IRFt==6
   irfname_temp=[strctident.hbartext_signres,strctident.hbartext_zerores,strctident.hbartext_magnres,strctident.hbartext_relmagnres,strctident.hbartext_FEVDres,strctident.hbartext_CorrelInstrumentShock,':::',' restrictions)'];
   irfname_temp=erase(irfname_temp,', :::'); % delete the last ,
   irfname=['impulse response functions (structural identification by IV ',strctident.Instrument,' & ',irfname_temp];
   end
   set(hf,'name',irfname);
set(hf,'name',irfname);
set(hf,'Color',[0.9 0.9 0.9]);
%get(0,'ScreenSize');
set(hf,'Units','normalized');
% unitsperplot=0.2*identified*n/6;
% if unitsperplot>1
% set(hf,'position',[0,0,1,1])
% else
% set(hf,'position',[0,0,0.2*identified,0.2*identified])
% end

%%%%% set(hf,'position',[250,0,1800,1200])
%set(hf, 'Color', [1,1,1]);
%fontsize=ceil(18*(6/n)); %if n=6 fontsize=18;
fontsize=11;

for jj=1:identified
    for ii=1:n
        count;
        subplot(n,identified,jj+identified*(ii-1));
        hold on
Xpatch=[(1:IRFperiods) (IRFperiods:-1:1)];
Ypatch=[irf_estimates{jj,ii}(1,:) fliplr(irf_estimates{jj,ii}(3,:))];
IRFpatch=patch(Xpatch,Ypatch,[1 0.5 0]);
set(IRFpatch,'facealpha',0.9);
set(IRFpatch,'edgecolor','none');
plot(irf_estimates{jj,ii}(2,:),'Color',[1 0 0],'LineWidth',3);
plot([1,IRFperiods],[0 0],'k--');
hold off
minband=min(irf_estimates{jj,ii}(1,:));
maxband=max(irf_estimates{jj,ii}(3,:));
space=maxband-minband;
Ymin=minband-0.2*space;
Ymax=maxband+0.2*space;
set(gca,'XLim',[1 IRFperiods],'YLim',[Ymin Ymax],'FontName','Times New Roman','Fontsize',fontsize);
% top labels
subplotcolumn=jj+identified*(ii-1);
if subplotcolumn <= identified
title(signreslabels{subplotcolumn,1},'FontWeight','normal','Interpreter','none','Fontsize',fontsize);
end
if count<=n
   ylabel(endo{namecount,1},'FontWeight','normal','Interpreter','none','Fontsize',fontsize);
   namecount=namecount+1;
end
count=count+1;
end
end
% top supertitle
ax=axes('Units','Normal','Position',[.11 .075 .85 .88],'Visible','off');
set(get(ax,'Title'),'Visible','on')
%%%%% title('Shock:','FontSize',11,'FontName','Times New Roman','FontWeight','normal');
% side supertitle
%%%%% ylabel('Response of:','FontSize',16,'FontName','Times New Roman','FontWeight','normal');
set(get(ax,'Ylabel'),'Visible','on')
    set(gcf,'PaperPositionMode','Auto')
    %fname=strcat(pref.datapath, '\results\');
    %saveas(gcf,[fname,PrintName],'epsc')
    %saveas(gcf,'IRFs','png');
% then display the results for D and gamma, if a structural decomposition was selected

if IRFt==2 || IRFt==3 || IRFt==4 || IRFt==5 || IRFt==6

filelocation=fullfile(pref.results_path, [pref.results_sub '.txt']);
fid=fopen(filelocation,'at');

fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');

svarinfo1=['D (structural decomposition matrix): posterior estimates'];
fprintf('%s\n',svarinfo1);
fprintf(fid,'%s\n',svarinfo1);

% recover D
D=reshape(D_estimates,n,n);
% calculate the (integer) length of the largest number in D, for formatting purpose
width=length(sprintf('%d',floor(max(abs(bear.vec(D))))));
% add a separator, a potential minus sign and three digits (total=5) to obtain the total space for each entry in the matrix
width=width+5;
for ii=1:n
temp=[];
   for jj=1:n
   % convert matrix entry into string
   number=num2str(D(ii,jj),'% .3f');
      % pad potential missing blanks
      while numel(number)<width
      number=[' ' number];
      end
   number=[number '  '];
   temp=[temp number];
   end
fprintf('%s\n',temp);
fprintf(fid,'%s\n',temp);
end

fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');

svarinfo2=['gamma (structural disturbances covariance matrix): posterior estimates'];
fprintf('%s\n',svarinfo2);
fprintf(fid,'%s\n',svarinfo2);

% recover gamma
gamma=reshape(gamma_estimates,n,n);
% calculate the (integer) length of the largest number in D, for formatting purpose
width=length(sprintf('%d',floor(max(abs(bear.vec(gamma))))));
% add a separator, a potential minus sign and three digits (total=5) to obtain the total space for each entry in the matrix
width=width+5;
for ii=1:n
temp=[];
   for jj=1:n
   % convert matrix entry into string
   number=num2str(gamma(ii,jj),'% .3f');
      % pad potential missing blanks
      while numel(number)<width
      number=[' ' number];
      end
   number=[number '  '];
   temp=[temp number];
   end
fprintf('%s\n',temp);
fprintf(fid,'%s\n',temp);
end

fclose(fid);
end

% finally, record results in excel
% retranspose the cell of records
irf_estimates=irf_estimates';
bear.data.excelrecord4fcn(identified, IRFperiods, IRFt, endo, strctident, irf_estimates, n, pref)


text:
function test_eigen_test01 ( )

%*****************************************************************************80
%
%% TEST_EIGEN_TEST01 tests R8SYMM_TEST.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 February 2012
%
%  Author:
%
%    John Burkardt
%
  n = 100;
  bin_num = 10;
  lambda_dev = 1.0;
  lambda_mean = 1.0;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST_EIGEN_TEST01\n' );
  fprintf ( 1, '  R8SYMM_TEST makes an arbitrary size symmetric matrix\n' );
  fprintf ( 1, '  with known eigenvalues and eigenvectors.\n' );

  [ a, q, lambda, seed ] = r8symm_test ( n, lambda_mean, lambda_dev, seed );
%
%  Get the eigenvalue range.
%
  lambda_min = min ( lambda(1:n) );
  lambda_max = max ( lambda(1:n) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  LAMBDA_MIN = %g\n', lambda_min );
  fprintf ( 1, '  LAMBDA_MAX = %g\n', lambda_max );
%
%  Bin the eigenvalues.
%
  [ bin, bin_limit ] = r8vec_bin ( n, lambda, bin_num, lambda_min, lambda_max );

  r8bin_print ( bin_num, bin, bin_limit, '  Lambda bins:' );

  if ( 0 )
    r8mat_print ( n, n, a, '  The matrix A:' );
  end

  if ( 0 )
    r8mat_print ( n, n, q, '  The eigenvector matrix Q:' );
  end

  aq = a * q;

  for j = 1 : n
    lambda2(j) = norm ( aq(1:n,j) );
  end

  r8vec2_print ( n, lambda, lambda2, '  LAMBDA versus column norms of A*Q:' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function DEM_coupled_oscillators
% Dual estimation of the Lorenz system: Cross-validation of Laplace schemes
%__________________________________________________________________________
% This routine illustrates the inversion of a loosely coupled oscillator
% model using generalised filtering. In this example, three regions are
% coupled in terms of their amplitude and phase in a hierarchical fashion.
% Data are generated under a particular set of parameters. The timeseries
% are then transformed using a Hilbert transform into the corresponding
% analytic signal. This then constitutes the data feature for subsequent
% inversion using generalised filtering; here, in four generalised
% coordinates of motion. By assuming fairly precise priors on the amplitude
% of random fluctuations one can recover the parameters and use the
% posterior density for subsequent Bayesian model comparison. In this
% example, we used Bayesian model reduction to assess the evidence for
% models with and without amplitude or phase coupling.
%
% The parameters and orders of this example have been optimised to provide
% proof of principle this sort of  model can be inverted using generalised
% filtering.  The sensitivity to these parameters and orders can be
% assessed numerically by editing the code.
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: DEM_coupled_oscillators.m 7679 2019-10-24 15:54:07Z spm $

% specify states and parameters
%==========================================================================
N     = 128;                             % number of time points
n     = 3;                               % number of sources (oscillators)
Hz    = 8;                               % characteristic frequency (Hz)
dt    = 1/64;                            % sampling interval (sec)

% model states (where hidden states comprise phase differences)
%--------------------------------------------------------------------------
x.r   = zeros(n,1);                      % amplitude
x.p   = zeros(n,1);                      % phase  differences
x.w   = zeros(1,1);                      % phase (common to sources)

% model parameters
%--------------------------------------------------------------------------
A     = spm_speye(n,n,-1);               % form of adjacency matrix
P.L   = diag(ones(n,1));                 % lead field (measurement mapping)
P.Ap  = A/16 - A'/32 - eye(n,n)/16;      % amplitude coupling
P.Ar  = A/16 + A'/4;                     % phase coupling
P.C   = sparse(1,1,1,n,1)/8;             % exogenous input to first source
P.r   = 1/8;                             % weak amplitude
P.w   = 2*pi*Hz*dt;                      % intrinsic frequency

% observation function (to generate timeseries)
%--------------------------------------------------------------------------
g = @(x,v,P) P.L*((x.r).*cos(x.p + x.w));

% equations of motion (simplified coupled oscillator model)
%--------------------------------------------------------------------------
f = @(x,v,P) [P.Ap*(x.r - P.r) + P.C*v;
              sum(P.Ar.*sin(bsxfun(@minus,x.p,x.p')),2) - P.C*v; ...
              P.w];

% causes or exogenous input (a Gaussian function of peristimulus time)
%--------------------------------------------------------------------------
U = exp(-((1:N) - N/2).^2/((N/8)^2));    % exogenous input
T = (1:N)*dt;                            % sample times (seconds)

% parameters for generalised filtering (see spm_LAP)
%--------------------------------------------------------------------------
E.n     = 4;                             % embedding dimension
E.d     = 1;                             % data embedding
E.nN    = 8;                             % number of iterations
E.s     = 1/2;                           % smoothness of fluctuations

% first level state space model
%--------------------------------------------------------------------------
M(1).E  = E;                             % filtering parameters
M(1).x  = x;                             % initial states
M(1).f  = f;                             % equations of motion
M(1).g  = g;                             % observation mapping
M(1).pE = P;                             % model parameters
M(1).V  = exp(12);                       % precision of observation noise
M(1).W  = exp(12);                       % precision of state noise

% second level - causes or exogenous forcing term
%--------------------------------------------------------------------------
M(2).v  = 0;                             % initial causes
M(2).V  = exp(16);                       % precision of exogenous causes

% create data with known parameters (P)
%==========================================================================
DEM = spm_DEM_generate(M,U,P);

% transform analytic signal to create a new data feature
%==========================================================================

% analytic signal (via Hilbert transform)
%--------------------------------------------------------------------------
Y  = spm_hilbert(full(DEM.Y)');          % analytic signal
Yr = abs(Y)';                            % amplitude
Yp = unwrap(angle(Y));                   % phase
Yp = Yp' - ones(n,1)*(1:N)*P.w;          % phase difference
Y  = [Yr; Yp];                           % analytic data feature

% show synthetic data,latent states and exogenous input
%--------------------------------------------------------------------------
spm_figure('GetWin','synthetic data');
spm_DEM_qU(DEM.pU);

subplot(4,2,2), plot(T,DEM.pU.x{1}(1:n,:)')
title('hidden amplitude','FontSize',16)
xlabel('time (seconds)'), spm_axis tight, box off

subplot(4,2,4), plot(T,DEM.pU.x{1}((1:n) + n,:)')
title('phase difference','FontSize',16)
xlabel('time (seconds)'), spm_axis tight, box off

subplot(4,2,6), plot(T,Yr)
title('response amplitude','FontSize',16)
xlabel('time (seconds)'), spm_axis tight, box off

subplot(4,2,8), plot(T,Yp)
title('unwrapped phase','FontSize',16)
xlabel('time (seconds)'), spm_axis tight, box off
drawnow

% Now try to recover model parameters from data features
%==========================================================================

% change observation function (g) to generate analytic signal
%--------------------------------------------------------------------------
g = @(x,v,P) [P.L*x.r; x.p];

% initialization of priors over parameters
%--------------------------------------------------------------------------
pE       = P;                            % prior parameters
pC       = spm_zeros(P);                 % prior variance

pE.Ar    = zeros(n,n);                   % set prior phase and amplitude
pE.Ap    = -speye(n,n)/16;               % coupling parameters to 0
pC.Ar    = (P.Ar ~= 0);                  % and set the prior variance to 1
pC.Ap    = (P.Ap ~= 0);

Vr       = ones(1,n)*8;                  % log precision of sampling noise
Vp       = ones(1,n)*8;                  % and state noise

% place new observation function and priors in generative model
%--------------------------------------------------------------------------
DEM.M(1).g  = g;
DEM.M(1).pE = pE;
DEM.M(1).pC = pC;
DEM.M(1).V  = exp([Vr Vp]);
DEM.M(1).W  = exp([Vr Vp 32]);           % use precise beliefs about time

% data and known input; removing initial time points to suppress artefacts
%--------------------------------------------------------------------------
DEM.Y = Y(:,8:end);
DEM.U = U(:,8:end);

% Inversion using generalised filtering
%==========================================================================
LAP   = spm_DEM(DEM);

% Show parameters
%--------------------------------------------------------------------------
spm_figure('GetWin','Parameters'); clf; spm_DEM_qP(LAP.qP,LAP.pP)
subplot(2,1,1),legend('mean','90% CI','Location','North'), legend(gca,'boxoff')
title('Estimated and true (black) parameters','FontSize',16)

% use Bayesian model reduction to test different hypotheses
%==========================================================================
model{1} = 'no coupling';
model{2} = 'no amplitude coupling';
model{3} = 'no phase coupling';
model{4} = 'Full model';

% apply precise shrinkage priors to off-diagonal coupling elements
%--------------------------------------------------------------------------
PC{1} = pC; PC{1}.Ar = diag(diag(pC.Ar)); PC{1}.Ap = diag(diag(pC.Ap));
PC{2} = pC; PC{2}.Ap = diag(diag(pC.Ap));
PC{3} = pC; PC{3}.Ar = diag(diag(pC.Ar));
PC{4} = pC;

%  evaluate the evidence for these new models or prior constraints
%--------------------------------------------------------------------------
qE    = LAP.qP.P{1};
qC    = LAP.qP.C;
pE    = LAP.M(1).pE;
pC    = LAP.M(1).pC;
for m = 1:numel(PC)
    rC     = diag(spm_vec(PC{m}));
    F(m,1) = spm_log_evidence(qE,qC,pE,pC,pE,rC);
end

% report marginal log likelihood or evidence
%--------------------------------------------------------------------------
F = F - min(F);

spm_figure('GetWin','Model Comparison');clf;
subplot(2,2,1), bar(F,'c')
title('Log evidence','FontSize',16)
xlabel(model), axis square, box off

subplot(2,2,2), bar(spm_softmax(F(:)),'c')
title('Probability','FontSize',16)
xlabel(model), axis square, box off

author: spm


text:
function [Neurons] = neural_gas(D,n,epochs,alpha0,lambda0)

%NEURAL_GAS Quantizes the data space using the neural gas algorithm.
%
% Neurons = neural_gas(D, n, epochs, [alpha0], [lambda0])
%
%   C = neural_gas(D,50,10);
%   sM = som_map_struct(sD);
%   sM.codebook = neural_gas(sD,size(sM.codebook,1),10);
%
%  Input and output arguments ([]'s are optional):
%   D          (matrix) the data matrix, size dlen x dim
%              (struct) a data struct
%   n          (scalar) the number of neurons
%   epochs     (scalar) the number of training epochs (the number of
%                       training steps is dlen*epochs)
%   [alpha0]   (scalar) initial step size, 0.5 by default
%   [lambda0]  (scalar) initial decay constant, n/2 by default
%
%   Neurons    (matrix) the neuron matrix, size n x dim
%
% See also SOM_MAKE, KMEANS.

% References:
%  T.M.Martinetz, S.G.Berkovich, and K.J.Schulten. "Neural-gas" network
%  for vector quantization and its application to time-series prediction.
%  IEEE Transactions on Neural Networks, 4(4):558-569, 1993.

% Contributed to SOM Toolbox vs2, February 2nd, 2000 by Juha Vesanto
% Copyright (c) by Juha Vesanto
% http://www.cis.hut.fi/projects/somtoolbox/

% juuso 101297 020200

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Check arguments and initialize

error(nargchk(3, 5, nargin));  % check the number of input arguments

if isstruct(D), D = D.data; end
[dlen,dim] = size(D);
Neurons = (rand(n,dim)-0.5)*10e-5; % small initial values
train_len = epochs*dlen;

if nargin<4 || isempty(alpha0) || isnan(alpha0), alpha0 = 0.5; end
if nargin<5 || isempty(lambda0) || isnan(lambda0), lambda0 = n/2; end

% random sample order
rand('state',sum(100*clock));
sample_inds = ceil(dlen*rand(train_len,1));

% lambda
lambda = lambda0 * (0.01/lambda0).^([0:(train_len-1)]/train_len);

% alpha
alpha = alpha0 * (0.005/alpha0).^([0:(train_len-1)]/train_len);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Action

for i=1:train_len,

  % sample vector
  x = D(sample_inds(i),:); % sample vector
  known = ~isnan(x);       % its known components
  X = x(ones(n,1),known);  % we'll need this

  % neighborhood ranking
  Dx = Neurons(:,known) - X;  % difference between vector and all map units
  [~, inds] = sort((Dx.^2)*known'); % 1-BMU, 2-BMU, etc.
  ranking(inds) = [0:(n-1)];
  h = exp(-ranking/lambda(i));
  H = h(ones(length(known),1),:)';

  % update
  Neurons = Neurons + alpha(i)*H.*(x(ones(n,1),known) - Neurons(:,known));

  % track
  fprintf(1,'%d / %d \r',i,train_len);
  if 0 && mod(i,50) == 0,
    hold off, plot3(D(:,1),D(:,2),D(:,3),'bo')
    hold on, plot3(Neurons(:,1),Neurons(:,2),Neurons(:,3),'r+')
    drawnow
  end
end

fprintf(1,'\n');

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

author: ilarinieminen

repo: SOM-Toolbox


text:
% removERP.m removes the averaged ERP from every single epoch of its corresponding epoched EEG dataset
%
% Input:
%
% EEG       - epoched EEG dataset
% ERP       - averaged ERP (ERPset)
%
% Output:
%
% EEG       - epoched EEG dataset
%
% IMPORTANT:
% - epoched EEG dataset must have assigned bins (e.g. processed by ERPLAB).
% - epoched EEG dataset and averaged ERP (ERPset) must have the same amount of channels, samples and sample rate.
% - the time-locked event must be at the same sample location in both epoched EEG dataset and averaged ERP (ERPset).
%
%
%
% *** This function is part of ERPLAB Toolbox ***
% Author: Javier Lopez-Calderon & Johanna Kreither
% Center for Mind and Brain
% University of California, Davis,
% Davis, CA
% July-August 2014

function EEG = removERP(EEG, ERP)

nepoch   = EEG.trials;
nchaneeg = EEG.nbchan;
npntseeg = EEG.pnts;
srateeeg = EEG.srate;
eegzerolatpos = find(EEG.times == 0,1,'first'); % catch zero-time locked code position,

nchanerp = ERP.nchan;
npntserp = ERP.pnts;
srateerp = ERP.srate;
erpzerolatpos = find(ERP.times == 0,1,'first'); % catch zero-time locked code position,

if nchaneeg~=nchanerp
      error('ERPLAB:error', 'EEG has %g channels and ERP has %g canales!', nchaneeg, nchanerp)
end
if npntseeg~=npntserp
      error('ERPLAB:error', 'EEG has %g points and ERP has %g points!', npntseeg, npntserp)
end
if srateeeg~=srateerp
      error('ERPLAB:error', 'EEG has a srate of %g sps and ERP has a srate of %g sps!', srateeeg, srateerp)
end
if eegzerolatpos~=erpzerolatpos
      error('ERPLAB:error', 'The time-locked event is not located at the same sample for EEG and ERP.')
end
for k=1:nepoch
      latenarray = EEG.epoch(k).eventlatency;
      biniarray  = EEG.epoch(k).eventbini;

      if iscell(latenarray)
            latenarray = cell2mat(latenarray);
      end

      indxtimelock = find(latenarray == 0,1,'first'); % catch zero-time locked code position,
      bini  = biniarray(indxtimelock);

      if isempty(bini)
            error('ERPLAB:error', 'Epoch %g does not habe a bin assigned (empty).', k)
      end
      if iscell(bini)
            bini = cell2mat(bini);
      end
      if bini<1
            error('ERPLAB:error', 'Epoch %g does not habe a bin assigned (-1).', k)
      end
      if length(bini)>1
            error('ERPLAB:error', 'Epoch %g has more than one bin assigned.', k)
      end

      EEG.data(:,:,k) = EEG.data(:,:,k) - ERP.bindata(:,:, bini);
end

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
function Mh = oprBuilderHmx(opr,k,Xunk,Mx,X,Xnrm,Yunk,My,Y,Ynrm,tol)
%+========================================================================+
%|                                                                        |
%|            OPENOPR - LIBRARY FOR SPECIFIC OPERATORS IN BEM             |
%|           openOpr is part of the GYPSILAB toolbox for Matlab           |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal (c) 2017-2019.                             |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : oprBuilderHmx.m                               |
%|    #    |   VERSION    : 0.61                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 05.09.2018                                    |
%| ( === ) |   SYNOPSIS   : Finite element builder with low-rank          |
%|  `---'  |                approximation for handle function             |
%+========================================================================+

% Initialize H-Matrix
Mh = hmx(Xunk,Yunk,tol);

% Admissibility
[isfar,Xdim,Ydim] = hmxFar(Mh);

% Compression for far distances
if isfar
    % ACA compression for quadrature matrix
    if strcmp(opr,'nxK')
        [A,B,flag] = oprACAv(opr,k,X,Xnrm,Y,Ynrm,tol);
    else
        [A,B,flag] = oprACA(opr,k,X,Xnrm,Y,Ynrm,tol);
    end

    % Finite element integration
    if flag
        % Integration
        if strcmp(opr,'H') || strcmp(opr,'T')
            for n = 1:length(Mx)
                Mx{n} = Mx{n} * A;
                My{n} = B * My{n};
            end
            A = cell2mat(Mx);
            B = cell2mat(My');
        elseif strcmp(opr,'nxK')
            A = [ Mx{1}*A(:,:,2) , - Mx{1}*A(:,:,3) , ...
                Mx{2}*A(:,:,3) , - Mx{2}*A(:,:,1) , ...
                Mx{3}*A(:,:,1) , - Mx{3}*A(:,:,2) ] ;
            B = [ B(:,:,2)*My{3} ; B(:,:,3)*My{2} ; ...
                B(:,:,3)*My{1} ; B(:,:,1)*My{3} ; ...
                B(:,:,1)*My{2} ; B(:,:,2)*My{1} ] ;
        else
            A = Mx * A;
            B = B * My;
        end

        % Recompression
        [A,B] = hmxQRSVD(A,B,tol);
    end

else
    flag = 0;
end

%%% Compression
if flag
    % Type
    Mh.typ = 1;

    % Low-rank
    Mh.dat = {A,B};

%%%% Full or sparse for smallest box (stopping criterion)
elseif (min(size(Mh)) < 100)
    % Type
    Mh.typ = 2;

    % Quadrature matrix
    if strcmp(opr,'nxK')
        Gxy = zeros(size(X,1),size(Y,1),3);
        for j = 1:size(Gxy,2)
            Gxy(:,j,:) = oprGreenKernel(opr,k,X,Xnrm,Y(j,:),Ynrm(j,:));
        end
    else
        Gxy = zeros(size(X,1),size(Y,1));
        for j = 1:size(Gxy,2)
            Gxy(:,j) = oprGreenKernel(opr,k,X,Xnrm,Y(j,:),Ynrm(j,:));
        end
    end

    % Finite element integration
    if strcmp(opr,'H') || strcmp(opr,'T')
        Mh.dat = zeros(size(Mh));
        for n = 1:length(Mx)
            Mh.dat = Mh.dat + Mx{n} * Gxy * My{n};
        end
    elseif strcmp(opr,'nxK')
        Mh.dat = Mx{1} * Gxy(:,:,2) * My{3} - Mx{1} * Gxy(:,:,3) * My{2} + ...
            Mx{2} * Gxy(:,:,3) * My{1} - Mx{2} * Gxy(:,:,1) * My{3} + ...
            Mx{3} * Gxy(:,:,1) * My{2} - Mx{3} * Gxy(:,:,2) * My{1} ;
    else
        Mh.dat = Mx * Gxy * My;
    end

%%% H-Matrix (recursion)
else
    % Type
    Mh.typ = 0;

    % Subdivision for X
    [I1,I2] = hmxSubdivide(Xunk,Xdim);
    Mh.row  = {I1 , I1 , I2 , I2};

    % Subdivision for Y
    [I1,I2] = hmxSubdivide(Yunk,Ydim);
    Mh.col  = {I1 , I2 , I1 , I2};

    % H-Matrix (recursion)
    for i = 1:4
        % Dof indices
        Ir = Mh.row{i};
        Ic = Mh.col{i};

        % Fem matrix subdivision and quadratures points indices
        [Mxchd,Ix] = oprSubdivideCell(Mx,Ir,'left');
        [Mychd,Iy] = oprSubdivideCell(My,Ic,'right');

        % Recursion
        Mh.chd{i} = oprBuilderHmx(opr,k,Xunk(Ir,:),Mxchd,X(Ix,:),Xnrm(Ix,:),...
            Yunk(Ic,:),Mychd,Y(Iy,:),Ynrm(Iy,:),tol);
    end

    % Fusion
    Mh = hmxFusion(Mh);
end
end

%         % Add H-Vector
%         Mv(Ix) = Mv(Ix) + tmp
%
%
% % Corrective Matrix-vector product for Stokes Stresslet
% if strcmp(opr,('Tc'))
%     % H-Matrix (recursion)
%     if (Mh.typ == 0)
%
%     % Compressed leaf
%     elseif (Mh.typ == 1)
%          TODO
% %         Mv = Mx\(A*sum(B,2));
%
%     % Full leaf
%     elseif (Mh.typ == 2)
%         Mv = sum(Gxy*My,2);
%     end
% end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [xcon,ycon] = consolidator13(x,y,aggregation_mode,tol,tolerance_style)
% consolidator: consolidate "replicates" in x, also aggregate corresponding y
% usage: [xcon,ycon] = consolidator13(x,y,aggregation_mode,tol,tolerance_style)
%
% arguments: (input)
%  x - rectangular array of data to be consolidated. If multiple (p)
%      columns, then each row of x is interpreted as a single point in a
%      p-dimensional space. (x may be character, in which case xcon will
%      be returned as a character array.)
%
%      x CANNOT be complex. If you do have complex data, split it into
%      real and imaginary components as columns of an array.
%
%      If x and y are both ROW vctors, they will be transposed and
%      treated as single columns.
%
%  y - outputs to be aggregated. If y is not supplied (is left empty)
%      then consolidator is similar to unique(x,'rows'), but with a
%      tolerance on which points are distinct. (y may be complex)
%
%      y MUST have the same number of rows as x.
%
%  aggregation_mode - (OPTIONAL) character flag, denotes the method
%      to be used in the aggregation of elements of y. If y was
%      left empty, then aggregation_mode is ignored. Contractions
%      are allowed and capitalization is ignored.
%
%      DEFAULT: 'mean'
%
%      'mean'  --> compute the average y for each consolidated x
%      'sum'   --> compute the sum of y for each consolidated x
%      'min'   --> compute the minimum y for each consolidated x
%      'max'   --> compute the maximum y for each consolidated x
%      'std'   --> compute the standard deviation of y for each x
%      'var'   --> compute the variance of y for each x
%      'prod'  --> compute the product of y for each consolidated x
%      'geom'  --> compute the geometric mean at each consolidated x
%      'harm'  --> compute the harmonic mean at each consolidated x
%      'median' --> compute the median at each consolidated x
%      'count' --> compute the # of consolidated elements 7at each x
%
%      Note: 'count' only returns 1 column, regardless
%      of the size of y. All other operations operate on
%      each column of y independently.
%
%  tol - (OPTIONAL) tolerance to identify replicate elements of x. If
%      x has multiple columns, then the same tolerance is applied to
%      all columns of x.
%
%      DEFAULT (if tolerance_style is 'absolute'): 0
%      DEFAULT (if tolerance_style is 'relative'): 1.e-12
%
%      Note: Tolerances are a complicated thing, especially in
%      higher dimensions and when trying to do it effficiently for
%      high numbers of data points. A complete interpoint distance
%      matrix would be the proper way to do t, but for many thousands
%      of points it would use too much memory. One also has the problem
%      of transitivity. If A is within tol of B, and B is within tol
%      of C, then is A near enough to C? If not, then how is the
%      consolidation to be done? So please view these tolerances
%      with some tolerance of your own.
%
%      Note: relative tolerances near 1 are not meaningful.
%
% tolerance_style - (OPTIONAL) character flag, denotes whether
%      the tolerance criterion is absolute over all columns or
%      a relative tolerance within a single column of x.
%
%      Note: DO NOT use a relative tolerance of 0, nor greater
%      than 1. Even a relative tolerance near 1 is not meaningful.
%
%      DEFAULT: 'absolute'
%
%      'relative' --> Tol is a relative tolerance between the
%         minimum and maximum of each column of x.
%
%      'absolute' --> look for differences of size tol in any
%         column of x.
%
%
% arguments: (output)
%  xcon - consolidated x. Replicates wthin the tolerance are removed.
%      if no y was specified, then consolidation is still done on x.
%
%  ycon - aggregated value as specified by the aggregation_mode.
%
%
% Example usages:
%
% Group means:
%  x = round(rand(1000,1)*5);
%  y = x+randn(size(x));
%  [xg,yg] = consolidator(x,y,'mean');
%  [xg,yg]
%  ans =
%         0    0.1668
%    1.0000    0.9678
%    2.0000    2.0829
%    3.0000    2.9688
%    4.0000    4.0491
%    5.0000    4.8852
%
% Group counts on x
%  x = round(randn(100000,1));
%  [xg,c] = consolidator(x,[],'count');
%  [xg,c]
%  ans =
%         -4          26
%         -3         633
%         -2        5926
%         -1       24391
%          0       38306
%          1       24156
%          2        5982
%          3         559
%          4          21
% Unique(x,'rows'), but with a tolerance
%  x = rand(100,2);
%  xc = consolidator(x,[],[],.05);
%  size(xc)
%  ans =
%      62     2
%
% See also: unique
%
% Author: John D'Errico
% e-mail address: woodchips@rochester.rr.com
% Release: 3
% Release date: 5/2/06

% is it a character array?
if ischar(x)
  charflag = 1;
  x=double(x);
else
  charflag = 0;
end

% check for/supply defaults
if (nargin<5) | isempty(tolerance_style)
  tolerance_style = 'absolute';
else
  tolerance_style = lower(tolerance_style);
end

if (nargin<4) | isempty(tol)
  switch tolerance_style
    case 'absolute'
      tol = 0;
    case 'relative'
      tol = 1.e-12;
    otherwise
      error 'tolerance_style must be either ''absolute'' or ''relative'''
  end
end
if (tol<0) & strcmp(tolerance_style,'absolute')
  error 'Absolute tolerance must be non-negative.'
elseif (tol<=0) & strcmp(tolerance_style,'relative')
  error 'Absolute tolerance must be non-negative.'
end
tol = tol*(1+10*eps);

if (nargin<3) | isempty(aggregation_mode)
  aggregation_mode = 'mean';
else
  aggregation_mode = lower(aggregation_mode);
end
valid_agmode = {'mean' 'sum' 'min' 'max' 'std' ...
  'var' 'prod' 'geom' 'harm' 'count' 'median'};
k = strmatch(aggregation_mode,valid_agmode);
if isempty(k)
  error 'Invalid aggregation_mode'
else
  aggregation_mode = valid_agmode{k};
end

% was y supplied, or empty?
[n,p] = size(x);
if (nargin<2) | isempty(y)
  y = zeros(n,0);
end
% check for mismatch between x and y
[junk,q] = size(y);
if n~=junk
  error 'y must have the same number of rows as x.'
end

% are both x and y row vectors?
if (n == 1)
  x=x';
  n = length(x);
  p = 1;

  if ~isempty(y)
    y=y';
  else
    y=zeros(n,0);
  end

  q = size(y,2);
end

if isempty(y)
  aggregation_mode = 'count';
end

% consolidate elements of x.
% first shift, scale, and then round up.
switch tolerance_style
  case 'absolute'
    if tol>0
      xhat = x - repmat(min(x,[],1),n,1)+tol*eps;
      xhat = ceil(xhat/tol);
    else
      xhat = x;
    end
  case 'relative'
    xhat = x - repmat(min(x,[],1),n,1)+tol*eps;
    xhat = ceil(xhat*diag(1 ./ (tol*max(xhat,[],1))));
end
[xhat,tags] = sortrows(xhat);
x=x(tags,:);
y=y(tags,:);

% count the replicates
iu = [1;any(diff(xhat),2)];
eb = cumsum(iu);
% count is the vector of counts for the consolidated
% x values
% Consolidator code was:  count=accumarray(eb,1).';
count = full(sparse(1,eb,1,1,eb(end)));

% ec is the expanded counts, i.e., counts for the
% unconsolidated x
ec = count(eb);

% special case for aggregation_mode of 'count',
% but we still need to aggregate (using the mean) on x
if strcmp(aggregation_mode,'count')
  ycon = count.';
  q = 0;
else
  ycon = zeros(length(count),q);
end

% loop over the different replicate counts, aggregate x and y
ucount = unique(count);
xcon = repmat(NaN,[length(count),p]);
for k=ucount
  if k==1
    xcon(count==1,:) = x(ec==1,:);
  else
    v=permute(x(ec==k,:),[3 2 1]);
    v=reshape(v,[p,k,prod(size(v))/p/k]);
    v=permute(v,[2 1 3]);
    xcon(count==k,:)=reshape(mean(v),p,prod(size(v))/p/size(v,1)).';
  end

  if q>0
    % aggregate y as specified
    if k==1
      switch aggregation_mode
        case {'std' 'var'}
          ycon(count==1,:) = 0;
        otherwise
          ycon(count==1,:) = y(ec==1,:);
      end
    else
      v=permute(y(ec==k,:),[3 2 1]);
      v=reshape(v,[q,k,prod(size(v))/q/k]);
      v=permute(v,[2 1 3]);

      nv = prod(size(v));
      n1 = size(v,1);

      switch aggregation_mode
        case 'mean'
          ycon(count==k,:)=reshape(mean(v),q,nv/n1/q).';
        case 'sum'
          ycon(count==k,:)=reshape(sum(v),q,nv/n1/q).';
        case 'min'
          ycon(count==k,:)=reshape(min(v),q,nv/n1/q).';
        case 'max'
          ycon(count==k,:)=reshape(max(v),q,nv/n1/q).';
        case 'std'
          ycon(count==k,:)=reshape(std(v),q,nv/n1/q).';
        case 'var'
          ycon(count==k,:)=reshape(var(v),q,nv/n1/q).';
        case 'median'
          ycon(count==k,:)=reshape(median(v),q,nv/n1/q).';
        case 'prod'
          ycon(count==k,:)=reshape(prod(v),q,nv/n1/q).';
        case 'harm'
          ycon(count==k,:)=reshape(1./mean(1./v),q,nv/n1/q).';
        case 'geom'
          ycon(count==k,:)=reshape(exp(mean(log(v))),q,nv/n1/q).';
%       case 'count'
%         % we already did count above
      end
    end
  end
end

% was it originally a character array?
if charflag
  xcon=char(xcon);
end

author: Sable

repo: mcbench-benchmarks


text:
function [err_norm,n] = InverseKinematics_LM_ver2(to, Target)
% Levenberg-Marquardt, Chan-Lawrence, Sugihara's modification
global uLINK

idx = FindRoute(to);
wn_pos = 1/0.3;
wn_ang = 1/(2*pi);
We = diag([wn_pos wn_pos wn_pos wn_ang wn_ang wn_ang]);
Wn = eye(length(idx));
We2 = We*We;

ForwardKinematics(1);
err = CalcVWerr(Target, uLINK(to));
Ek = err'*We*err;

for n = 1:10
  J  = CalcJacobian(idx);
  Jh = J'*We*J + Wn*(Ek + 0.002);  %Hk + wn

  gerr = J'*We*err;    %gk
  dq   = Jh \ gerr;    %new

  MoveJoints(idx, dq);
  ForwardKinematics(1);
  err = CalcVWerr(Target, uLINK(to));
  Ek2 = err'*We*err;
  if Ek2 < 1E-12
      break;
  elseif Ek2 < Ek
      Ek = Ek2;
  else
      MoveJoints(idx, -dq);  % revert
      ForwardKinematics(1);
      break,
  end
end

err_norm = norm(err);

author: s-kajita

repo: IntroductionToHumanoidRobotics


text:
function z=hyperg(a,b,c,x,n)
% HYPERGEOMETRIC2F1 Computes the hypergeometric function
% using a series expansion:
%
%    f(a,b;c;x)=
%
%    1 + [ab/1!c]x + [a(a+1)b(b+1)/2!c(c+1)]x^2 +
%    [a(a+1)(a+2)b(b+1)(b+2)/3!c(c+1)(c+2)]x^3 + ...
%
% The series is expanded to n terms
%
% This function solves the Gaussian Hypergeometric Differential Equation:
%
%     x(1-x)y'' + {c-(a+b+1)x}y' - aby = 0
%
% The Hypergeometric function converges only for:
% |x| < 1
% c != 0, -1, -2, -3, ...
%
%
% Comments to:
% Diego Garcia   - d.garcia@ieee.org
% Chuck Mongiovi - mongiovi@fast.net
% June 14, 2002

if nargin ~= 5
    error('Usage: hypergeometric2f1(a,b,c,x,n) --> Wrong number of arguments')
end

if (n <= 0 | n ~= floor(n))
    error('Usage: hypergeometric2f1(a,b,c,x,n) --> n has to be a positive integer')
end

% if (abs(x) > 1)
%    z=min(0.99,x);
%    return;
%    error('Usage: hypergeometric2f1(a,b,c,x,n) --> |x| has to be less than 1')
% end

if (c <= 0 & c == floor(c))
    error('Usage: hypergeometric2f1(a,b,c,x,n) --> c != 0, -1, -2, -3, ...')
end

z = 0;
m = 0;
while (m<n)
    if (m == 0)
        delta = 1;
    else
        delta = delta .* x .* (a + (m - 1)) .* (b + (m-1)) ./ m ./ (c + (m-1));
   end
   z = z + delta;
   m = m + 1;
end

author: jtkim-kaist

repo: Speech-enhancement


text:
function animate_results(results, record, folder_name)
import iris.thirdParty.polytopes.*;

if nargin < 2
  record = false;
else
  if nargin < 3
    folder_name = ['videos/', datestr(now,'yyyy-mm-dd_HH.MM.SS')];
  end
end

dim = length(results.start);
C = results.e_history{1}.C;
d = results.e_history{1}.d;
A_bounds = results.p_history{1}.A(end-2*dim+1:end,:);
b_bounds = results.p_history{1}.b(end-2*dim+1:end);
bound_pts = lcon2vert(A_bounds, b_bounds);
lb = min(bound_pts)';
ub = max(bound_pts)';

if record
  frame = 1;
  system(sprintf('mkdir -p %s', folder_name));
  w = VideoWriter([folder_name, '/', 'animation']);
  w.FrameRate = 5;
  w.open();
  save([folder_name, '/', 'results'], 'results');

  delay = 1/w.FrameRate;
  delay1 = 1;
  delay2 = 2;
  gif_fname = [folder_name, '/', 'animation.gif'];
  loops = 65525;

  draw([], [], C, d, results.obstacles, lb, ub, results);
  h = gcf;
  movegui(h);
  w.writeVideo(getframe(h));
  print(gcf, sprintf('%s/%d_a', folder_name, 0), '-dpdf');
  img = getframe();
  [M, c_map]= rgb2ind(img.cdata,256);
  imwrite(M,c_map,[gif_fname],'gif','LoopCount',loops,'DelayTime',delay1);
end

figure(2);
clf;
for j = 1:length(results.p_history)
  A = results.p_history{j}.A;
  b = results.p_history{j}.b;
  draw(A, b, C, d, results.obstacles, lb, ub, results);
  if record
    h = gcf;
    movegui(h);
    w.writeVideo(getframe(h));
    print(gcf, sprintf('%s/%d_a', folder_name, j), '-dpdf');
    img = getframe();
    [M, c_map]= rgb2ind(img.cdata,256);
    imwrite(M,c_map,[gif_fname],'gif','WriteMode','append','DelayTime',delay)
  end
%   pause(0.1);
  C = results.e_history{j+1}.C;
  d = results.e_history{j+1}.d;
  draw(A, b, C, d, results.obstacles, lb, ub, results);
  if record
    h = gcf;
    movegui(h);
    w.writeVideo(getframe(h));
    print(gcf, sprintf('%s/%d_b', folder_name, j), '-dpdf', '-painters');
    img = getframe();
    [M, c_map]= rgb2ind(img.cdata,256);
    if j==length(results.p_history)
      imwrite(M,c_map,[gif_fname],'gif','WriteMode','append','DelayTime',delay2);
    else
      imwrite(M,c_map,[gif_fname],'gif','WriteMode','append','DelayTime',delay);
    end
  end
%   pause(0.25);
end

if record
  w.close();
end
end

function draw(A, b, C, d, obstacles, lb, ub, results)
  import iris.drawing.*;
  dim = length(results.start);
  if dim == 2
      draw_2d(A,b,C,d,obstacles,lb,ub);
      plot(results.start(1), results.start(2), 'go', 'MarkerFaceColor', 'g', 'MarkerSize', 8);
  elseif dim == 3
      draw_3d(A,b,C,d,obstacles,lb,ub);
      plot3(results.start(1), results.start(2), results.start(3), 'go', ...
        'MarkerFaceColor', 'g', 'MarkerSize', 15);
  else
    draw_Nd(A,b,C,d,obstacles,lb,ub);
  end
end

author: rdeits

repo: iris-distro

save_path: github-repos/MATLAB/rdeits-iris-distro


text:
function [ dist ] = dist_latlong( lat1, long1, lat2, long2, ref_lat, ref_long )
%DIST_LATLONG calculates distance between two specified points in meters
%accurate only in the (wide) area around the geodetic reference point ref_lat/_long

    [x1, y1] = latlong2xy( lat1, long1, ref_lat, ref_long );
    [x2, y2] = latlong2xy( lat2, long2, ref_lat, ref_long );

    dist = 1000 * sqrt( (x1-x2)^2 + (y1-y2)^2 );
end

repo: tdoa-evaluation-rtlsdr


text:
function [boxes, scores] = run_gop(im, num_proposals)

  if ~isdeployed
    old_path = path;
    addpath('/path/to/geodesic/gop_matlab/matlab');
  end

  interp_data = load('/path/to/gop_num_candidates_interpolation_data.mat');
  interp_data.num_candidates = [-1; interp_data.num_candidates; 1e10];
  interp_data.all_params = [interp_data.all_params(1,:); interp_data.all_params; interp_data.all_params(end,:)];
  params = interp1(interp_data.num_candidates, interp_data.all_params, num_proposals);
  assert(~any(isnan(params)));
  N_S = params(1);
  N_T = params(2);

  boxes = gop_wrapper(im, N_S, N_T);
  scores = [];

  if ~isdeployed
    path(old_path);
  end
end

function boxes = gop_wrapper(im, N_S, N_T)
% reimplementation of eval_box.py
  gop_mex( 'setDetector', 'MultiScaleStructuredForest("/path/to/geodesic/gop_matlab/data/sf.dat")' );
  data_path = '/path/to/geodesic/gop_matlab/data/';
  max_iou = 0.7;
  args = {};
  if N_S > 3
    args = [args {'seed', [data_path 'seed_final.dat']}];
  end
  if N_T >= 5
    max_iou = 0.8;
  end
  if N_T > 6
    max_iou = 0.9;
  end
  if N_S > 40
    args = [args {'unary', N_S, N_T, 'seedUnary()', 'backgroundUnary({0,1,15})', ...
                  'unary', 0, N_T,   'seedUnary()', 'backgroundUnary({0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15})', 0.1, 1}];
  elseif N_S > 1
    args = [args {'unary', N_S, N_T, 'seedUnary()', 'backgroundUnary({0,15})'}];
  else
    args = [args {'unary', N_S, N_T, 'seedUnary()', 'backgroundUnary({15})'}];
  end

  args = [{'max_iou', max_iou}, args];
  p = Proposal(args{:});
  tic();
  os = OverSegmentation( im );
  props = p.propose( os );
  boxes = double(os.maskToBox( props ));
  toc();
end

author: hosang

repo: detection-proposals


text:
clear all
close all

disp('Synthetic data from nonlinear ODE model');
disp('defined in Ramsay et al. (2007)');
disp('based on Van der Pol oscillator and which');
disp('reduces to Fitzhugh-Nagumo for certain parameters');

sigma_e=0.01;
[M,U] = mci_ramsay_struct(sigma_e);

% Use higher integration tolerances than by default
M.reltol=1e-3;
M.abstol=1e-5;

rand_init=0;
if rand_init
    P=spm_normrnd(M.vpE,M.pC,1);
else
    P=[log(0.2) log(0.2)]';
end

tic;
Y = mci_ramsay_gen(P,M,U);
toc

mci_plot_outputs(M,Y);

%mcmc.inference='vl';
%mcmc.inference='langevin';
%mcmc.maxits=16;
%mcmc.verbose=1;

mcmc.inference='ais';
mcmc.anneal='geometric';
mcmc.prop='lmc';
mcmc.nprop=1;
mcmc.J=16;
mcmc.maxits=8;
mcmc.rec_traj=1;

post = spm_mci_post (mcmc,M,U,Y,P);

if ~strcmp(mcmc.inference,'vl')
    spm_mci_diag(post);

    stats = spm_mci_mvnpost (post,'ESS')
    stats = spm_mci_mvnpost (post,'thinning')

    for j=1:length(P),
        spm_mci_quantiles (post,j,0);
    end
end

load ramsay-surface
figure
surf(S.x,S.y,L);
xlabel(S.name{1});
ylabel(S.name{2});

figure;
imagesc(S.pxy(1,:),S.pxy(2,:),L);
axis xy
hold on
xlabel(S.name{1});
ylabel(S.name{2});
hold on
ms=10;
plot(M.pE(1),M.pE(2),'wo','MarkerSize',ms);
plot(P(1),P(2),'w+','MarkerSize',ms);
if ~strcmp(mcmc.inference,'vl')
    j=post.ind;
    plot(post.P(1,j),post.P(2,j),'w.');
end
plot(post.Ep(1),post.Ep(2),'kx','MarkerSize',ms);

if strcmp(mcmc.inference,'ais')
    % plot individual trajectories
    figure
    xlim([-2 0]);
    ylim([-2 0]);
    hold on
    for i=1:mcmc.maxits,
        xx=squeeze(post.traj(i,1,:));
        yy=squeeze(post.traj(i,2,:));
        plot(xx,yy,'k-');
        plot(xx(end),yy(end),'kx','MarkerSize',ms);
        xlabel(S.name{1});
        ylabel(S.name{2});
        disp('Press space to see more trajectories');
        pause

    end
end


text:
function pool_params = load_pooling_params(pool_params_def, varargin)
% load_pooling_params(pool_params_def,...): it initializes the struct that
% contains the adaptive region max pooling parameters and the parameters
% related to the region type
%
% INPUTS:
% 1) pool_params_def: string with path to the configuration file that
% contains the adaptive region max pooling parameters
% The rest input arguments are given in the form of Name,Value pair
% arguments and are related to the region type:
% 'scale_inner': scalar value with the scaling factor of the inner rectangle
% of the region. In case this value is 0 then actually no inner rectangle
% is being used
% 'scale_outer': scalar value with the scaling factor of the outer rectangle
% of the region.
% 'half_bbox': intiger value in the range [1,2,3,4]. If this parameter is set
% to 1, 2, 3, or 4 then each bounding box will be reshaped to its left,
% right, top, or bottom half part correspondingly. This action is performed
% prior to scaling the box according to the scale_inner and scale_outer
% params. If this parameter is missing or if it is empty then the action of
% taking the half part of bounding box is NOT performed.
%
% OUTPUT:
% 1) pool_params: struct that contains the adaptive region max pooling
% parameters and the parameters related to the region type
%
% This file is part of the code that implements the following ICCV2015 accepted paper:
% title: "Object detection via a multi-region & semantic segmentation-aware CNN model"
% authors: Spyros Gidaris, Nikos Komodakis
% institution: Universite Paris Est, Ecole des Ponts ParisTech
% Technical report: http://arxiv.org/abs/1505.01749
% code: https://github.com/gidariss/mrcnn-object-detection
%
% Part of the code in this file comes from the SPP-Net code:
% https://github.com/ShaoqingRen/SPP_net
%
% AUTORIGHTS
% --------------------------------------------------------
% Copyright (c) 2015 Spyros Gidaris
%
% "Object detection via a multi-region & semantic segmentation-aware CNN model"
% Technical report: http://arxiv.org/abs/1505.01749
% Licensed under The MIT License [see LICENSE for details]
% ---------------------------------------------------------
% Copyright (c) 2014, Shaoqing Ren
%
% This file is part of the SPP code and is available
% under the terms of the Simplified BSD License provided in
% LICENSE. Please retain this notice and LICENSE if you use
% this file (or any portion of it) in your project.
% ---------------------------------------------------------

ip = inputParser;
ip.addParamValue('scale_inner',   [], @isnumeric);
ip.addParamValue('scale_outer',   [], @isnumeric);
ip.addParamValue('half_bbox',     [], @isnumeric);

ip.parse(varargin{:});
opts = ip.Results;

%% Read the adaptive region max pooling parameters from the configuration file
[~, ~, ext] = fileparts(pool_params_def);
if isempty(ext), pool_params_def = [pool_params_def, '.m']; end
assert(exist(pool_params_def, 'file') ~= 0);

cur_dir = pwd; % change folder to avoid too long path for eval()
[pool_def_dir, pool_def_file] = fileparts(pool_params_def);

cd(pool_def_dir);
pool_params = eval(pool_def_file);
cd(cur_dir);

%% Set the parameters related to the region type
pool_params.scale_inner = opts.scale_inner;
pool_params.scale_outer = opts.scale_outer;
pool_params.half_bbox   = opts.half_bbox;
end

author: gidariss

repo: LocNet

save_path: github-repos/MATLAB/gidariss-LocNet


text:
function [logval val]=mndensity(x,mu,sigma,k)

% function [logval val]=mndensity(x,mu,sigma,k)
% computes the density of the multivariate normal distribution
% uses the log of the density (a.2.1), split into several parts for computation ease
% inputs:  - vector 'x': the argument of the density function
%          - vector 'mu': the mean of the distribution
%          - matrix 'sigma': the covariance of the distribution
%          - integer 'k': the dimension of the argument vector
% outputs: - scalar 'logval': the log of the density value
%          - scalar 'val': the density value

temp(1,1)=-(k/2)*log(2*pi);
temp(2,1)=-0.5*log(det(sigma));
temp(3,1)=-0.5*(x-mu)'*(sigma\(x-mu));

% temp(3,1)=-0.5*(x-mu)'*inv(sigma)*(x-mu);

logval=sum(temp);
val=exp(logval);


text:
function area = triangleArea3d(pt1, pt2, pt3)
%TRIANGLEAREA3D Area of a 3D triangle.
%
%   AREA = triangleArea3d(P1, P2, P3)
%   Computes area of the 3D triangle whose vertices are given by P1, P2 and
%   P3. Each vertex is either a 1-by-3 row vector, or an array with 3
%   columns, each column representing coordinate of a vertex.
%   The result AREA has as many rows as the number of rows of the largest
%   input array.
%   Compared to polygonArea3d, this function is assumed to be faster, as it
%   does not requires iteration over vertices. Moreover, it can be used to
%   computes the area of several triangles simultaneously.
%
%   AREA = triangleArea3d(PTS)
%   Concatenates vertex coordinates in a 3-by-3 array. Each row of the
%   array contains coordinates of one vertex.
%
%
%   Example
%   triangleArea3d([10 10 10], [30 10 10], [10 40 10])
%   ans =
%       300
%
%   See also
%   polygons3d, polygonArea3d
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2011-08-23,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

% if data is given as one array, split vertices
if nargin == 1
    pt2 = pt1(2,:);
    pt3 = pt1(3,:);
    pt1 = pt1(1,:);
end

% compute individual vectors
v12 = bsxfun(@minus, pt2, pt1);
v13 = bsxfun(@minus, pt3, pt1);

% compute area from cross product
area = vectorNorm3d(cross(v12, v13, 2)) / 2;

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%--- help for generic/report ---
%
%  Assigns the elements of interest to a rise_report.report object
%
%    - display the default inputs::
%
%         report(rise.empty(0))
%
%    - assign the reported elements in rep_items to destination_root::
%
%         report(obj,destination_root,rep_items)
%
%    - assign varargin to obj before doing the rest::
%
%         report(obj,destination_root,rep_items,varargin)
%
%  Args:
%
%     obj (rise | dsge):
%     destination_root (rise_report.report): handle for the actual report
%     rep_items (char | cellstr): list of desired items to report. This list
%        can only include
%
%        - endogenous
%        - exogenous
%        - observables
%        - parameters
%        - solution
%        - estimation
%        - estimation_statistics
%        - equations
%        - code
%
%  Returns:
%     :
%
%     - none (Update the report in place)
%
%

author: jmaih


function [idxa,idxb,flag] = gen_train_sample_kissme(label_train, cam_train)

% generate ground truth pairs for training
uni_label = unique(label_train);
idxa = []; % index of the first image in a pair
idxb = []; % index of the second image in a pair
flag = []; % indicate whether two images are of the same identity
for n = 1:length(uni_label)
    curr_label = uni_label(n);
    pos = find(label_train == uni_label(n));
    if length(pos) == 1
        pos = [pos, pos];
    end
    comb = nchoosek(pos,2);
    idxa = [idxa; comb(:, 1)];
    idxb = [idxb; comb(:, 2)];
end
% remove pairs from the same camera
cam1 = cam_train(idxa);
cam2 = cam_train(idxb);
Eq_pos = find(cam1 == cam2);
diff_pos = setdiff(1:length(idxa), Eq_pos);
idxa = idxa(diff_pos);
idxb = idxb(diff_pos);
nPos = length(idxa);
flag = [flag; ones(nPos, 1)];

% generate negative training pairs
nTrainImg = length(label_train);
rand_pos = ceil(rand(150000, 2).*nTrainImg);
ID1 = label_train(rand_pos(:, 1));
ID2 = label_train(rand_pos(:, 2));
Eq_pos = find(ID1 == ID2);
diff_pos = setdiff(1:150000, Eq_pos); % remove pairs of the same identity
rand_pos = rand_pos(diff_pos, :);
cam1 = cam_train(rand_pos(:, 1));
cam2 = cam_train(rand_pos(:, 2));
Eq_pos = find(cam1 == cam2);
diff_pos = setdiff(1:length(rand_pos), Eq_pos);% remove pairs of the same camera

%%%% training image pairs and their ground truth labels %%%%%%%%
idxa = [idxa; rand_pos(diff_pos(1:nPos), 1)];
idxb = [idxb; rand_pos(diff_pos(1:nPos), 2)];
flag = [flag; zeros(nPos, 1)];


text:
function [obj,objs] = dcc_inference_objective(parameters,data,dataAsym,m,l,n,univariate) %#ok<INUSL>
% Objective function used by DCC, CCC and related multivariate volatility models to perform
% inference.  It is not used (and cannot be used) to estimate parameters.
%
% USAGE:
%  [OBJ,OBJS] = dcc_inference_objective(PARAMETERS,DATA,DATAASYM,M,L,N,UNIVARIAT)
%
% INPUTS:
%   PARAMETERS   - Vector of ADCC parameters including possibly volatility and intercepts
%   DATA         - A T by K matrix of zero mean residuals -OR-
%                    K by K by T array of covariance estimators (e.g. realized covariance)
%   DATAASYM     - [OPTIONAL] K by K by T array of asymmetric covariance estimators only needed if
%                    DATA is 3-dimensional and O>0 or L>0
%   M            - Order of symmetric innovations in DCC model
%   L            - Order of asymmetric innovations in ADCC model
%   N            - Order of lagged correlation in DCC model
%   UNIVARIATE   - Cell array of structures containing information needed to compute volatilities.
%
% OUTPUTS:
%   OBJ          - "As if" objective for estimating correlation intercept parameters
%   OBJS         - T by 1 vector of "as if" objectives for estimating correlation intercept parameters
%
% COMMENTS:
%
% See also DCC, CCC_MVGARCH

% Copyright: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 1    Date: 4/13/2012

[k,~,T] = size(data);
% Parse parameters
count = 0;
for i=1:k
    u = univariate{i};
    count = count + u.p+u.o+u.q+1;
end
garchParameters = parameters(1:count);
offset = count;
% R is next
count = k*(k-1)/2;
R = corr_ivech(parameters(offset + (1:count)));
offset = offset + count;
if l>0
    count = k*(k+1)/2;
    N = ivech(parameters(offset+(1:count)));
end

H = dcc_reconstruct_variance(garchParameters,univariate);
stdData = zeros(k,k,T);
stdDataAsym = zeros(k,k,T);
for t=1:T
    h = sqrt(H(t,:));
    stdData(:,:,t) = data(:,:,t)./(h'*h);
    stdDataAsym(:,:,t) = dataAsym(:,:,t)./(h'*h);
end

scales = diag(mean(stdData,3));
objs = zeros(T,1);
for j=1:k-1 % Cols
    for i=j+1:k % Rows
        scale = sqrt(scales(i)*scales(j));
        errors = squeeze(stdData(i,j,:))/scale - R(i,j);
        objs = objs + 0.5*(errors.^2);
    end
end

if l>0
    for j=1:k
        for i= j:k
            errors = squeeze(stdDataAsym(i,j,:)) - N(i,j);
            objs = objs + 0.5*(errors.^2);
        end
    end
end

obj = sum(objs);

author: bashtage

repo: mfe-toolbox


text:
function uI = edgeinterpolate2(exactu,node,edge,face,face2edge)
%% EDGEINTERPOLATE2 interpolate to the quadratic (1st type) edge finite element space.
%
% uI = edgeinterpolate(u,node,edge) interpolates a given function
% u into the lowesr order edge finite element spaces. The coefficient
% is given by the line integral int_e u*t ds. Simpson rule is used to
% evaluate this line integral.
%
% The input u could be a funtional handel or an array of length N (linear
% element) or N+NE (quadratic element).
%
% Example
%
%   node = [-1,-1,-1; 1,-1,-1; 1,1,-1; -1,1,-1; -1,-1,1; 1,-1,1; 1,1,1; -1,1,1];
%   elem = [1,2,3,7; 1,6,2,7; 1,5,6,7; 1,8,5,7; 1,4,8,7; 1,3,4,7];
%   maxIt = 3;
%   HcurlErr = zeros(maxIt,1);
%   L2Err = zeros(maxIt,1);
%   N = zeros(maxIt,1);
%   for k = 1:maxIt
%       [node,elem] = uniformbisect3(node,elem);
%       [elem2dof,T] = dof3NE2(elem);
%       pde = Maxwelldata2;
%       uI = edgeinterpolate2(pde.exactu,node,T.edge,T.face,T.face2edge);
%       HcurlErr(k) = getHcurlerror3ND2(node,elem,pde.curlu,uI);
%       L2Err(k) = getL2error3ND2(node,elem,pde.exactu,uI);
%       N(k) = length(uI);
%   end
%   figure(1)
%   r1 = showrate(N,HcurlErr,1,'r-+');
%   hold on
%   r2 = showrate(N,L2Err,1,'b-+');
%   legend('||u-u_I||_{curl}',['N^{' num2str(r1) '}'],...
%          '||u-u_I||',['N^{' num2str(r2) '}'],'LOCATION','Best');
%
% See also edgeinterpolate, edgeinterpolate1
%
% Maxwell doc Section: Dirichlet
% boundary condition
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

%% Evaluate function at vertices and middle points
if isnumeric(exactu)
    uQ = exactu;
else
    mid = (node(edge(:,1),:) + node(edge(:,2),:))/2;
    uQ = exactu([node; mid]);
end

%% Edge dof coefficients
N = size(node,1);	NE = size(edge,1);    NF = size(face,1);
edgeVec = node(edge(:,2),:)-node(edge(:,1),:);
uI = zeros(2*(NE+NF),1);
uI(1:NE,1) = dot(edgeVec,(uQ(edge(:,1),:)+uQ(edge(:,2),:)+4*uQ(N+1:N+NE,:))/6,2);
uI(NE+1:2*NE,1) = dot(edgeVec,0.5*(uQ(edge(:,1),:)-uQ(edge(:,2),:)),2);

%% Face dof coefficients
eik = node(face(:,3),:)-node(face(:,1),:);
eij = node(face(:,2),:)-node(face(:,1),:);
uquadpts = uQ(N+face2edge(:,1),:)+uQ(N+face2edge(:,2),:)+uQ(N+face2edge(:,3),:);
lf = [4*dot(eik,uquadpts,2) ...
      4*dot(eij,uquadpts,2)];
face2edgeDofValue = [uI(face2edge(:,1:3)) uI(face2edge(:,1:3)+NE)];
localMatrix = [4 8 4 -4  0 4; ...
               8 4 -4 0 -4 4]';
lf = lf - face2edgeDofValue*localMatrix;
uI((2*NE+1):end) = [2*lf(:,1) - lf(:,2); ...
                    2*lf(:,2) - lf(:,1)]/3;

repo: ifem


text:
%% Example: Dense SDP relaxation for quartic optimization over the sphere

clc; clear; close all; restoredefaultpath; % start clean

mosekpath = '../../mosek';
pgdpath   = '../STRIDE';
sdpnalpath  = '../../SDPNAL+v1.0';
manoptpath  = '../manopt';
utilspath   = '../utils';
addpath(genpath(utilspath))
addpath(genpath(manoptpath))
addpath(genpath('./solvers'))
addpath(genpath('../spotless')) % Use spotless for defining polynomials
addpath('../SDPRelaxations') % implementations for SDP relaxation

%% Generate random binary quadratic program
d       = 30; % d variables
x       = msspoly('x',d); % symbolic decision variables using SPOTLESS
v       = monomials(x,0:4);
c       = randn(length(v),1);
f       = c'*v; % objective function
h       = sum(x.^2) - 1; % equality constraints

%% Relax BQP into an SDP
problem.vars            = x;
problem.objective       = f;
problem.equality        = h;
kappa                   = 2; % relaxation order
[SDP,info]              = dense_sdp_relax(problem,kappa);
SDP.M       = 3; % upper bound on the trace of the moment matrix
% need the following for fast computation in the local search method
info.v      = msspoly2degcoeff(info.v);
info.f      = msspoly2degcoeff(info.f);
info.J      = msspoly2degcoeff(info.J);

%{
%% Solve using MOSEK, should be slow
prob       = convert_sedumi2mosek(SDP.sedumi.At,...
                                  SDP.sedumi.b,...
                                  SDP.sedumi.c,...
                                  SDP.sedumi.K);
addpath(genpath(mosekpath))
[~,res]    = mosekopt('minimize info',prob);
[Xopt,yopt,Sopt,obj] = recover_mosek_sol_blk(res,SDP.blk);
rmpath(genpath(mosekpath))

figure; bar(eig(Xopt{1})); % if rank = 1, then relaxation is exact/tight
%}

%% solve using stride
addpath(genpath(pgdpath))

pgdopts.pgdStepSize     = 10;
pgdopts.SDPNALpath      = sdpnalpath;
pgdopts.tolADMM         = 10e-5;
pgdopts.phase1          = 1;
pgdopts.rrOpt           = 1:3;
pgdopts.rrFunName       = 'local_search_q4s';
pgdopts.rrPar           = info;
pgdopts.maxiterLBFGS    = 1000;
pgdopts.maxiterSGS      = 300;
pgdopts.tolLBFGS        = 1e-12;
pgdopts.tolPGD          = 1e-8;

[outPGD,Xopt,yopt,Sopt]     = PGDSDP(SDP.blk, SDP.At, SDP.b, SDP.C, [], pgdopts);
time_pgd                    = outPGD.totaltime;
% round solutions and check optimality certificate
res = get_performance_q4s(Xopt,yopt,Sopt,SDP,info,pgdpath);

%% helper functions
function s = msspoly2degcoeff(f)
[~,degmat,coeff,~] = decomp(f);
s.degmat = degmat';
s.coefficient = coeff;
end

author: MIT-SPARK

repo: CertifiablyRobustPerception


text:
% saves images of individual estimated components....
function varout = check_profiles(sep5, offset5, pathname, prename, niter)
p.path = pathname;

for kk5 = 1: length(offset5)
    for ll5=1 : length(sep5)
        p.namedir = [prename num2str(100*sep5(ll5)) 'offset_' num2str(offset5(kk5))];
        cd ([p.path p.namedir])
        load ([p.namedir '-iter_1'])
        %         ims(psf);
        %         SaveImageFULL('psf', 'pf');
        is = imarfun('imsum',dpixc_ind);
        prof = double(is(:,15));
        varout(kk5, ll5,:) = prof;
% % %         for mm5=1:niter
% % %             varout = bg;
% % %             imstiled(icapixNMF{mm5}(:,:,1:2));
% % %             close all
% % %         end
    end
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [dat,stats] = hewma_extract_voxel(EXPT,coords)
% Extracts data for one voxel and runs hewma2 on it.
%
% :Usage:
% ::
%
%     function [dat,stats] = hewma_extract_voxel(EXPT,coords)
%
% Uses plot option in hewma2 to create plots.

dat = [];
stats = [];

persistent timeseriesfile
persistent lamfile
persistent files
persistent vfiles

zsl = num2str(coords(3));

mm_coords = spm_orthviews('Pos')';
fprintf(1, 'mm coordinates are (x, y, z) = (%3.0f, %3.0f, %3.0f)\t\t Voxel coords: (%3.0f, %3.0f, %3.0f)\n', mm_coords, coords);

if isempty(files) || isempty(vfiles)
    mydir = spm_get(-1,'*','Select Directory above individual subject directories',pwd);
end

if isempty(files)

    % EXPT = getfunctnames2(EXPT,['z_slice' zsl '.mat'],'tmp');
    EXPT = getfunctnames2(EXPT,['z_slice' zsl '.mat'], 'tmp', [], mydir);

    files = str2mat(EXPT.tmp{:});
end

if isempty(vfiles)
    %EXPT = getfunctnames2(EXPT,['var_slice' zsl '.mat'],'tmp');
    EXPT = getfunctnames2(EXPT,['var_slice' zsl '.mat'], 'tmp', [], mydir);

    vfiles = str2mat(EXPT.tmp{:});
end

if isempty(files) || all(files(:) == ' ')
    disp('Cannot find z_slice files for subjects.  Check that path names are correct in EXPT structure.');
    return
end

if isempty(vfiles) || all(vfiles(:) == ' ')
    disp('Cannot find var_slice files for subjects.  Check that path names are correct in EXPT structure.');
    return
end

if isempty(timeseriesfile)

    timeseriesfile = spm_get(1,'*mat','Select hewma_timeseries.mat containing img dims',pwd);

end

  load(timeseriesfile, 'xdim', 'ydim');
% % try
% %     load hewma0001/hewma_timeseries xdim ydim
% % catch
% %     if ~(exist('hewma_timeseries') == 2)
% %             file = spm_get(1,'*mat','Select hewma_timeseries.mat',pwd);
% %             [dd,ff,ee] = fileparts(file);
% %             cd(dd)
% %             load hewma_timeseries xdim ydim
% %     end
% % end

if isempty(lamfile)
    lamfile = spm_get(1,'*mat','Select lam.mat file containing lambda parameter for one subject',pwd);
end

load(lamfile, 'lam');

% % try
% %     load([EXPT.subjects{1} filesep 'lam'])    % load lambda param, assume same for all subjects
% % catch
% %     error('Cannot find lam.mat.  Start in dir above individual ewma results dirs.')
% % end

fprintf(1,'Loading data, subject ');

for i = 1:length(EXPT.subjects)
    fprintf(1,'%3.0f ',i);

    % load EWMA stat
    load(deblank(files(i,:)));
    zdat = full(zdat);
    T = size(zdat,2);
    zdat = reshape(zdat,xdim,ydim,T);

    try
        dat(i,:) = squeeze(zdat(coords(1),coords(2),:));
    catch
        disp('Warning!  data for this subject does not match others in size.  ')
        dat(i,:) = squeeze(zdat(coords(1),coords(2),1:size(dat,2)));
    end

    % load variance
    load(deblank(vfiles(i,:)));
    vardat = full(vardat);
    vardat = reshape(vardat,xdim,ydim,T);

    try
        vdat(i,:) = squeeze(vardat(coords(1),coords(2),:));
    catch
        vdat(i,:) = squeeze(vardat(coords(1),coords(2),1:size(vdat,2)));
    end

end
fprintf(1,'\n');

if isfield(EXPT,'cov');
    mycov = EXPT.cov;
end
if ~isempty(mycov), mycov = mycov(:,1); end

[p,tm,Zcor,sbmean,Zpop,tvals,sb,stats] = hewma2(dat, vdat, lam,1,1,mycov);

%figure;plot(dat');
%hold on; plot(mean(dat),'k','LineWidth',2);
%title('EWMA statistics with mean');

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [e, g] = gp_eg(w, gp, x, y, varargin)
%GP_EG  Evaluate the energy function (un-normalized negative marginal
%       log posterior) and its gradient
%
%  Description
%    [E, G] = GP_EG(W, GP, X, Y, OPTIONS) takes a Gaussian process
%    structure GP together with a matrix X of input vectors and a
%    matrix Y of targets, and evaluates the energy function E and
%    its gradient G. Each row of X corresponds to one input vector
%    and each row of Y corresponds to one target vector.
%
%    The energy is minus log posterior cost function:
%        E = EDATA + EPRIOR
%          = - log p(Y|X, th) - log p(th),
%    where th represents the parameters (lengthScale,
%    magnSigma2...), X is inputs and Y is observations (regression)
%    or latent values (non-Gaussian likelihood).
%
%    OPTIONS is optional parameter-value pair
%      z - optional observed quantity in triplet (x_i,y_i,z_i)
%          Some likelihoods may use this. For example, in case of
%          Poisson likelihood we have z_i=E_i, that is, expected
%          value for ith case.
%
%  See also
%    GP_E, GP_G
%

% Copyright (c) 2010 Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

% Single function for some optimization routines, no need for mydeal...
e=gp_e(w, gp, x, y, varargin{:});
if nargout>1
  if isnan(e)
    g=NaN;
  else
    g=gp_g(w, gp, x, y, varargin{:});
  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% LOTKA-VOLTERRA system

clear all, close all, clc
figpath = '../FIGURES/NOT_USED/'; mkdir(figpath)
datapath = '../DATA/';
addpath('../utils');

%% Load Models
InputSignalTypeModel = 'sphs'; % prbs; chirp; noise; sine2; sphs; mixed

% DelayDMDc
ModelTypeDMDc = 'DMDc';
load(fullfile(datapath,['EX_LOTKA_SI_',ModelTypeDMDc,'_',InputSignalTypeModel,'.mat'])) % Model: DelayDMDc , DMDc
Models.DelayDMDc = Model;

% SINDYc
load(fullfile(datapath,['EX_LOTKA_SI_SINDYc','_',InputSignalTypeModel,'.mat'])) % Model
Models.SINDYc = Model;

% NARX
load(fullfile(datapath,['EX_LOTKA_SI_NARX','_',InputSignalTypeModel,'.mat'])) % Model
Models.NARX = Model;

%% Get training data
InputSignalType = InputSignalTypeModel;% prbs; chirp; noise; sine2; sphs; mixed
Ndelay = Models.DelayDMDc.Ndelay;
ONLY_TRAINING_LENGTH = 1;
getTrainingData

%% Get validation data
InputSignalType = 'sine2';% prbs; chirp; noise; sine2; sphs; mixed
Ndelay = Models.DelayDMDc.Ndelay;
getValidationData

%% FIGURE 1:  // Model comparison + Validation

% Reference
[tA,xA] = ode45(@(t,x)lotkacontrol(t,x,forcing(x,t),a,b,d,g),tspanV,x(end,1:2),options);   % true model

% DelayDMDc / DMDc
% Model
if Ndelay == 1
    x0      = [x(end,1:2)];
    Hunew   = [u(end),u_valid(1:end-1)];
    [xB,tB] = lsim(Models.DelayDMDc.sys,Hunew,tspanV,[x0-[xmean]]');
elseif Ndelay > 1
    x0      = [x(end-Ndelay+1,1:2),x(end,1:2)];
    Hunew   = [ u(end-Ndelay+1:end),u_valid(1:end-Ndelay);
        u(end),u_valid(1:end-1)];
    [xB,tB] = lsim(Models.DelayDMDc.sys,Hunew,tspanV,[x0-[xmean]]');
    xB = xB(:,3:4); xB = xB + repmat(xmean,[size(xB,1) 1]);
end
xB = xB + repmat(xmean,[length(tB) 1]);

% SINDYc
[tC,xC]=ode45(@(t,x)sparseGalerkinControl(t,x,forcing(x,t),Models.SINDYc.Xi(:,1:2),Models.SINDYc.polyorder,Models.SINDYc.usesine),tspanV,x(end,1:2),options);  % approximate

% NARX
Udummy = [u(end),u_valid(1:end)];
Hdummy = zeros(2,size(Udummy,2));
Hdummy(:,1) = x(end,1:2)';
[Us,Ui,Si] = preparets(Models.NARX.net,con2seq(Udummy),{},con2seq(Hdummy));
xD = Models.NARX.net(Us,Ui,Si); % Predict on validation data
xD = cell2mat(xD)'; xD = [x0;xD(1:end-1,:)];

%%
clear ph
h = figure;
subplot(2,1,1), box on, hold on
plot([tB(1) tB(1)],[0 150],'--k')
plot(t,x(:,1),'Color',[.4 .4 .4],'LineWidth',1.5);
plot(tA,xA(:,1),'k','LineWidth',1.5);
plot(tB,xB(:,1),'r-','LineWidth',1.5);
plot(tC,xC(:,1),'g--','LineWidth',1.5);
plot(tA(1:end),xD(:,1),'-.','Color',[0.7,0.7,1],'LineWidth',1.5);
grid on
ylim([0 150])
% xlabel('Time','FontSize',13)
ylabel('Prey, x_1','FontSize',13)
set(gca,'FontSize',13)
subplot(2,1,2), box on, hold on
plot([tB(1) tB(1)],[0 60],'--k')
ph(1) = plot(t,x(:,2),'Color',[.4 .4 .4],'LineWidth',1.5);
ph(2) = plot(tA,xA(:,2),'k','LineWidth',1.5);
ph(3) = plot(tB,xB(:,2),'r-','LineWidth',1.5);
ph(4) = plot(tC,xC(:,2),'g--','LineWidth',1.5);
% ph(5) = plot(tB,xD(:,2),'-.','Color',[0.7,0.7,1],'LineWidth',1.5);
ph(5) = plot(tA(1:end),xD(:,2),'-.','Color',[0.7,0.7,1],'LineWidth',1.5);
l1=legend(ph,'Training','Validation',ModelTypeDMDc,'SINDYc','NARX');
set(l1,'Location','NorthWest')
grid on
ylim([0 60])
ylabel('Predator, x_2','FontSize',13)
set(gca,'FontSize',13)
xlabel('Time','FontSize',13)
set(gca,'FontSize',13)

set(h,'Units','Inches');
set(gcf,'Position',[1 1 6. 5.5])
pos = get(h,'Position');
set(h,'PaperPositionMode','Auto','PaperSize',[pos(3), pos(4)])
print(h,'-painters','-depsc2',  '-loose','-cmyk', [figpath,'EX_LOTKA_SI_Comparison_Validation_',InputSignalTypeModel,'_',InputSignalType,'.eps'],'-r0');

%% Actuation signal
clear ph
figure,box on, hold on,
ccolors = get(gca,'colororder');
plot([tA(1),tA(1)],[-15 260],':','Color',[0.4,0.4,0.4],'LineWidth',1.5)
plot(t,u,'-k','LineWidth',1);
plot(t_valid,u_valid,'-k','LineWidth',1);
grid off
ylim([min([u u_valid])+0.05*min([u u_valid]) max([u u_valid])+0.05*max([u u_valid])])
xlim([0 200])
xlabel('Time')
ylabel('Input')
set(gca,'LineWidth',1, 'FontSize',14)
set(gcf,'Position',[100 100 300 200])
set(gcf,'PaperPositionMode','auto')
print('-depsc2', '-loose','-cmyk', [figpath,'EX_LOTKA_SI_Comparison_Validation_',InputSignalTypeModel,'_',InputSignalType,'_Actuation.eps']);

%% MPC
% Parameters MPC
options = optimoptions('fmincon','Algorithm','sqp','Display','none', ...
    'MaxIterations',100);
N  = 10;                        % Control / prediction horizon (number of iterations)
Duration = 100;                 % Run for 'Duration' time units
Ton = 0;                        % Time units when control is turned on
Nvar = 2;
getMPCparams

x0n=xA(end,:)';                 % Initial condition
xref1 = [g/d;a/b];              % critical point
Tcontrol = tA(end);             % Time offset to combine training, prediction and control phase

% Parameters Models
ModelCollection = {ModelTypeDMDc, 'SINDYc', 'NARX'};
Nmodels = length(ModelCollection);

% Parameters True Model
p.a = a;
p.b = b;
p.d = d;
p.g = g;

%% Load Results
clear Results
load(fullfile(datapath,['MPC_Results_AllModels_',InputSignalTypeModel,'_',InputSignalType,'.mat']),'Results')

%% Unforced
[tU,xU] = ode45(@(t,x)lotkacontrol(t,x,0,a,b,d,g),Results(1).t,xA(end,1:2),options);   % true model

%% Show results
MODEL_COLOR = 1;
Ts = Models.DelayDMDc.dt;
ct = (Duration/Ts);
for iM = 1:Nmodels
    select_model = ModelCollection{iM};
    xHistory = Results(iM).x;
    uHistory = Results(iM).u;
    tHistory = Results(iM).t;
    VIZ_SI_Validation_MPC
end

%% Get cost // evaluate performance
% Data2 = Datal
Costs = zeros(length(Results(1).u),3,Nmodels);

for iM = 1:Nmodels
    [J] = evalObjectiveFCN(Results(iM).u,Results(iM).x,Results(iM).xref,diag(Q),R,Ru);
    Costs(:,1,iM) = J; %Costs(:,2,iM) = Js; Costs(:,3,iM) = Ju;
end

%% Show cost
clear ph
figure, hold on, box on
ph(1) = plot(Results(1).t,cumsum(Results(1).J),'-r','LineWidth',2);
ph(2) = plot(Results(1).t,cumsum(Results(2).J),'--g','LineWidth',2);
ph(3) = plot(Results(1).t,cumsum(Results(3).J),'-.','Color',[0.7,0.7,1],'LineWidth',2);

xlim([200 xlimval(2)])
l1=legend(ph,ModelCollection);
set(l1,'Location','SouthEast')
% ylim([0 5.5*10^4])
ylabel('Cost','FontSize',14)
xlabel('Time','FontSize',14)
set(gca,'LineWidth',1, 'FontSize',14)
set(gcf,'Position',[100 100 300 200])
set(gcf,'PaperPositionMode','auto')
print('-depsc2', '-loose','-cmyk', [figpath,'EX_LOTKA_SI_Comparison_Validation_',InputSignalTypeModel,'_',InputSignalType,'_Cost.eps']);

clear ph
figure, hold on, box on
ph(1) = plot(Results(1).t,cumsum(Costs(:,1,1)),'-r','LineWidth',2);
ph(2) = plot(Results(1).t,cumsum(Costs(:,1,2)),'--g','LineWidth',2);
ph(3) = plot(Results(1).t,cumsum(Costs(:,1,3)),'-.','Color',[0.7,0.7,1],'LineWidth',2);

xlim([200 xlimval(2)])
l1=legend(ph,ModelCollection);
set(l1,'Location','SouthEast')
% ylim([0 5.5*10^4])
ylabel('Cost','FontSize',14)
xlabel('Time','FontSize',14)
set(gca,'LineWidth',1, 'FontSize',14)
set(gcf,'Position',[100 100 300 200])
set(gcf,'PaperPositionMode','auto')
%%
iJ = 3;
clear ph
figure, hold on, box on
ph(1) = plot(Results(1).t,cumsum(Costs(:,iJ,1)),'-r','LineWidth',2);
ph(2) = plot(Results(1).t,cumsum(Costs(:,iJ,2)),'--g','LineWidth',2);
ph(3) = plot(Results(1).t,cumsum(Costs(:,iJ,3)),'-.','Color',[0.7,0.7,1],'LineWidth',2);

l1=legend(ph,ModelCollection);
set(l1,'Location','NorthEast')
ylabel('Cost','FontSize',14)
xlabel('Time','FontSize',14)
set(gca,'LineWidth',1, 'FontSize',14)
set(gcf,'Position',[100 100 300 200])
set(gcf,'PaperPositionMode','auto')

author: eurika-kaiser

repo: SINDY-MPC


text:
function J = jacobian( F )
%JACOBIAN   Jacobian determinant of a CHEBFUN2V.
%   J = JACOBIAN(F) computes the determinant of the Jacobian matrix associated
%   to the vector-valued CHEBFUN2V F. The CHEBFUN2V must have two components.
%
%   Note we return the determinant of the Jacobian matrix and not the Jacobian
%   matrix itself.
%
% See also CHEBFUN2/GRADIENT.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Empty check.
if ( isempty(F) )
    J = [];
    return
end

if ( F.nComponents == 3 )
    error('CHEBFUN:CHEBFUN2V:jacobian:notSquare', ...
        'Jacobian matrix is not square.')
end

% Determinant formula:
Fx = diff( F, 1, 2 );
Fy = diff( F, 1, 1 );

% Jacobian:
Fxc = Fx.components;
Fyc = Fy.components;
J = Fxc{1} .* Fyc{2} - Fyc{1} .* Fxc{2};

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% CONFIGROBOT initializes parameters specific of a particular robot
%             (e.g., icubGazeboSim)

%% --- Initialization ---

Config.ON_GAZEBO         = true;
ROBOT_DOF                = 23;
Config.GRAV_ACC          = 9.81;

%%%% INITIAL CONDITIONS %%%%
Config.init_jointPos = [ 0.0,0.0,0.0, ...
                        -30.0,30.0,0.0,45.0, ...
                        -30.0,30.0,0.0,45.0, ...                   %
                         0.0,0.0,0.0,0.0,0.0,0.0, ...
                         0.0,0.0,0.0,0.0,0.0,0.0]*pi/180;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Robot configuration for WBToolbox
WBTConfigRobot           = WBToolbox.Configuration;
WBTConfigRobot.RobotName = 'icubSim';
WBTConfigRobot.UrdfFile  = 'model.urdf';
WBTConfigRobot.LocalName = 'WBT';

% Controlboards and joints list. Each joint is associated to the corresponding controlboard
WBTConfigRobot.ControlBoardsNames     = {'torso','left_arm','right_arm','left_leg','right_leg'};
WBTConfigRobot.ControlledJoints       = [];
Config.numOfJointsForEachControlboard = [];

ControlBoards                                        = struct();
ControlBoards.(WBTConfigRobot.ControlBoardsNames{1}) = {'torso_pitch','torso_roll','torso_yaw'};
ControlBoards.(WBTConfigRobot.ControlBoardsNames{2}) = {'l_shoulder_pitch','l_shoulder_roll','l_shoulder_yaw','l_elbow'};
ControlBoards.(WBTConfigRobot.ControlBoardsNames{3}) = {'r_shoulder_pitch','r_shoulder_roll','r_shoulder_yaw','r_elbow'};
ControlBoards.(WBTConfigRobot.ControlBoardsNames{4}) = {'l_hip_pitch','l_hip_roll','l_hip_yaw','l_knee','l_ankle_pitch','l_ankle_roll'};
ControlBoards.(WBTConfigRobot.ControlBoardsNames{5}) = {'r_hip_pitch','r_hip_roll','r_hip_yaw','r_knee','r_ankle_pitch','r_ankle_roll'};

for n = 1:length(WBTConfigRobot.ControlBoardsNames)

    WBTConfigRobot.ControlledJoints       = [WBTConfigRobot.ControlledJoints, ...
                                             ControlBoards.(WBTConfigRobot.ControlBoardsNames{n})];
    Config.numOfJointsForEachControlboard = [Config.numOfJointsForEachControlboard; length(ControlBoards.(WBTConfigRobot.ControlBoardsNames{n}))];
end

% Frames list
Frames.BASE       = 'root_link';
Frames.IMU        = 'imu_frame';
Frames.LEFT_FOOT  = 'l_sole';
Frames.RIGHT_FOOT = 'r_sole';
Frames.COM        = 'com';

% Config.SATURATE_TORQUE_DERIVATIVE: if true, the derivative of the control
% input is saturated. In this way, it is possible to reduce high frequency
% oscillations and discontinuities in the control input.
Config.SATURATE_TORQUE_DERIVATIVE         = false;

% if TRUE, the controller will STOP if the joints hit the joints limits
% and/or if the (unsigned) difference between two consecutive joints
% encoders measurements is greater than a given threshold.
Config.EMERGENCY_STOP_WITH_JOINTS_LIMITS  = false;
Config.EMERGENCY_STOP_WITH_ENCODER_SPIKES = false;

% Config.USE_MOTOR_REFLECTED_INERTIA: if set to true, motors reflected
% inertias are included in the system mass matrix. If
% Config.INCLUDE_COUPLING is true, then the coupling effect (some joints
% motion is the result of more than one motor motion) is taken into account.
% Config.INCLUDE_HARMONIC_DRIVE_INERTIA is true, then the harmonic drive
% reflected inertia is also considered
Config.USE_MOTOR_REFLECTED_INERTIA    = false;
Config.INCLUDE_COUPLING               = false;
Config.INCLUDE_HARMONIC_DRIVE_INERTIA = false;

% Config.USE_IMU4EST_BASE: if set to false, the base frame is estimated by
% assuming that either the left or the right foot stay stuck on the ground.
% Which foot the controller uses depends on the contact forces acting on it.
% If set to true, the base orientation is estimated by using the IMU, while
% the base position by assuming that the origin of either the right or the
% left foot do not move.
Config.USE_IMU4EST_BASE = false;

% Config.YAW_IMU_FILTER when the flag Config.USE_IMU4EST_BASE = true, then
% the orientation of the floating base is estimated as explained above. However,
% the foot is usually perpendicular to gravity when the robot stands on flat
% surfaces, and rotation about the gravity axis may be affected by the IMU drift
% in estimating this angle. Hence, when either of the flags Config.YAW_IMU_FILTER
% is set to true, then the yaw angle of the contact foot is ignored and kept
% equal to the initial value.
Config.FILTER_IMU_YAW   = true;

% Config.CORRECT_NECK_IMU: when set equal to true, the kinematics from the
% IMU and the contact foot is corrected by using the neck angles. If it set
% equal to false, recall that the neck is assumed to be in (0,0,0). Valid
% ONLY while using the ICUB HEAD IMU!
Config.CORRECT_NECK_IMU = false;

% Config.USE_QP_SOLVER: if set to true, a QP solver is used to account for
% inequality constraints of contact wrenches.
Config.USE_QP_SOLVER    = true;

% Ports name list
Ports.WRENCH_LEFT_FOOT  = '/wholeBodyDynamics/left_leg/cartesianEndEffectorWrench:o';
Ports.WRENCH_RIGHT_FOOT = '/wholeBodyDynamics/right_leg/cartesianEndEffectorWrench:o';
Ports.IMU               = ['/' WBTConfigRobot.RobotName '/inertial'];
Ports.NECK_POS          = ['/' WBTConfigRobot.RobotName '/head/state:o'];

% Ports dimensions
Ports.NECK_POS_PORT_SIZE         = 3;
Ports.IMU_PORT_SIZE              = 12;
Ports.IMU_PORT_ORIENTATION_INDEX = [1,2,3];
Ports.WRENCH_PORT_SIZE           = 6;

author: robotology

repo: whole-body-controllers


text:
%SETRNGSEED  Sets state of default random number generator
%
%     cv.setRNGSeed(seed)
%
% ## Input
% * __seed__ integer, new state for default random number generator.
%
% The function cv.setRNGSeed sets state of default random number generator to
% custom value.
%
% See also: rng, rand, randn
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function f=dctresample(f,L,dim)
%DCTRESAMPLE   Resample signal using Fourier interpolation
%   Usage:  h=dctresample(f,L);
%           h=dctresample(f,L,dim);
%
%   `dctresample(f,L)` returns a discrete cosine interpolation of the signal *f*
%   to length *L*. If the function is applied to a matrix, it will apply
%   to each column.
%
%   `dctresample(f,L,dim)` does the same along dimension *dim*.
%
%   If the input signal is not a periodic signal (or close to), this method
%   will give much better results than |fftresample| at the endpoints, as
%   this method assumes than the signal is even a the endpoints.
%
%   The algorithm uses a DCT type iii.
%
%   See also:  fftresample, middlepad, dctiii

%   AUTHOR: Peter L. Søndergaard

% ------- Checking of input --------------------
complainif_argnonotinrange(nargin,2,3,mfilename);

if nargin<3
  dim=[];
end;

[f,L,Ls,W,dim,permutedsize,order]=assert_sigreshape_pre(f,L,dim,'DCTRESAMPLE');

wasreal=isreal(f);

% The 'dim=1' below have been added to avoid dct and middlepad being
% smart about choosing the dimension.
f=dctiii(postpad(dctii(f,[],1),L))*sqrt(L/Ls);

f=assert_sigreshape_post(f,dim,permutedsize,order);

if wasreal
  f=real(f);
end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [Aeq beq]= getAbeq(n_seg, n_order, waypoints, ts, start_cond, end_cond)
    n_all_poly = n_seg*(n_order+1);% the number of all polynomial coefficients
    %#####################################################
    % p,v,a,j constraint in start,
    Aeq_start = zeros(4, n_all_poly);
    beq_start = zeros(4, 1);
    % STEP 2.1: write expression of Aeq_start and beq_start
    Aeq_start(1:1:4, 1:1:n_order+1) = [0, 0, 0, 0, 0, 0, 0, 1;
                                       0, 0, 0, 0, 0, 0, 1, 0;
                                       0, 0, 0, 0, 0, 2, 0, 0;
                                       0, 0, 0, 0, 6, 0, 0, 0];
    beq_start = start_cond';% p,v,a,j

    %#####################################################
    % p,v,a constraint in end
    Aeq_end = zeros(4, n_all_poly);
    beq_end = zeros(4, 1);
    % STEP 2.2: write expression of Aeq_end and beq_end
    T = ts(size(ts,1));% time of the last trajectory
    Aeq_end(1:1:4, n_all_poly-n_order:1:n_all_poly) = [ T^7,     T^6,    T^5,    T^4,   T^3, T^2, T, 1;
                                                      7*T^6,   6*T^5,  5*T^4,  4*T^3, 3*T^2, 2*T, 1, 0;
                                                     42*T^5,  30*T^4, 20*T^3, 12*T^2,   6*T,   2, 0, 0;
                                                    210*T^4, 120*T^3, 60*T^2,   24*T,     6,   0, 0, 0];
    beq_end = end_cond';% p,v,a,j

    %#####################################################
    % position constrain in all middle waypoints
    Aeq_wp = zeros(n_seg-1, n_all_poly);
    beq_wp = zeros(n_seg-1, 1);
    % STEP 2.3: write expression of Aeq_wp and beq_wp
    for midwp_index = 1:n_seg-1
        index = 1 + 8 * (midwp_index - 1);
        T = ts(midwp_index);
        Aeq_wp(midwp_index,index:index+7) = [T^7, T^6, T^5, T^4, T^3, T^2, T, 1];% the end of previous segment
        % Aeq_wp(midwp_index,index+8:index+8+7) = [0, 0, 0, 0, 0, 0, 0, 1];% the begin of next segment
    end
    beq_wp = waypoints(2:n_seg,1);

    %#####################################################
    % position continuity constrain between each 2 segments
    Aeq_con_p = zeros(n_seg-1, n_all_poly);
    beq_con_p = zeros(n_seg-1, 1);
    % STEP 2.4: write expression of Aeq_con_p and beq_con_p
    for con_p_index = 1:n_seg-1
        index = 1 + 8 * (con_p_index - 1);
        T = ts(con_p_index);
        Aeq_con_p(con_p_index,index:index+7) = [T^7, T^6, T^5, T^4, T^3, T^2, T, 1];% the end of previous segment
        Aeq_con_p(con_p_index,index+8:index+8+7) = [0, 0, 0, 0, 0, 0, 0, -1];% the begin of next segment
    end
    % beq_con_p is a zero vector

    %#####################################################
    % velocity continuity constrain between each 2 segments
    Aeq_con_v = zeros(n_seg-1, n_all_poly);
    beq_con_v = zeros(n_seg-1, 1);
    % STEP 2.5: write expression of Aeq_con_v and beq_con_v
    for con_v_index = 1:n_seg-1
        index = 1 + 8 * (con_v_index - 1);
        T = ts(con_v_index);
        Aeq_con_v(con_v_index,index:index+7) = [7*T^6, 6*T^5, 5*T^4, 4*T^3, 3*T^2, 2*T, 1, 0];% the end of previous segment
        Aeq_con_v(con_v_index,index+8:index+8+7) = [0, 0, 0, 0, 0, 0, -1, 0];% the begin of next segment
    end
    % beq_con_v is a zero vector

    %#####################################################
    % acceleration continuity constrain between each 2 segments
    Aeq_con_a = zeros(n_seg-1, n_all_poly);
    beq_con_a = zeros(n_seg-1, 1);
    % STEP 2.6: write expression of Aeq_con_a and beq_con_a
    for con_a_index = 1:n_seg-1
        index = 1 + 8 * (con_a_index - 1);
        T = ts(con_a_index);
        Aeq_con_a(con_a_index,index:index+7) = [42*T^5, 30*T^4, 20*T^3, 12*T^2, 6*T, 2, 0, 0];% the end of previous segment
        Aeq_con_a(con_a_index,index+8:index+8+7) = [0, 0, 0, 0, 0, -2, 0, 0];% the begin of next segment
    end
    % beq_con_a is a zero vector

    %#####################################################
    % jerk continuity constrain between each 2 segments
    Aeq_con_j = zeros(n_seg-1, n_all_poly);
    beq_con_j = zeros(n_seg-1, 1);
    % STEP 2.7: write expression of Aeq_con_j and beq_con_j
    for con_j_index = 1:n_seg-1
        index = 1 + 8 * (con_j_index - 1);
        T = ts(con_j_index);
        Aeq_con_j(con_j_index,index:index+7) = [210*T^4, 120*T^3, 60*T^2, 24*T, 6, 0, 0, 0];% the end of previous segment
        Aeq_con_j(con_j_index,index+8:index+8+7) = [0, 0, 0, 0, -6, 0, 0, 0];% the begin of next segment
    end
    % beq_con_j is a zero vector

    %#####################################################
    % combine all components to form Aeq and beq
    Aeq_con = [Aeq_con_p; Aeq_con_v; Aeq_con_a; Aeq_con_j];
    beq_con = [beq_con_p; beq_con_v; beq_con_a; beq_con_j];
    Aeq = [Aeq_start; Aeq_end; Aeq_wp; Aeq_con];
    beq = [beq_start; beq_end; beq_wp; beq_con];
end

author: Mesywang

repo: Motion-Planning-Algorithms


text:
function hh = nmdsfig_fill(varargin)
% Purpose: to take information about objects in multidimensional space
% and draw colored contours around them.
% Used within nmdsfig.m
%
% Usage:
%
% 1) If c is a structure from cluster_nmdsfig, which is compatible with
% nmdsfig_tools, then:
% ::
%
%    hh = nmdsfig_fill(c)
%
% Fields used are classes, groupSpace, and colors (see code for
% details)
%
% 2) Pass in arguments directly:
% ::
%
%    hh = nmdsfig_fill(classes, positions, colors)
%
% classes is n x 1 vector of group assignments (or all ones for one
% group)
%
% positions is n x 2 matrix of x, y coordinates
% colors is a cell containing as many colors as classes, {'r' 'g' 'b' ...}
% or {[1 0 0] [0 1 0] [0 0 1] ...}
%
% :Examples:
% ::
%
%    load nmdsfig_output
%    hh = nmdsfig_fill(c)
%    set(findobj('Type','Line'), 'Color', 'k')
%
% ..
%    by tor wager, feb 07
%
%    Other notes:
%    fills area around coordinates of regions in a class (group) in a color
%    using spline interpolation and other stuff.
%    designed for cluster imaging in nmdsfig figures.
%    Uses fill_area_around_points, within which
%    Many choices can be set that control how fills are done.
%    tor recommends .2 for borderscale...
% ..

hh = [];
if nargin == 0, help nmdsfig_fill, return, end

if isstruct(varargin{1})
    c = varargin{1};
    clas = c.ClusterSolution.classes;
    positions = c.GroupSpace;

    if isfield(c, 'colors')
        colors = c.colors;
    else
        colors = {'ro' 'go' 'bo' 'yo' 'co' 'mo' 'ko' 'r^' 'g^' 'b^' 'y^' 'c^' 'm^' 'k^'};
    end
else
    clas = varargin{1};
    positions = varargin{2};
    colors = varargin{3};
end

   for i = 1:max(clas)

       coords = positions(clas == i, 1:2);

       if ischar(colors{i})
        h = fill_area_around_points(coords(:,1), coords(:,2), .2, colors{i}(1));
       else
           h = fill_area_around_points(coords(:,1), coords(:,2), .2, colors{i});
       end

       if ~isempty(h) && all(ishandle(h))
           delete(h(2)) % line

            hh(i) = h(1);
       end

   end

   drawnow

end


text:
% ----------------- Add the path of DeepMIMO function --------------------%
addpath('DeepMIMO_functions')

% -------------------- DeepMIMO Dataset Generation -----------------------%
% Load Dataset Parameters
dataset_params = read_params('parameters.m');
[DeepMIMO_dataset, dataset_params] = DeepMIMO_generator(dataset_params);

% -------------------------- Output Examples -----------------------------%
% DeepMIMO_dataset{i}.user{j}.channel % Channel between BS i - User j
% %  (# of User antennas) x (# of BS antennas) x (# of OFDM subcarriers)
%
% DeepMIMO_dataset{i}.user{j}.params % Parameters of the channel (paths)
% DeepMIMO_dataset{i}.user{j}.LoS_status % Indicator of LoS path existence
% %     | 1: LoS exists | 0: NLoS only | -1: No paths (Blockage)|
%
% DeepMIMO_dataset{i}.user{j}.loc % Location of User j
% DeepMIMO_dataset{i}.loc % Location of BS i
%
% % BS-BS channels are generated only if (params.enable_BSchannels == 1):
% DeepMIMO_dataset{i}.basestation{j}.channel % Channel between BS i - BS j
% DeepMIMO_dataset{i}.basestation{j}.loc
% DeepMIMO_dataset{i}.basestation{j}.LoS_status
%
% % Recall that the size of the channel vector was given by
% % (# of User antennas) x (# of BS antennas) x (# of OFDM subcarriers)
% % Each of the first two channel matrix dimensions follows a certain
% % reshaping sequence that can be obtained by the following
% % 'antennamap' vector: Each entry is 3 integers in the form of
% % 'xyz' where each representing the antenna number in x, y, z directions
% antennamap = antenna_channel_map(params.num_ant_BS(1), ...
%                                  params.num_ant_BS(2), ...
%                                  params.num_ant_BS(3), 1);
%
% -------------------------- Dynamic Scenario ----------------------------%
%
% DeepMIMO_dataset{f}{i}.user{j}.channel % Scene f - BS i - User j
% % Every other command applies as before with the addition of scene ID
% params{f} % Parameters of Scene f
%
% ------------------------------------------------------------------------%

author: DeepMIMO

repo: DeepMIMO-matlab