text:
function selected_points = test_terrain()
% NOTEST

import iris.terrain_grid.*;
import iris.inflate_region;
import iris.thirdParty.polytopes.*;
import iris.cspace.*;
import iris.drawing.*;
lcmgl = drake.util.BotLCMGLClient(lcm.lcm.LCM.getSingleton(), 'segmentation');
% load('example_feas_map.mat','Q');
load('data/example_heights.mat','heights','px2world');
load('data/example_state.mat', 'x0');

px2world(1,end) = px2world(1,end) - sum(px2world(1,1:3)); % stupid matlab 1-indexing...
px2world(2,end) = px2world(2,end) - sum(px2world(2,1:3));
resample = 2;
mag = 2^(resample-1);
heights = interp2(heights, (resample-1));
px2world = px2world * [1/mag 0 0 (1-1/mag); 0 1/mag 0 (1-1/mag ); 0 0 1 0; 0 0 0 1];

world2px = inv(px2world);
world2px_2x3 = world2px(1:2,[1,2,4]);
px2world_2x3 = px2world(1:2, [1,2,4]);

Q = imfilter(heights, [1, -1]) - (0.06/mag) > 0;
Q = Q | imfilter(heights, [-1, 1]) - (0.06/mag) > 0;
Q = Q | imfilter(heights, [1; -1]) - (0.06/mag) > 0;
Q = Q | imfilter(heights, [-1; 1]) - (0.06/mag) > 0;
Q(isnan(heights)) = 1;

% grid = Q(85:125,25:85);
grid = ~Q;
lb = [180;70;-pi];
ub = [245; 160; pi];
A_bounds = [-1,0,0;
            0,-1,0;
            0,0,-1;
            1,0,0;
            0,1,0;
            0,0,1];
b_bounds = [-lb;ub];

bot = world2px_2x3 * 0.9*bsxfun(@minus, [-0.0820,-0.0820, 0.1780, 0.1780;
                                         0.0624,-0.0624, 0.0624,-0.0624;
                                         0,0,0,0], [0.0480; 0; 0;]);
% obstacles = segment_grid(~grid(lb(1):ub(1),lb(2):ub(2)));
% for j = 1:length(obstacles)
%   obstacles{j} = bsxfun(@plus, obstacles{j}, lb(1:2)-1);
% end
% obstacles = cspace3(obstacles, bot, 4);
selected_points = [];

while true
  figure(6)
  imshow(grid, 'InitialMagnification', 'fit')
  [c,r] = ginput(1);
  if ~isempty(c)
    selected_points(:,end+1) = [r;c];
  end
  if isempty(c)
    break
  end
  c = round(c);
  r = round(r);
  black_edges = [];
  [black_edges(1,:), black_edges(2,:)] = ind2sub(size(grid), find(component_boundary(grid, [r;c])));
  obs_mask = all(bsxfun(@minus, A_bounds([1,2,4,5],1:2) * black_edges, b_bounds([1,2,4,5])) <= max(max(bot)));
  obstacles = mat2cell(black_edges(:,obs_mask) , 2, ones(1,sum(obs_mask)));
  obstacles = cspace3(obstacles, bot, 4);

  [A,b,C,d,results] = inflate_region(obstacles, A_bounds, b_bounds, [r;c;0]);
%   animate_results(results);

  [inner_poly, outer_poly] = project_c_space_region(A,b);

  % px_xy = [V(:,2)'; V(:,1)'];
  if ~isempty(inner_poly)
    px_xy = [inner_poly(2,:); inner_poly(1,:)];
    world_xy = px2world_2x3 * [px_xy; ones(1,size(px_xy,2))];
    z = ones(size(px_xy,2)) * heights(r,c) + 0.03;
    world_xyz = [world_xy; z];
    lcmgl.glColor3f(0,1,0);
    lcmgl.glLineWidth(10);
    lcmgl.glBegin(lcmgl.LCMGL_LINES);
    for j = 1:size(world_xy,2)-1
      lcmgl.glVertex3d(world_xyz(1,j),world_xyz(2,j),world_xyz(3,j));
      lcmgl.glVertex3d(world_xyz(1,j+1),world_xyz(2,j+1),world_xyz(3,j+1));
    end
    lcmgl.glVertex3d(world_xyz(1,end),world_xyz(2,end),world_xyz(3,end));
    lcmgl.glVertex3d(world_xyz(1,1),world_xyz(2,1),world_xyz(3,1));
    lcmgl.glEnd();
  end

  px_xy = [outer_poly(2,:); outer_poly(1,:)];
  world_xy = px2world_2x3 * [px_xy; ones(1,size(px_xy,2))];
  z = ones(size(px_xy,2)) * heights(r,c) + 0.03;
  world_xyz = [world_xy; z];
  lcmgl.glColor3f(1,1,0);
  lcmgl.glLineWidth(10);
  lcmgl.glBegin(lcmgl.LCMGL_LINES);
  for j = 1:size(world_xy,2)-1
    lcmgl.glVertex3d(world_xyz(1,j),world_xyz(2,j),world_xyz(3,j));
    lcmgl.glVertex3d(world_xyz(1,j+1),world_xyz(2,j+1),world_xyz(3,j+1));
  end
  lcmgl.glVertex3d(world_xyz(1,end),world_xyz(2,end),world_xyz(3,end));
  lcmgl.glVertex3d(world_xyz(1,1),world_xyz(2,1),world_xyz(3,1));
  lcmgl.glEnd();

  lcmgl.glColor3f(0,1,0);
  lcmgl.sphere([(px2world_2x3 * [c;r;1])', heights(r,c)], 0.05, 20, 20);

%   lcmgl.switchBuffers();
%   plot_lcm_points([px_xy', px_z'],repmat([1,0,0],length(px_z),1),101,'region1',3,1);
end
lcmgl.switchBuffers();
end

author: rdeits

repo: iris-distro

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


text:
function [u,s,U_r,U_s,U_pk,U_pd,U_l]  = projAhmPntIntoOmniCamOnRob(Rf, Sf, Spk, Spd, l)

% PROJAHMPNTINTOOMNICAMONROB Project Ahm pnt into omnidirectional camera on robot.
%    [U,S] = PROJAHMPNTINTOOMNICAMONROB(RF, SF, SPK, SPD, L) projects 3D
%    anchored homogeneous points into a omni-cam mounted on a robot,
%    providing also the non-measurable depth. The input parameters are:
%       RF : robot frame
%       SF : omni-cam sensor frame in robot
%       SPK: omni-cam intrinsic parameters [xc yc c d e]'
%       SPD: radial distortion polynom [a0 a1 a2 ...]'
%       L  : 3D anchored homog. point [x y z vx vy vz rho]'
%    The output parameters are:
%       U  : 2D pixel [u v]'
%       S  : non-measurable depth
%
%    The function accepts an ahm points matrix L = [L1 ... Ln] as input.
%    In this case, it returns a pixels matrix  U = [U1 ... Un] and a depths
%    row-vector S = [S1 ... Sn].
%
%    [U,S,U_R,U_S,U_K,U_D,U_L] = ... gives also the jacobians of the
%    observation U wrt all input parameters. Note that this only works for
%    single points.
%
%    See also PINHOLE, TOFRAME, PROJEUCPNTINTOPINHOLEONROB.
%

%   Copyright 2012 Grigory Abuladze @ ASL-vision
%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

if nargout <= 2  % No Jacobians requested

    p     = ahm2euc(l);
    [u,s] = projEucPntIntoOmniCamOnRob(Rf, Sf, Spk, Spd, p);

else            % Jacobians requested

    if size(l,2) == 1

        % function calls
        [p,P_l]                     = ahm2euc(l);
        [u,s,U_r,U_s,U_pk,U_pd,U_p] = projEucPntIntoOmniCamOnRob(Rf, Sf, Spk, Spd, p);

        % chain rule
        U_l = U_p*P_l;

    else
        error('??? Jacobians not available for multiple AHM points.')

    end

end

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

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

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function f = testFile(d1,d2,d3,d4,d5,d6,dth1,dth2,dth3,dth4,dth5,dth6,g,l1,l2,l3,l4,l5,m1,m2,m3,m4,m5,m6,th1,th2,th3,th4,th5,th6)
%TESTFILE
%    F = TESTFILE(D1,D2,D3,D4,D5,D6,DTH1,DTH2,DTH3,DTH4,DTH5,DTH6,G,L1,L2,L3,L4,L5,M1,M2,M3,M4,M5,M6,TH1,TH2,TH3,TH4,TH5,TH6)

%    This function was generated by the Symbolic Math Toolbox version 6.0.
%    04-Nov-2014 12:56:21

t2 = cos(th1);
t3 = dth2.^2;
t4 = th1-th2;
t5 = sin(t4);
t6 = dth3.^2;
t7 = th1-th3;
t8 = sin(t7);
t9 = dth4.^2;
t10 = th1-th4;
t11 = sin(t10);
t12 = dth5.^2;
t13 = th1-th5;
t14 = sin(t13);
t15 = cos(th2);
t16 = dth6.^2;
t17 = dth1.^2;
t18 = th2-th3;
t19 = sin(t18);
t20 = th2-th4;
t21 = sin(t20);
t22 = th2-th5;
t23 = sin(t22);
t24 = cos(th3);
t25 = th3-th4;
t26 = sin(t25);
t27 = th3-th5;
t28 = sin(t27);
t29 = cos(th4);
t30 = th4-th5;
t31 = sin(t30);
t32 = cos(th5);
t33 = th1-th6;
t34 = sin(t33);
t35 = th2-th6;
t36 = sin(t35);
t37 = th3-th6;
t38 = sin(t37);
t39 = th4-th6;
t40 = sin(t39);
t41 = th5-th6;
t42 = sin(t41);
f = [-d1.*g.*m1.*t2-g.*l1.*m2.*t2-g.*l1.*m3.*t2-g.*l1.*m4.*t2-g.*l1.*m5.*t2-g.*l1.*m6.*t2-d2.*l1.*m2.*t3.*t5-d3.*l1.*m3.*t6.*t8-d4.*l1.*m4.*t9.*t11-d5.*l1.*m5.*t12.*t14-d6.*l1.*m6.*t16.*t34-l1.*l2.*m3.*t3.*t5-l1.*l2.*m4.*t3.*t5-l1.*l2.*m5.*t3.*t5-l1.*l2.*m6.*t3.*t5-l1.*l3.*m4.*t6.*t8-l1.*l3.*m5.*t6.*t8-l1.*l3.*m6.*t6.*t8-l1.*l4.*m5.*t9.*t11-l1.*l4.*m6.*t9.*t11-l1.*l5.*m6.*t12.*t14;-d2.*g.*m2.*t15-g.*l2.*m3.*t15-g.*l2.*m4.*t15-g.*l2.*m5.*t15-g.*l2.*m6.*t15+d2.*l1.*m2.*t5.*t17-d3.*l2.*m3.*t6.*t19-d4.*l2.*m4.*t9.*t21-d5.*l2.*m5.*t12.*t23-d6.*l2.*m6.*t16.*t36+l1.*l2.*m3.*t5.*t17+l1.*l2.*m4.*t5.*t17+l1.*l2.*m5.*t5.*t17+l1.*l2.*m6.*t5.*t17-l2.*l3.*m4.*t6.*t19-l2.*l3.*m5.*t6.*t19-l2.*l3.*m6.*t6.*t19-l2.*l4.*m5.*t9.*t21-l2.*l4.*m6.*t9.*t21-l2.*l5.*m6.*t12.*t23;-d3.*g.*m3.*t24-g.*l3.*m4.*t24-g.*l3.*m5.*t24-g.*l3.*m6.*t24+d3.*l2.*m3.*t3.*t19+d3.*l1.*m3.*t8.*t17-d4.*l3.*m4.*t9.*t26-d5.*l3.*m5.*t12.*t28-d6.*l3.*m6.*t16.*t38+l2.*l3.*m4.*t3.*t19+l2.*l3.*m5.*t3.*t19+l1.*l3.*m4.*t8.*t17+l2.*l3.*m6.*t3.*t19+l1.*l3.*m5.*t8.*t17+l1.*l3.*m6.*t8.*t17-l3.*l4.*m5.*t9.*t26-l3.*l4.*m6.*t9.*t26-l3.*l5.*m6.*t12.*t28;-d4.*g.*m4.*t29-g.*l4.*m5.*t29-g.*l4.*m6.*t29+d4.*l2.*m4.*t3.*t21+d4.*l1.*m4.*t11.*t17+d4.*l3.*m4.*t6.*t26-d5.*l4.*m5.*t12.*t31-d6.*l4.*m6.*t16.*t40+l2.*l4.*m5.*t3.*t21+l2.*l4.*m6.*t3.*t21+l1.*l4.*m5.*t11.*t17+l1.*l4.*m6.*t11.*t17+l3.*l4.*m5.*t6.*t26+l3.*l4.*m6.*t6.*t26-l4.*l5.*m6.*t12.*t31;-d5.*g.*m5.*t32-g.*l5.*m6.*t32+d5.*l2.*m5.*t3.*t23+d5.*l1.*m5.*t14.*t17+d5.*l3.*m5.*t6.*t28+d5.*l4.*m5.*t9.*t31-d6.*l5.*m6.*t16.*t42+l2.*l5.*m6.*t3.*t23+l1.*l5.*m6.*t14.*t17+l3.*l5.*m6.*t6.*t28+l4.*l5.*m6.*t9.*t31;d6.*m6.*(-g.*cos(th6)+l2.*t3.*t36+l3.*t6.*t38+l1.*t17.*t34+l4.*t9.*t40+l5.*t12.*t42)];

author: MatthewPeterKelly

repo: dscTutorials


text:
function [f]= spm_fx_dem_salience(x,v,P)
% returns the flow for visual search
% FORMAT [f]= spm_fx_dem_salience(x,v,P)
%
% x    - hidden states:
%   o(1) - oculomotor angle
%   o(2) - oculomotor angle
%   x(1) - relative amplitude of visual hypothesis 1
%   x(2) - relative amplitude of visual hypothesis 2
%   x(3) - ...
%
% v    - hidden causes - attracting location
% P    - parameters
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_fx_dem_salience.m 4862 2012-08-24 19:21:27Z karl $

% intialise flow (to ensure fields are aligned)
%--------------------------------------------------------------------------
f    = x;

% motion of oculomotor angles (attracted to target)
%==========================================================================
f.o  = (v - x.o)/4;

% motion of hypothesis states (with competition)
%==========================================================================
f.x = 1 - sum(exp(x.x)) - x.x/128;

author: spm


text:
function TDD = divdiff(X, Y)
%
% TDD = divdiff(X, Y)
%
% DIVDIFF
%   Newton's Method for Divided Differences.
%
%   The following formula is solved:
%       Pn(x) = f(x0) + f[x0,x1](x-x0) + f[x0,x1,x2](x-x0)(x-x1) + ...
%           + f[x0,x1,..,xn](x-x0)(x-x1)..(x-x[n-1])
%       where f[x0,x1] = (f(x1-f(x0))/(x1-x0)
%             f[x0,x1,..,xn] = (f[x1,..,xn]-f[x0,..,x_[n-1]])/(xn-x0)
%
% NOTE: f^(n+1)(csi)/(n+1)! aprox. = f[x0,x1,..,xn,x_[n+1]]
%
% Input::
%	X = [ x0 x1 .. xn ] - object vector
%	Y = [ y0 y1 .. yn ] - image vector
%
% Output:
%   TDD - table of divided differences
%
% Example:
%   TDD = divdiff( [ 1.35 1.37 1.40 1.45 ], [ .1303 .1367 .1461 .1614 ])
%
% Author:	Tashi Ravach
% Version:	1.0
% Date:     22/05/2006
%

    if nargin ~= 2
        error('divdiff: invalid input parameters');
    end

    [ p, m ] = size(X); % m points, polynomial order <= m-1
    if p ~= 1 || p ~=size(Y, 1) || m ~= size(Y, 2)
        error('divdiff: input vectors must have the same dimension');
    end

    TDD = zeros(m, m);
    TDD(:, 1) = Y';
    for j = 2 : m
        for i = 1 : (m - j + 1)
            TDD(i,j) = (TDD(i + 1, j - 1) - TDD(i, j - 1)) / (X(i + j - 1) - X(i));
        end
    end

end

author: Sable

repo: mcbench-benchmarks


text:
function n = p29_n ( )

%*****************************************************************************80
%
%% P29_N returns the number of variables for problem 29.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 January 2001
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer N.  If N is positive, N represents the
%    only legal value for N for this problem.  If N is
%    negative, the absolute value of N is the least legal
%    value of N, but other values are allowable.
%
  n = 2;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function test_one ( )

%*****************************************************************************80
%
%% TEST_ONE performs interpolation with SPINTERP.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 January 2013
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST_ONE:\n' );
  fprintf ( 1, '  An example of the use of SPINTERP\n' );
  fprintf ( 1, '  for sparse interpolation in multiple dimensions.\n' );

  addpath ( '../spinterp' );

  options = spset ( ...
    'GridType', 'Chebyshev', ...
    'MinDepth', 11 );

  m = 2;

  box = [ 0.0, 1.0; ...
          0.0, 1.0 ];
%
%  Set up the sparse grid structure for the function.
%
  c = spvals ( @f_one, m, box, options );
%
%  C is a structure, so this is the easiest way to
%  print it out and examine it:
%
  c
%
%  Pick some sample points in the unit square
%  and compare the interpolant Z to the function F_ONE(X,Y).
%
  sample_num = 100;

  x = rand ( 1, sample_num );
  y = rand ( 1, sample_num );
  z = spinterp ( c, x, y );

  error = max ( abs ( z - f_one ( x, y ) ) );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Maximum interpolation error was %g\n', error );
%
%  Plot the function and the interpolant.
%
  subplot ( 1, 2, 1 );
  ezmesh ( @f_one, [ 0.0, 1.0 ] );
  title ( 'f(x,y) = sin(x) + cos(y)' );
  subplot ( 1, 2, 2 );
  ezmesh ( @(x,y) spinterp ( c, x, y ), [ 0.0, 1.0 ] );
  title ( 'Sparse grid interpolant' );

  rmpath ( '../spinterp' )
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST_ONE:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function [ vX ] = LevinsonRecursion( mT, vY )
% ----------------------------------------------------------------------------------------------------- %
%[ vX ] = LevinsonRecursion( mT, vY )
% Solves Linear System mT * vX = vY where the Matrix 'mT' is a Toeplitz
% Matrix using Levinson Recursion improved speed (Complexity of N^2 instead
% of N^3). Input:
%   - mT    -   Input Matrix.
%               Matrix with a Toeplitz Structure.
%               Structure: Matrix (N x N).
%               Type: 'Single' / 'Double'.
%               Range: (-inf, inf).
%   - vY    -   Input Vector.
%               Structure: Column Vector (N x 1).
%               Type: 'Single' / 'Double'.
%               Range: (-inf, inf).
% Output:
%   - vX    -   Output Vector.
%               The solution of the system mT * vX = vY.
%               Structure: Column Vector (N x 1).
%               Type: 'Single' / 'Double'.
%               Range: (-inf, inf).
% References
%   1.  Levinson Recursion (Wikipedia) - https://en.wikipedia.org/wiki/Levinson_recursion
% Remarks:
%   1.  Joint work with Yair Shemer.
%   2.  This implementation supports only rectangular matrices.
% TODO:
%   1.  A
% Release Notes:
%   -   1.0.000     15/02/2017  Royi Avital
%       *   First realease version.
% ----------------------------------------------------------------------------------------------------- %

numRows = size(mT, 1); %<! Square Matrix

vF = zeros([numRows, 1]); %<! Forward Vector
vB = zeros([numRows, 1]); %<! Backward Vector
vX = zeros([numRows, 1]);

% Initialized vF, vB and vX;
vF(1) = 1 / mT(1, 1);
vB(1) = 1 / mT(1, 1);
vX(1) = vY(1) / mT(1, 1);

for iRow = 2 : numRows
    % Calculate the epsilons (error expressions)
    epsF = mT(iRow, 1:(iRow - 1)) * vF(1:(iRow - 1));
    epsB = mT(1, 2:iRow) * vB(1:(iRow - 1));
    epsX = mT(iRow, 1:(iRow - 1)) * vX(1:(iRow - 1));

    % Given vF, vB and vX with size n-1 find those vectors with size n
    vFPrev      = vF;
    scalingFctr = 1 / (1 - (epsF * epsB));

    vF(1:iRow) = scalingFctr * vF(1:iRow) - (epsF * scalingFctr) * [0; vB(1:(iRow - 1))];
    vB(1:iRow) = scalingFctr * [0; vB(1:(iRow - 1))] - (epsB * scalingFctr) * vFPrev(1:iRow);
    vX(1:iRow) = vX(1:iRow) + (vY(iRow) - epsX) * vB(1:iRow);
end

end

author: RoyiAvital

repo: StackExchangeCodes


% Applies a frequency shift to time domain samples
%
% @param samples Input samples (must be row or column vector)
% @param sample_rate Sample rate (in Hz) that the frequency offset should be relative to
% @param freq_offset Frequency offset (in Hz) to apply
% @return samples Frequency shifted time domain samples
function [samples] = apply_freq_offset_time_domain(samples, sample_rate, freq_offset)
    assert(isrow(samples) || iscolumn(samples), "Samples must be a row or column vector");

    if (isrow(samples))
        samples = samples .* exp(1j * 2 * pi * freq_offset / sample_rate * [0:length(samples)-1]);
    else
        samples = samples .* exp(1j * 2 * pi * freq_offset / sample_rate * [0:length(samples)-1].');
    end
end


function mse = get_mse(V_equal,W_equal,Weighted)
%compute mse or wmse

global Vn H Ns;

H_equal = W_equal' * H * V_equal;
E_matrix = H_equal * H_equal' - H_equal - H_equal' + eye(Ns) + Vn * W_equal'*W_equal;

%normal MSE
if nargin == 2
    mse = trace(E_matrix);
end

if nargin == 3
    mse = trace(Weighted * E_matrix);
end


text:
function s=unsplit_string(c, sep)
% 'zips' a cell with strings and a seperator symbol; the inverse of
% SPLIT_STRING
%
% S=UNSPLIT_STRING(C,SEP) takes a cell with strings C and a seperator
% string SEP, and returns the string [C{1} SEP C{2} SEP ... SEP C{end}].
%
% It is assumed that C and SEP are strings with exactly one row.
%
% NNO Jan 2010

if ~iscellstr(c)
    error('Expected cell string input');
end

n=numel(c);
if n==0
    s='';
    return
end

% insert separator string
parts=cell(1,2*n-1);
for k=1:n
    parts{k*2-1}=c{k};
    if k<n
        parts{k*2}=sep;
    end
end

% concatenate
s=sprintf('%s',parts{:});

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% Convert bin-epoched EEG to a bin-organized EEG structure
% may be useful for subsequent decoding scripts
% axs 2020
%
% Inputs:
%    EEG - an EEG set, which must have been processed with ERPLAB, with
%         valid EventList, and be bin-epoched EEG data
%    write_file_path - optional full path to write file to, string
%
% Outputs:
%  binorgEEG - reorganisation of the EEG data, around the ERPLAB bins
%              with the data in binorgEEG.binwise_data(bin_n).data(elec,timepoints,trial_n)
function binorgEEG = binepEEG_to_binorgEEG(EEG, write_file_path)

dim_data = numel(size(EEG.data));
if dim_data ~= 3
    error('binepEEG_to_binorgEEG works on bin-epoched EEG. Please ensure data is epoched (not continuous) and has bin labels');
end

% Prepare info for averager call
artif = 1; stderror = 1; excbound = 1;

% Call ERP Averager subfunction, populating the epoch_list

[ERP2, EVENTLISTi, countbiORI, countbinINV, countbinOK, countflags, workfname,epoch_list] = averager(EEG, artif, stderror, excbound);

% We now have epoch_list, which has a list of 'good' bepoch indecies to use
% epoch_list(1).good_bep_indx(:) has those trials for bin 1
% epoch_list(2).good_bep_indx(:) has those trials for bin 2, etc
nbin = ERP2.nbin;

% Prepare binorg data structure
for bin = 1:nbin

    binorgEEG.binwise_data(bin).data = EEG.data(:,:,epoch_list(bin).good_bep_indx);
    binorgEEG.n_trials_this_bin(bin) = numel(epoch_list(bin).good_bep_indx);

end

% Now, the data from each trial is saved, reorganized by bin, such that:
% for bin X
%  binorgEEG.binwise_data(X).data(:,:,1) gives all elecs and all
%  timepoints for the 1st trial of bin X

if ischar(write_file_path)
    % If a string of a full write file path is provided in arg, save it there
    save(write_file_path,'binorgEEG')

elseif write_file_path == 1
    % UI file pick when write_file_path == 1
    [file_name, file_path] = uiputfile('*.mat','Please pick a path to save the Bin-Organized EEG data');
    write_file_path = [file_path file_name];
    save(write_file_path,'binorgEEG');
end

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
function error=matlab2tsai_error(X,camera,uv_d,xy_matlab_N)

% JMM Montiel. Sept 2005
%
% Error between the Tsia undistortion model and the
%   Swaminathan model
% Input
%   k1      -Swaminathan distortion model
%   camera  -Camera parameters. camera.k1 is irrelevant
%   uv_d    -Distorted point coordinates
%   uv_u    -Undistorted coordinates according to Tsai model

k1=X(1);
k2=X(2);
f=X(3);
cx = X(4);
cy = X(5);

camera.f = f;
camera.Cx = cx;
camera.Cy = cy;
camera.k1 = k1;
camera.k2 = k2;

xy_tsai_N = undistort_and_normalize_fm( uv_d, camera );
if(1)
    cla;
    plot(xy_matlab_N(1,:),xy_matlab_N(2,:),'+r');
    hold on
    plot(xy_tsai_N(1,:),xy_tsai_N(2,:),'+g');
    xlabel('red, matlab.    green Tsai with 1 radial distorion parameter')
end
error=[xy_tsai_N(1,:)-xy_matlab_N(1,:);xy_tsai_N(2,:)-xy_matlab_N(2,:)];

repo: ekfmonoslam


text:
%%% Extraction of the final intrinsic and extrinsic paramaters:

check_active_images;

if ~exist('solution_error')
   solution_error = zeros(6*n_ima + 15,1);
end;

fc = solution(1:2);%***
cc = solution(3:4);%***
alpha_c = solution(5);%***
kc = solution(6:10);%***

fc_error = solution_error(1:2);
cc_error = solution_error(3:4);
alpha_c_error = solution_error(5);
kc_error = solution_error(6:10);

% Calibration matrix:

KK = [fc(1) fc(1)*alpha_c cc(1);0 fc(2) cc(2); 0 0 1];
inv_KK = inv(KK);

% Extract the extrinsic paramters, and recomputer the collineations

for kk = 1:n_ima,

   if active_images(kk),

      omckk = solution(15+6*(kk-1) + 1:15+6*(kk-1) + 3);%***
      Tckk = solution(15+6*(kk-1) + 4:15+6*(kk-1) + 6);%***

      omckk_error = solution_error(15+6*(kk-1) + 1:15+6*(kk-1) + 3);
      Tckk_error = solution_error(15+6*(kk-1) + 4:15+6*(kk-1) + 6);

   	Rckk = rodrigues(omckk);

   	Hkk = KK * [Rckk(:,1) Rckk(:,2) Tckk];

   	Hkk = Hkk / Hkk(3,3);

   else

      omckk = NaN*ones(3,1);
      Tckk = NaN*ones(3,1);
      Rckk = NaN*ones(3,3);
      Hkk = NaN*ones(3,3);
      omckk_error = NaN*ones(3,1);
      Tckk_error = NaN*ones(3,1);

   end;

   eval(['omc_' num2str(kk) ' = omckk;']);
   eval(['Rc_' num2str(kk) ' = Rckk;']);
   eval(['Tc_' num2str(kk) ' = Tckk;']);
   eval(['H_' num2str(kk) '= Hkk;']);
   eval(['omc_error_' num2str(kk) ' = omckk_error;']);
   eval(['Tc_error_' num2str(kk) ' = Tckk_error;']);

end;


text:
function result = lasso_cv(y, x, k, options)
% -------------------------------------------------------------------------
% function result = lasso_cv(y, x, k, options)
% -------------------------------------------------------------------------
% PURPOSE: fits the parameters of a linear model by using the lasso and
%          k-folds cross validation
% -------------------------------------------------------------------------
% INPUTS:
% y:         the values of the dependent variable
% x:         the values of the explanatory variables - constants will be
%            ignored
% k:         determine number of cross validation folds
% n_div:
% alignment:
% cv_assignment: allows that the division within groups to be defined by
%                the user. cv_assignment is a nx1 vector where n is the
%                number of available observations. The second columns
%                defines which observations belongs to each group.
%                If this
%                parameter is used, k is ignored and set to be equal to the
%                number of unique elements of cv_assignment;
% -------------------------------------------------------------------------
% OUTPUTS:
% a structure containing:
% all the fields contained in the strucutre returned by the LASSO function
% (see the lasso function help)  plus the following fields:
% result.CV.intercept    : intercept for the model chosen by CV
% result.CV.betas        : coefficients for the model selected by CV
% result.CV.fit          : Fitted values for the cross-validated s
% result.CV.residuals    : Residuals for the cross-validated s
% result.CV.SSR          : Sum of squared residuals for the cross-validated
%                          estimate
% result.CV.MSE          : minimum cross-validated MSE
% result.CV.cv_assignment: a matrix determining which observations where
%                          assigned to which CV groups
% result.CV.k            : number of CV folds
% -------------------------------------------------------------------------
% WARNING:
% unless cv_assignment is provided, this function invokes ASSIGN_CV.
% ASSIGN_CV calls randperm affecting the state of the random number gen.
% -------------------------------------------------------------------------
% Author: Guilherme V. Rocha
%         Department of Statistics
%         University of California, Berkeley
%         gvrocha@stat.berkeley.edu, gvrocha@gmail.com
% 2006/09
% -------------------------------------------------------------------------
% See also: LASSO, ASSIGN_CV, RANDPERM

if(nargin < 4)
  options = [];
end;

if(~isfield(options, 'alignment'))
  options.alignment = 'normalized_penalty';
end;
if(~isfield(options, 'cv_assignment'))
  assign_cv_flag = 1;
else
  assign_cv_flag = 0;
  cv_assignment  = options.cv_assignment;
  k              = length(unique(cv_assignment));
end;
if(~isfield(options, 'num_div'))
  options.num_div =100;
end;
n_div = options.num_div;
alignment = options.alignment;

n_obs           = size(y, 1);
max_index       = -Inf;
if(assign_cv_flag)
  cv_assignment = assign_cv(n_obs, k);
end;

for i = 1:k
  fit_indexes   = find(cv_assignment~=i);
  y_fit         = y(fit_indexes);
  x_fit         = x(fit_indexes,:);
  res(i)        = lasso_rocha(y_fit, x_fit, options);
  switch lower(alignment)
    case 'normalized_penalty'
      max_index = max(max_index, max(res(i).npenalty));
    case 'penalty'
      max_index = max(max_index, max(res(i).penalty));
    case 'lambda'
      max_index = max(max_index, max(res(i).lambda));
  end;
end;

index_slices = 0:max_index/n_div:max_index;

for i = 1:k
  interpol    = lasso_coefficients(res(i), index_slices, alignment);
  cv_indexes  = find(cv_assignment==i);
  y_cv        = y(cv_indexes);
  x_cv        = x(cv_indexes,:);

  y_pred      = lasso_predict(res(i), x_cv, index_slices, options);
  residuals   = repmat(y_cv, 1, size(y_pred, 2)) - y_pred;
  if(size(residuals, 1)>1)
    MSE(i,:)    = mean(residuals.^2);
  else
    MSE(i,:)    = residuals.^2;
  end;
end;

[min_MSE, best_cv_MSE] = min(mean(MSE));
best_index             = index_slices(best_cv_MSE);
result                 = lasso_rocha(y, x, options);

result.CV                = lasso_coefficients(result, best_index, alignment);
result.CV.cv_index       = best_cv_MSE;
result.CV.cv_assignment  = cv_assignment;
result.CV.k              = k;
result.CV.min_MSE        = min_MSE;
result.CV.MSEs           = MSE;
result.CV.path           = lasso_coefficients(result, index_slices, alignment);
result.options           = options;
rm_fields                = {'ny', 'nx', 'xtx', 'xty', 'beta', 'intercept', 'penalty', 'xtrs'};
result.CV.fold_results   = rmfield(res, intersect(fieldnames(res), rm_fields));

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
classdef CoSaMP < handle
    % Implements CoSaMP-MMV algorithm for sparse recovery

    properties
        % These properties can be configured before running cosamp
        % Default threshold
        errorNormThreshold = 1e-6;
        % Maximum number of iterations for approximation
        MaxIters
        Verbose = false
        % Indicates if matching pursuit should be used for identification
        % UseMPIdentification = false
        % Indicates if least squares should be run on final support
        LSOnFinalSupport = false
        % The norm to be chosen for rows
        P
        % Indicates that residuals should be orthogonalized for rank awareness
        RankAwareResidual = false
    end

    properties(SetAccess=private)
        % The dictionary
        Dict
        % Ambient signal dimensions
        N
        % Number of atoms in dictionary
        D
        % Sparsity level of representations (may be negative)
        K
        % Result of a solver
        result
    end

    methods
        function self  = CoSaMP(Dict, K, P, options)
            if nargin < 3
                % By default we apply l_1 norm on rows
                P = 1;
            end
            if P ~= 1 &&  P ~= 2
                error('Only l_1 and l_2 norms are supported.');
            end
            self.P = P;
            % We assume that all the columns in dictionary are normalized.
            if isa(Dict, 'spx.dict.Operator')
                self.Dict = Dict;
            elseif ismatrix(Dict)
                self.Dict = spx.dict.MatrixOperator(Dict);
            else
                error('Unsupported operator.');
            end
            [self.N, self.D] = size(Dict);
            if nargin < 2
                % No sparsity level has been pre-specified.
                % We make an estimate of sparsity level
                % based on phase transition analysis by Donoho
                K = round((self.N / 2) * log (self.D));
            end
            self.K = K;
            % Maximum number of iterations
            maxIter = 30;
            self.MaxIters = maxIter;
            if nargin >= 4
                % options have been specified
                if isfield(options, 'RankAwareResidual')
                    self.RankAwareResidual = options.RankAwareResidual;
                end
            end
        end

        function result  = solve(self,Y)
            n = self.N;
            d = self.D;
            k = self.K;
            % The number of signals being approximated.
            s = size(Y, 2);
            dict = self.Dict;
            % Current estimate
            solution_nz_mat = zeros(k,s);
            % Current residual
            residual_mat = Y;
            % Number of iterations
            iterations = 0;
            y_norm = norm(Y, 'fro');
            old_residual_norm = 1;
            min_residual_norm = old_residual_norm;
            errorNormThreshold = self.errorNormThreshold;
            maxIterations = self.MaxIters;
            result.halted_on_max_iter = false;
            result.halted_on_residual_norm = false;
            result.halted_on_norm_change = false;
            result.halted_on_support_change = false;
            % K indices for current support
            current_support = [];
            while true
                iterations  = iterations+1;
                % We identify the support for largest 2K entries
                extra = 2*k;
                % if isempty(current_support)
                %     extra = k + extra;
                % end
                if self.RankAwareResidual
                    residual_mat = orth(residual_mat);
                end
                % Compute the proxy
                proxy_mat = dict.apply_ctranspose(residual_mat);
                % identify the largest entries
                largest_2k  = spx.pursuit.joint.CoSaMP.largest_k_l2(proxy_mat, extra);
                % We now compute our support
                % build an array holding up to 3 K largest indices
                support_3k = false(d, 1);
                % the older K columns
                support_3k(current_support) = true;
                % New 2K columns
                support_3k(largest_2k) = true; % T <= 3K
                % We pickup corresponding columns from Phi
                subdict = dict.columns(support_3k); % MxT
                % We compute signal estimate over these columns
                % B_subdict =  linsolve(subdict, Y); % TxM * MxS
                B_subdict = spx.pursuit.joint.CoSaMP.least_squares_on_support(subdict, Y);
                % sort them in descending row norm order
                pruned_largest_k =  spx.pursuit.joint.CoSaMP.largest_k_l2(B_subdict, k);
                % keep only first k rows as new solution
                solution_nz_mat = B_subdict(pruned_largest_k, :);
                old_support = current_support;
                % update current support by choosing largest k indices
                support_3k = find(support_3k == 1);
                current_support = support_3k(pruned_largest_k);
                % compute the measurement estimate
                y_estimate_matrix = dict.apply_columns(...
                    solution_nz_mat, current_support);
                % compute the new residual
                residual_mat = Y - y_estimate_matrix;
                % compute the residual norm
                residual_norm = norm(residual_mat, 'fro');
                residual_norm = residual_norm / y_norm;
                % how many indices have been added and/or removed
                support_change = length(union(old_support,current_support)) - length(intersect(old_support,current_support));
                % if residual_norm > 1.2*min_residual_norm
                %     % We are diverging
                %     break;
                % end
                if residual_norm < old_residual_norm
                    improvement = (old_residual_norm - residual_norm) / old_residual_norm;
                    % if improvement < 1e-8
                    %     % No improvement in this iteration
                    %     result.halted_on_norm_change = true;
                    %     break;
                    % end
                end
                if(residual_norm < errorNormThreshold)
                    result.halted_on_residual_norm = true;
                    break;
                end
                if iterations >= maxIterations
                    % Too many iterations we are going nowhere
                    result.halted_on_max_iter = true;
                    break;
                end
                if ~isempty(old_support) && support_change == 0
                    result.halted_on_support_change = true;
                    break;
                end
                % TODO add support for detection of no change
                % in support
                old_residual_norm = residual_norm;
                min_residual_norm = min(residual_norm, min_residual_norm);
            end
            % CoSaMP is done
            if  self.LSOnFinalSupport
                % This is used only in high coherence situations.
                % The least squares estimate is not good enough
                % as it allows more distribution of energy on to other 2K indices
                % solve another least squares problem.
                subdict = dict.columns(current_support);
                solution_nz_mat = linsolve(subdict, Y);
            end
            % copy the solution
            result.Z = zeros(d, s);
            result.Z(current_support, :) = solution_nz_mat;
            % copy the measurement residual
            result.R = residual_mat;
            result.iterations = iterations;
            result.support = current_support;
        end

    end

    methods(Static)
        function largest_indices  = largest_thresholding(dict, R, count)
            % We create the proxy of current residual
            E = dict.apply_ctranspose(R);
            tmp = abs(e);
            [~, indices] = sort(tmp, 'descend');
            % We identify the support for largest count entries
            largest_indices  = indices(1:count);
        end

        function largest_indices = largest_k_l1(data_matrix, K)
            sums = spx.norm.norms_l1_rw(data_matrix);
            [~, indices] = sort(sums, 'descend');
            % We identify the support for largest K entries
            largest_indices  = indices(1:K);
        end

        function largest_indices = largest_k_l2(data_matrix, K)
            sums = spx.norm.norms_l2_rw(data_matrix);
            [~, indices] = sort(sums, 'descend');
            % We identify the support for largest K entries
            largest_indices  = indices(1:K);
        end

        function largest_indices  = largest_mp(dict, r, count)
            % Uses matching pursuit to identify largest indices.
            % We create the proxy of current residual
            d = size(dict,2);
            indices = false(1, d);
            k = 0;
            while k < count
                products = dict.apply_ctranspose(r);
                abs_products = abs(products);
                % Find the highest inner product
                [~, index] = max(abs_products);
                % Add this index to support
                if ~indices(index)
                    % we have discovered a new atom.
                    k = k + 1;
                    indices(index) = true;
                end
                % pick up the coefficient of this inner product
                coeff = products(index);
                % update residual
                r = r - coeff * dict.column(index);
            end
            % We identify the support for largest 2K entries
            largest_indices  = find(indices);
        end

        function X = least_squares_on_support(subdict, Y)
            [U,s,V] = csvd(subdict);
            num_signals = size(Y, 2);
            num_cols = size(subdict, 2);
            X = zeros(num_cols, num_signals);
            x = zeros(num_cols, 1);
            for ns=1:num_signals
                y = Y(:, ns);
                normBound = 2*norm(y);
                [alpha2,lambda,maxIterReached] = lsqi(U,s,V,y,normBound);
                if maxIterReached
                    % lsqi did not converge; use cvx (slower) instead
                    [aa,bb] = size(subdict);
                    cvx_begin quiet
                      variable alpha3(bb) complex;
                      minimize norm(y - subdict*alpha3)
                      subject to
                      norm(alpha3) <= normBound;
                    cvx_end
                    x = alpha3;
                else
                    x = alpha2;
                end
                X(:, ns) = x;
            end
        end
    end
end

author: indigits

repo: sparse-plex


text:
function check_WC_params(par)

if isfield(par,'features') && strcmp('features','wav')

    if isfield(par,'scales')&&  isfield(par,'par.w_post') &&  isfield(par,'par.w_pre')
         L = wmaxlev(par.w_pre+par.w_post,'haar');
         if L<scale
             error('[par.scale] exceeds maximum wavelet decomposition level for a waveform of length [par.w_pre+par.w_post]')
         end
    end
    if isfield(par,'par.w_post') &&  isfield(par,'par.w_pre')
        x = log2(par.w_pre+par.w_post);
        if floor(x)~=x
            error('The length [par.w_pre+par.w_post] should be a power of 2')
        end
    end
end
end

author: csn-le


text:
% ir_mex_build_mri
% run matlab's "mex" command to "compile" the mri-related code
% into mex files.

dir_current = pwd;
dir_mri = path_find_dir('mri');
cd(dir_mri)

% compile mex files using ir_mex_fun
fun = ir_mex_fun();
fun('exp_xform_mex.c')
fun('mri_exp_mult_mex.c')

cd(dir_current)

% "test" mex files by running them (each should display usage)
try
	exp_xform_mex
	mri_exp_mult_mex
catch
    disp(['WARNING! Error compiling mex files in ', mfilename])
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
%DEFINITIONS_POLYGONS  One-line description here, please.
%
%   output = definitions_polygons(input)
%
%   Example
%   definitions_polygons
%
%   See also
%

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% INRAE - BIA Research Unit - BIBS Platform (Nantes)
% Created: 2022-02-14,    using Matlab 9.10.0.1739362 (R2021a) Update 5
% Copyright 2022 INRAE.

% create polygon
poly = [

% set up display
figure(1); clf; hold on;
axis equal; axis([20 230 70 230]);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
% SEGMENT_GRAPH Given a sparse, square graph of edges weights segment the nodes
% of the graph into connected sub-components using the greedy merge-based
% method of "Graph Based Image Segmentation".
%
% C = segment_graph(A)
% C = segment_graph(A,'ParameterName',ParameterValue, ...)
%
% Inputs:
%   A  #A by #A sparse, square matrix of edge weights
%   Optional:
%     'Threshold' followed by "C" threshold to use (paper writes that this
%       roughly corresponds to minimum size, though it's really just adding a
%       weight of size/C to components. In any case, increasing this will tend
%       to produce larger segments.
%     'MinSize' followed by the minimum size of an output component. This
%       constraint is enforced as a _post process_.
% Output:
%   C  #A by 1 list of component ids
%
% Example:
%   [V,F] = load_mesh('~/Dropbox/models/Cosmic blobs/Model9.off');
%   A = adjacency_dihedral_angle_matrix(V,F);
%   [AI,AJ,AV] = find(A);
%   A = sparse(AI,AJ,exp(abs(pi-abs(AV-pi))),size(A,1),size(A,2));
%   L = -(A - diag(sum(A,2)));
%   C = segment_graph(L,'Threshold',500,'MinSize',20);
%   tsurf(F,V,'CData',C);
%   colormap(cbrewer('Set1',(max(C))));
%   view(2);
%   axis equal;

author: alecjacobson

repo: gptoolbox


text:
function [QC, dat] = canlab_qc_metrics1(epi_names, mask_name, varargin)
% Calculate quality control metrics on a 4-D EPI Analyze or Nifti file
%
% Standard CANlab usage is to enter a single 4-D ravol* for one run, and
% the brain mask implicit_mask.img created in canlab_preproc.m
%
% :Inputs:
%
%   **epi_names:**
%        Names of (usually one) 4-D EPI file, in cell or string, full path
%
%   **mask_name:**
%        Name of brain mask image, string, full path
%              IF EMPTY: Uses implict masking (better) and calculates
%              ghost/signal outside mask
%
% :Optional Inputs:
%
%   **noplot:**
%        skip plots
%
%   **noverbose:**
%        skip output printout to screen
%
%   **printfile:**
%        followed by name of file to print output to, full path
%
%   **noheader:**
%        suppress printing of header (var names) in output
%
% :Missing values and basic info:
%
%   **num_images:**
%        number of images
%
%   **missing_vox:**
%        Voxels in mask with NaN values or zero values at every time point
%
%   **missing_images:**
%        Images with NaN values or zero values at every voxel
%
%
%   **missing_values:**
%        NaN or zero values in valid images / voxels. Zeros could
%        be interpreted as values of zero in analysis, causing artifacts in results
%        if these are actually invalid values.
%
% Missing voxels will often be ignored in analyses in most software, but
% Missing images/values could cause problems
%
% :Basic signal to noise:
%
%
%   **perc-mean-ghost:**
%        mean signal outside the mask / mean total signal
%
%
%   **mean_snr:**
%        mean Cohen's d (signal/noise, SNR) across time (temporal SNR) within the mask.
%    Mean divided by standard deviation across time at each voxel, averaged.  Higher is better.
%
%
%   **snr_inhomogeneity:**
%        standard deviation of SNR within the mask. Lower is better.
%
%
%   **snr_inhomogeneity95:**
%        95% confidence range for SNR within the mask. Lower is better.
%
%
%   **rms_successive_diffs:**
%        Essentially a high-pass filtered version of SNR,
%        expressed as a fraction of the overall mean and averaged across voxels. Lower is better.
%
%
%   **rms_successive_diffs_inhomogeneity:**
%        standard deviation of the above across voxels. Lower is better.
%
% :Left-right asymmetry:
%
%
%   **signal_rms_asymmetry:**
%        Voxel-wise left/right root mean square asymmetry in mean signal across time, expressed as
%        a fraction of the mean SNR.  Reflects both gross inhomogeneity and noise.  Lower is better.
%
%
%   **signal_hemispheric_asymmetry:**
%        Root mean square difference between left and
%        right hemispheres, expressed as a fraction of the grand mean signal across time.  Reflects
%        gross inhomogeneity.  Lower is better.
%
%
%   **snr_rms_asymmetry:**
%        Voxel-wise left/right root mean square asymmetry in SNR, expressed as
%        a fraction of the mean SNR.  Reflects both gross inhomogeneity and noise.  Lower is better.
%
%
%   **snr_hemispheric_asymmetry:**
%        Root mean square difference between left and
%        right hemispheres, expressed as a fraction of the mean SNR.  Reflects
%        gross inhomogeneity.  Lower is better.
%
% :Examples:
% ::
%
%    %SETUP
%    mydir{1} = '/Users/tor/Documents/Tor_Documents/Coursework_and_Teaching/psyc7215/Data/UM_Face_House/060518mw/Functional/Preprocessed/run_01';
%    wcard = 'rarun*img';
%    epi_names = filenames(fullfile(mydir{1}, wcard), 'absolute');
%
%    maskdir = '/Users/tor/Documents/Tor_Documents/Coursework_and_Teaching/psyc7215/Data/UM_Face_House/060518mw';
%    mask = 'implicit_mask.img';
%    mask_name = fullfile(maskdir, maskname);
%
%    %RUN
%    QC = canlab_qc_metrics1(epi_names, mask_name);
%
%    QC = canlab_qc_metrics1(epi_names, mask_name, 'noplot', 'printfile', 'test_qc.txt');
%    QC = canlab_qc_metrics1(epi_names, mask_name, 'noplot', 'printfile', 'test_qc.txt', 'noheader');
%

doplot = 1;
verbose = 1;
QC = [];
printfile = [];
noheader = 0;

for i = 1:length(varargin)

    if ischar(varargin{i})

        switch varargin{i}
            % functional commands
            case 'noplot', doplot = 0;
            case 'noverbose', doverbose = 0;
            case 'noheader', noheader = 1;
            case 'printfile', printfile = varargin{i+1}; varargin{i+1} = [];

            otherwise, warning(['Unknown input string option:' varargin{i}]);
        end
    end
end

if iscell(epi_names), epi_names = strvcat(epi_names{:}); end

% idstr = epi_names';
% idstr = idstr(:);
% if length(idstr) > 100, idstr = idstr(1:100); end
idstr = epi_names;
if size(idstr,2) > 50, idstr = spm_file(idstr(1,:),'short25'); end

QC = struct('filenames', epi_names(:), 'idstr', idstr);
%all_outputs = {'filenames'};
all_outputs = {};

% Get implicit mask, if no mask is entered
if nargin < 2 || isempty(mask_name)

    mask_name = fmri_mask_image(epi_names); % full mask with all valid voxels

%     maskobj = fmri_mask_image(epi_names, 'implicit');
%     indx = maskobj.volInfo.image_indx;
%     indx(maskobj.volInfo.wh_inmask(maskobj.removed_voxels)) = false;

end

% Load data
% ------------------------------------------
dat = fmri_data(epi_names, mask_name);

% Implicit mask
% ------------------------------------------
if isa(mask_name, 'image_vector')

    fprintf('Calculating implicit mask: ');
    m = mean(dat);
    [dummy, dummy, nvox, is_inmask] = fmri_mask_thresh_canlab(m); %to use dip-based mask: fmri_mask_thresh_canlab(m, [],'dip');
    fprintf('%3.0f voxels in, %3.0f voxels out of mask.\n', nvox, sum(~is_inmask));

    out_of_mask_dat = dat.dat(~is_inmask, :);

    dat = remove_empty(dat, ~is_inmask, []);

    perc_mean_ghost = mean(out_of_mask_dat(:)) ./ mean(dat.dat(:));

end

% do this here; we will use in plot
snr = mean(dat.dat') ./ std(dat.dat');

% Plot data
% ------------------------------------------
if doplot
    plot(dat);

    h = findobj('Tag', 'fmri data matrix');
    figure(h);
    subplot(2, 3, 6);
    [h, x] = hist(snr, 100);
    h = h ./ sum(h);
    plot(x, h, 'k', 'LineWidth', 3);

    title('SNR distribution (PDF)');
end

% Missing values
% ------------------------------------------

outnames = {'num_images', 'vox_vol', 'missing_vox', 'missing_images', 'missing_values'};
if exist('perc_mean_ghost', 'var'), outnames = [outnames 'perc_mean_ghost']; end

all_outputs = [all_outputs outnames];

dat = remove_empty(dat);

vox_dims = diag(dat.volInfo.mat(1:3, 1:3));
vox_vol = prod(abs(vox_dims));
num_images = size(dat.removed_images, 1);
missing_vox = sum(dat.removed_voxels);
missing_images = sum(dat.removed_images);

missing_values = sum(all(isnan(dat.dat) | dat.dat == 0, 2));

for i = 1:length(outnames)

    eval(['myvar = ' outnames{i} ';']);
    QC.(outnames{i}) = myvar;

end

% SNR
% ------------------------------------------
outnames = {'mean_snr', 'mean_snr_per_mm', 'snr_inhomogeneity', 'snr_inhomogeneity95' 'rms_successive_diffs' 'rms_successive_diffs_inhomogeneity'};
all_outputs = [all_outputs outnames];

%snr = mean(dat.dat') ./ std(dat.dat');
mean_snr = mean(snr);
mean_snr_per_mm = mean_snr ./ vox_vol;

snr_inhomogeneity = std(snr);
snr_inhomogeneity95 = prctile(snr, 95) - prctile(snr, 5);

rmssd = ( diff(dat.dat') .^ 2 ./ size(dat.dat, 2) ) .^ .5; % rms at each voxel
rms_successive_diffs = mean(rmssd);

rms_successive_diffs_inhomogeneity = std(rmssd);

for i = 1:length(outnames)

    eval(['myvar = ' outnames{i} ';']);
    QC.(outnames{i}) = myvar;

end

% Asymmetry calculation
% ------------------------------------------
outnames = {'signal_rms_asymmetry', 'signal_hemispheric_asymmetry', 'snr_rms_asymmetry', 'snr_hemispheric_asymmetry'};
all_outputs = [all_outputs outnames];

dat = replace_empty(dat);
m = mean(dat.dat');
snr = m ./ std(dat.dat');

[signal_rms_asymmetry, signal_hemispheric_asymmetry] = calculate_asymmetry(dat, m);

[snr_rms_asymmetry, snr_hemispheric_asymmetry] = calculate_asymmetry(dat, snr);

for i = 1:length(outnames)

    eval(['myvar = ' outnames{i} ';']);
    QC.(outnames{i}) = myvar;

end

% Print
% ------------------------------------------
for i = 1:length(all_outputs)
    out(1, i) = mean(QC.(all_outputs{i}));
end

if ~isempty(printfile)
    diary(printfile);
end

if verbose || printfile
    if noheader
        % no header row
        print_matrix(out, [], {QC.idstr});
    else
        print_matrix(out, all_outputs, {QC.idstr})
    end

end

if ~isempty(printfile)
    diary off
end

end % main function

function [rms_asymmetry, hemispheric_asymmetry] = calculate_asymmetry(dat, snr)

xyz = round(voxel2mm(dat.volInfo.xyzlist', dat.volInfo.mat)');
isleft = xyz(:, 1) < 0;
isright = xyz(:, 1) > 0;
xyz(:, 1) = abs(xyz(:, 1));
[u1, indx1] = unique(xyz, 'first', 'rows');
[u2, indx2] = unique(xyz, 'last', 'rows');

if any(any(u1 - u2)), disp('Bug in asymmetry algorithm'); end
wh_unique = indx1 ~= indx2; % these have two distinct voxels with the same mm coordinates after rounding and abs x

lrdat = [snr(indx1(wh_unique)); snr(indx2(wh_unique))];

rms_asymmetry = sqrt( nansum(diff(lrdat) .^ 2) ./ size(lrdat, 2) ) ./ nanmean(lrdat(:));

tmpsnr = snr;
tmpsnr(~isleft) = NaN;
leftsnr = tmpsnr([indx1(wh_unique); indx2(wh_unique)]);
leftsnr = nanmean(leftsnr);

tmpsnr = snr;
tmpsnr(~isright) = NaN;
rightsnr = tmpsnr([indx1(wh_unique); indx2(wh_unique)]);
rightsnr = nanmean(rightsnr);

hemispheric_asymmetry = sqrt((leftsnr - rightsnr)^2 ./ 2) ./ sum([leftsnr rightsnr]);

end


text:
classdef EulerAnglesDataClass < TimeSeriesDataBaseClass

    %% Public 'read-only' properties
    properties (SetAccess = private)
        FileNameAppendage = '_EulerAngles.csv';
        Phi = [];
        Theta = [];
        Psi = [];
    end

    %% Public methods
    methods (Access = public)
        function obj = EulerAnglesDataClass(varargin)
            fileNamePrefix = varargin{1};
            for i = 2:2:nargin
                if  strcmp(varargin{i}, 'SampleRate'), obj.SampleRate = varargin{i+1};
                else error('Invalid argument.');
                end
            end
            data = obj.ImportCSVnumeric(fileNamePrefix);
            obj.Phi = data(:,2);
            obj.Theta = data(:,3);
            obj.Psi = data(:,4);
            obj.SampleRate = obj.SampleRate;    % call set method to create time vector
        end
        function fig = Plot(obj)
            if(obj.NumPackets == 0)
                error('No data to plot.');
            else
                if(isempty(obj.Time))
                    time = 1:obj.NumPackets;
                else
                    time = obj.Time;
                end
                fig = figure('Name', obj.CreateFigName());
                hold on;
                plot(time, obj.Phi, 'r');
                plot(time, obj.Theta, 'g');
                plot(time, obj.Psi, 'b');
                title('Euler angles');
                xlabel(obj.TimeAxis);
                ylabel('Angle (degrees)');
                legend('\phi', '\theta', '\psi');
                hold off;
            end
        end
    end
end

author: xioTechnologies

repo: Gait-Tracking-With-x-IMU


text:
function pass = test_coeffs_vals( pref )

% Grab some preferences
if ( nargin == 0 )
    pref = chebfunpref();
end
tol = 1e3*pref.techPrefs.chebfuneps;
j = 1;
%grab some coeffs
u = randnfun2(.3);
v = randnfun2(.4);
f = [u;v];

% test coeffs2:
[x, y] = coeffs2(f);
pass(j) = norm(chebfun2(x, 'coeffs')-u, inf) < tol;
pass(j+1) = norm(chebfun2(y, 'coeffs')-v, inf) < tol;
j = j+2;

%parameters
[x, y] = coeffs2(f, 50, 60);
pass(j) = norm(x-coeffs2(u, 50, 60), inf) < tol;
pass(j+1) = norm(y-coeffs2(v, 50, 60), inf) < tol;
j = j+2;

% test coeffs2chebfun2v:
f2 = chebfun2v.coeffs2chebfun2v(coeffs2(u), coeffs2(v));
pass(j) = norm(f - f2) < tol;
j = j+1;

%test coeffs2vals:
[x, y] = coeffs2(f2);
[u, v] = chebfun2v.coeffs2vals(x,y);
pass(j) = norm(chebfun2.coeffs2vals(x)-u, 'inf') < tol;
pass(j+1) = norm(chebfun2.coeffs2vals(y)-v, 'inf')< tol;
j= j+2;

%testvals2coeffs:
[a,b] = chebfun2v.vals2coeffs(u,v);
pass(j) = norm(x - a, 'inf');
pass(j+1) = norm(y - b, 'inf');

if (nargout > 0)
    pass = all(pass(:));
end
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [I,J,K]=MRcart2imNeg(Xc,Yc,Zc,v,OR,rn,cn)

%--------------------------------------------------------------------------
% function [I,J,K]=MRcart2im(X,Y,Z,v,OR,r,c)
%
% This function calculates the image coordinates I,J,K of the cartesian
% coordinates X,Y,Z based on the voxel dimensions v, the origin OR (cartesian
% coorindates of the point I=1, J=1, K=1), the offsets r anc c.
%
% v=N_info.PixelSpacing;
% v(3)=N_info.SliceThickness;
% OR=N_info.ImagePositionPatient;
% r=N_info.ImageOrientationPatient(4:6);
% c=N_info.ImageOrientationPatient(1:3);
%
% Based on:
% http://cmic.cs.ucl.ac.uk/fileadmin/cmic/Documents/DavidAtkinson/DICOM_6up.pdf
% and the file "TranTransform matrix between two dicom image coordinates"
% from the matlab central file exchange written by Alper Yaman.
%--------------------------------------------------------------------------

reForm=~isvector(Xc);

%Switch convention
X=Yc; Y=Xc; Z=Zc;

size_M=size(Xc);

% r => row direction vector
r=vecnormalize(cn);

% c => column direction vector
c=vecnormalize(rn);

%Determine s => slice direction vector N.B. defined as cross(c,r) not
%cross(r,c)
s=cross(c',r'); %Determine s => slice direction vector
s=vecnormalize(s);

X=X(:)';Y=Y(:)';Z=Z(:)';

%Translation
T  = [1 0 0 OR(1);...
      0 1 0 OR(2);...
      0 0 1 OR(3);...
      0 0 0 1];
%Rotation
R  = [r(1) c(1) s(1) 0;...
      r(2) c(2) s(2) 0;...
      r(3) c(3) s(3) 0;...
      0    0    0    1];
%Scaling
S  = [v(1) 0    0    0;...
      0    v(2) 0    0;...
      0    0    v(3) 0;...
      0    0    0    1];

T0 = eye(size(T));

M  = T * R * S * T0; %The transformation matrix

XYZ=ones(4,length(X));
XYZ(1,:)=X;
XYZ(2,:)=Y;
XYZ(3,:)=Z;

IJK=(M\XYZ)';

I=IJK(:,1)+1; J=IJK(:,2)+1; K=IJK(:,3)+1;

if reForm %If the input is not a vector reshape
    I=reshape(I,size_M);
    J=reshape(J,size_M);
    K=reshape(K,size_M);
end

end

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function n = mep2dep(n)
   % convert a matlab date to an epoch date
   n = (n - 719529) * 86400;
end

author: geoscience-community-codes

repo: GISMO


text:
%ASAGLOB_AR_COND ARMA/MA to AR order selection conditioning
%   Declaring the variable ASAGLOB_AR_COND in the base workspace, for
%   example by typing at the command prompt,
%     ASAglob_ar_cond = 1;
%   will enable the ARMASA function ARMASEL to read out the value of
%   ASAGLOB_AR_COND. ASAGLOB_AR_COND controls whether the MA and ARMA
%   estimation routines, used in ARMASEL, will condition the selection of
%   AR orders to the order that has been previously selected by SIG2AR.
%   More specifically, by choosing a confined set of candidate AR orders
%   and assigning ASAGLOB_AR_COND = 1, orders of internally used AR models
%   in MA and ARMA estimation will be determined on the basis of the
%   specific selection of an AR model order from that confined set. By
%   leaving ASAGLOB_AR_COND unspecified or assigning ASAGLOB_AR_COND = 0,
%   a standard AR order selection procedure is repeated if necessary, to
%   supply the MA and ARMA estimators with that order.

disp('  Variable ASAGLOB_AR_COND is currently undefined.')
disp('  Type ''help ASAglob_ar_cond'' for more information.')


text:
function [Eft, Varft, lpyt, Eyt, Varyt] = gpia_loopred(gp_array, x, y, varargin)
%GPIA_LOOPRED  Leave-one-out predictions with Gaussian Process IA approximation
%
%  Description
%    [EFT, VARFT, LPYT, EYT, VARYT] = GPIA_LOOPRED(RECGP, X, Y)
%    takes a cell array of GP structures together with matrix X of
%    input vectors, matrix X of training inputs and vector Y of
%    training targets, and evaluates the leave-one-out predictive
%    distribution at inputs X and returns means EFT and variances
%    VARFT of latent variables, the logarithm of the predictive
%    densities PYT, and the predictive means EYT and variances
%    VARYT of observations at input locations X.
%
%    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.
%      is     - defines if importance sampling weighting is used 'on'
%               (default). If set to 'off', integration points and
%               weights for the full posterior are used.
%
%    Given Gaussian likelihood or non-Gaussian likelihood and
%    latent method Laplace or EP, LOO-posterior given
%    hyperparameters is computed analytically or with analytic
%    approximation and LOO-posterior of the hyperparameters is
%    approximated using importance sampling reweighted integration
%    points from gp_ia. Optionally integration weights for full data
%    posterior of hyperparameters can be used (by setting option
%    'is' to 'off').
%
%  References:
%    Aki Vehtari and Jouko Lampinen (2002). Bayesian model
%    assessment and comparison using cross-validation predictive
%    densities. Neural Computation, 14(10):2439-2468.
%
%  See also
%   GP_LOOPRED, GP_IA, GP_PRED
%
% Copyright (c) 2009 Ville Pietil�inen
% Copyright (c) 2009-2010 Jarno Vanhatalo
% Copyright (c) 2010,2012 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.

  ip=inputParser;
  ip.FunctionName = 'GPIA_LOOPRED';
  ip.addRequired('gp_array',@iscell);
  ip.addRequired('x', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
  ip.addRequired('y', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
  ip.addParamValue('z', [], @(x) isreal(x) && all(isfinite(x(:))))
  ip.addParamValue('is', 'on', @(x) ismember(x,{'on' 'off'}))
  ip.parse(gp_array, x, y, varargin{:});
  z=ip.Results.z;
  is=ip.Results.is;

  nGP = numel(gp_array);
  n=size(x,1);

  P_TH = zeros(1,nGP);
  Efts = zeros(n,nGP);
  Varfts = zeros(n,nGP);
  lpyts = zeros(n,nGP);
  if nargout > 3
    Eyts = zeros(n,nGP);
    Varyts = zeros(n,nGP);
  end

  for i1=1:nGP
    Gp=gp_array{i1};
    P_TH(1,i1) = Gp.ia_weight;
    % compute leave-one-out predictions for each hyperparameter sample
    % if latent method is MCMC, then these samples are for latent values, too
    if nargout <= 3
      [Efts(:,i1), Varfts(:,i1), lpyts(:,i1)] = gp_loopred(Gp, x, y, 'z', z);
    else
      [Efts(:,i1), Varfts(:,i1), lpyts(:,i1), Eyts_temp, Varyts_temp] = gp_loopred(Gp, x, y, 'z', z);
      if ~isempty(Eyts_temp)
        Eyts(:,i1) = Eyts_temp;
      end
      if ~isempty(Varyts_temp)
        Varyts(:,i1) = Varyts_temp;
      end
    end
  end

  if isequal(is,'off')
    P_TH=repmat(P_TH,n,1);
  else
    % log importance sampling weights
    lw=-lpyts;
    % normalize weights
    for i2=1:n
      % this works even when lw have large magnitudes
      lw(i2,:)=lw(i2,:)-sumlogs(lw(i2,:));
    end
    % importance sampling weights
    w=exp(lw);
    % reweight ia weights
    P_TH=bsxfun(@times,P_TH,w);
    P_TH=bsxfun(@rdivide,P_TH,sum(P_TH,2));
    % check the effective sample size
    m_eff=1./sum(P_TH.^2,2);
    if min(m_eff)<nGP/5
      warning(sprintf('For %d folds the effective sample size in IS is less than m/5',sum(m_eff<(nGP/5))))
    end
    %fprintf('nGP=%d, min(neff)=%.0f, min(neff)/nGP=%.2f\n',nGP, min(1./sum(P_TH.^2,2)), min(1./sum(P_TH.^2,2))./nGP)
    % PSIS
    if nGP>=200
        % (new) default sample size for is_normal and is_t is 200
        % CCD with nParam>=12 has at least 281 points
        [lw,pk] = psislw(log(P_TH'),10);
        P_TH=exp(lw');
        % check whether the variance and mean of the raw importance ratios is finite
        % PSIS weights have always finite variance and mean, but if raw importance
        % ratios have infinite variance the convergence to true value is
        % slower and if raw importance ratios have non-existing mean the the
        % estimate can't converge to true value
        vkn1=sum(pk>=0.5&pk<0.7);
        vkn2=sum(pk>=0.7&pk<1);
        vkn3=sum(pk>=1);
        n=numel(pk);
        if vkn1>0
            warning('%d (%.0f%%) PSIS Pareto k estimates between 0.5 and .7',vkn1,vkn1/n*100)
        end
        if vkn2>0
            warning('%d (%.0f%%) PSIS Pareto k estimates between 0.7 and 1',vkn2,vkn2/n*100)
        end
        if vkn3>0
            warning('%d (%.0f%%) PSIS Pareto k estimates greater than 1',vkn3,vkn3/n*100)
        end
    end
  end

  % compute combined predictions
  Eft = sum(Efts.*P_TH, 2);
  Varft = sum(Varfts.*P_TH,2) + sum(bsxfun(@minus,Efts,Eft).^2.*P_TH,2);
  if nargout > 2
    lpyt = log(sum(exp(lpyts+log(P_TH)),2));
    if nargout > 3
      Eyt = sum(Eyts.*P_TH,2);
      Varyt = sum(Varyts.*P_TH,2) + sum(bsxfun(@minus,Eyts,Eyt).^2.*P_TH, 2);
    end
  end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

function phiJump = cf_lognormjump(u,a,b,lambda,T)
% log characteristic function for lognormal jumps
    phiJump = lambda*T*(-a*u*1i + (exp(u*1i*log(1.0+a) ...
        +0.5*b*b*u*1i.*(u*1i-1.0))-1.0));
end