text:
% FORMAT:
%
% sampval = ms2sample(timems, fs, rounding, offset)
%
% INPUTS:
%
% timems      - time value in milliseconds
% fs          - sample rate
%
% Optional inputs:
%
% rounding    - 1 means round the result {default}, 0 means do not round
% offset      - offset time in samples (you can use ms2sample recursively here). Default 0
%
% OUTPUT:
%
% sampval     - time in samples
%
% EXAMPLES:
%
% 1) For a time serie recorded from 0 to 5 secs, at fs=500 sps, get the sample index at the time 674 ms.
%
%>> sampval = ms2sample(674, 500)
%
%sampval =
%
%   337
%
% 2) For a time serie recorded from -1 to 5 secs, at fs=500 sps, get the sample index at the time 674 ms.
%
%>> sampval = ms2sample(674, 500, 1, ms2sample(1000, 500))
%
%sampval =
%
%   837
%
%
% Author: Javier Lopez-Calderon
% Center for Mind and Brain
% University of California, Davis,
% Davis, CA
% 2013

function sampval = ms2sample(timems, fs, rounding, offset)
if nargin<1
        help ms2sample
        return
end
if nargin<4
        offset = 0;
end
if nargin<3
        rounding = 1;
end
if nargin<2
        error('Two inputs are requiered at least.')
end
sampval = offset + timems*fs/1000;
if rounding
        sampval = round(sampval);
end

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
%% Example
% Using iiwa as a plotter, the user can mount:
% 1- A pen, in such a case the robot is going to draw
% 2- A laser, in such a case the robot is going to cut the object

% This script utilizes the real time control of the KST
% where the iiwa7R800 is specified, if you have iiwa14R820 uncomment the
% code-line:
%

% Note that for acheiving a precise drawing do the following steps:
% A- On the Sunrise Code of the KST project:
%         1- Open the file (MatlabToolboxServer.java)
%         2- Locate the function:
%             double getTheDisplacment(double dj)
%             {
%                 double   a=0.07;
%                 double b=a*0.75;
%                 double exponenet=-Math.pow(dj/b, 2);
%                 return Math.signum(dj)*a*(1-Math.exp(exponenet));
%
%             }
%         3- Change it into
%             double getTheDisplacment(double dj)
%             {
%                 return dj;
%             }
% B- Check the pdf file in the folder (Tips and tricks) of the github repo,
% and reconfigure your PC according to the recomendations of the file.
% C- Use a PC with good specifications. The PC used for running this script
% is Corei7-6850k processor @3.6GHZ with  32 G RAM.

% After performing the previous steps, synchronise the Sunrise project of
% the KST into the controller.
% Then start the MatlabToolboxServer application on the KUKA iiwa controller
% Then run this script using Matlab

% Copyright: Mohammad SAFEEA, 26th of April 2018

close all,clear all;clc;
%% Add KST to Matlab path
kst_Path=getTheKSTDirectory(pwd);
addpath(kst_Path);

warning('off')

%% Create the robot object
ip='172.31.1.147'; % The IP of the controller
arg1=KST.LBR7R800; % choose the robot iiwa7R800
% arg1=KST.LBR14R820; % un comment this line if you have iiwa14R820
arg2=KST.Medien_Flansch_Touch_pneumatisch; % choose the type of flange
Tef_flange=eye(4); % transofrm matrix of EEF with respect to flange
Tef_flange(3,4)=0.293; % length of the tep of the pen from the surface of the flange
iiwa=KST(ip,arg1,arg2,Tef_flange); % create the object

%% Start a connection with the server
flag=iiwa.net_establishConnection();
if flag==0
  disp('Error unable to connected to the controller');
  return;
else
  disp('Successfully connected to the controller');
end
pause(0.5);

%% Go to an initial configuration
jPos={0,0,0,-pi/2,0,pi/2,0};
relVel=0.15;
iiwa.movePTPJointSpace(jPos, relVel); % move to initial configuration
%% Get Cartesian position of EEF (homing position)
fprintf('Current cartesian position:')
eefpos=iiwa.getEEFPos();
X=eefpos{1};
Y=eefpos{2};
Z=eefpos{3};
eefposDist=eefpos;
disp(eefpos);
%% Motion parameters
deltaZ=5; % Small clearance from the plotting-surface (mm)
z_sheet= 81.2; % the Z coordinates of the plotting-surface (mm)
zUp=z_sheet+deltaZ;
vFast=25/1000; % fast motion velocity, m/sec
vSlow=20/1000; % slow motion velocity, m/sec
%% Read the Cad file
fileName='kst.plt'; % this is the name of the cad file,
% in case you have the cad file in another folder the user can use the
% total path of the file, the file shall be in plt format
[plotFlag,corArray]=loadPltFileFun(fileName); % load the file
%% Move robot down to a posiiton close to the sheet
% fast motion
eefposDist{1}=corArray(1,1)+X;
eefposDist{2}=corArray(2,1)+Y;
eefposDist{3}=zUp;
iiwa.movePTPLineEEF(eefposDist, 100);
%% Start direct servo in Joint space
% get joints angles of robot
jPos  = iiwa.getJointsPos();
% calculate current position of flange point of the robot
qs=zeros(7,1);
for i=1:7
    qs(i)=jPos{i};
end

TefTool=eye(4);
[T0,~]=iiwa.directKinematics(qs); % Transformation matrix at the tip of the pin (TCP)
x=T0(1,4);
y=T0(2,4);

zUp1=T0(3,4);
zDown1=T0(3,4)-deltaZ/1000;

Tt=T0;
n=max(size(corArray));
%% scale and displace the drawing in the plane
scaleX=5;
scaleY=10;
dispY=-1.7;
dispX=-0.22;
corArray(1,:)=corArray(1,:)*scaleX+dispX;
corArray(2,:)=corArray(2,:)*scaleY+dispY;

for i=1:n
    corArray(1,i)=corArray(1,i)+x;
    corArray(2,i)=corArray(2,i)+y;
end
% start the direct servo
 iiwa.realTime_startDirectServoJoints();
disp('Starting direct servo joint space');
%% Start control loop
tic;
vz=0.003;
for i=1:n-1
    if(plotFlag(i)==0)
        [Tt,qs]=goDown(iiwa,qs,TefTool,vz,zDown1);
        v=vSlow;
    else
        [Tt,qs]=goUp(iiwa,qs,TefTool,vz,zUp1);
        v=vFast;
    end
    x0=x;
    x1=corArray(1,i);
    y0=y;
    y1=corArray(2,i);
    wattar=((x0-x1)*(x0-x1)+(y0-y1)*(y0-y1))^0.5;
    if wattar==0
        continue;
    end
    vx=v*(x1-x0)/wattar;
    vy=v*(y1-y0)/wattar;
    trajectoryTime0=toc;
    transmittionTime0=trajectoryTime0;
    motionFlag=1;
    while motionFlag
        deltaT=toc-trajectoryTime0;
        x=x0+vx*deltaT;
        y=y0+vy*deltaT;
        %
        wattar1=((x0-x)*(x0-x)+(y0-y)*(y0-y))^0.5;
        if wattar1>wattar
            motionFlag=0;
        end
        if(toc-transmittionTime0)>0.003
            Tt(1,4)=x;
            Tt(2,4)=y;
            [ qs ] = iiwa.gen_InverseKinematics( qs, Tt, 10,0.1 );
            transmittionTime0=toc;
            for k=1:7
                jPos{k}=qs(k);
            end
            iiwa.sendJointsPositionsf(jPos);
        end

    end
end
[Tt,qs]=goUp(iiwa,qs,TefTool,vz,zUp1);
%% Turn off direct servo control
iiwa.realTime_stopDirectServoJoints();
%% Move robot back to home position
% fast motion
eefposDist{1}=X;
eefposDist{2}=Y;
eefposDist{3}=Z;
iiwa.movePTPLineEEF(eefposDist, 100);

iiwa.net_turnOffServer();

repo: KST-Kuka-Sunrise-Toolbox


text:
function [L,dL] = vpbndloss(theta,vp,thetabnd,TolCon)
%VPLOSS Variational parameter loss function for soft optimization bounds.

compute_grad = nargout > 1;     % Compute gradient only if requested

K = vp.K;
D = vp.D;

if vp.optimize_mu
    mu = theta(1:K*D);
    idx_start = K*D;
else
    mu = vp.mu(:)';
    idx_start = 0;
end
if vp.optimize_sigma
    lnsigma = theta(idx_start+(1:K));
    idx_start = idx_start + K;
else
    lnsigma = log(vp.sigma(:));
end
if vp.optimize_lambda
    lnlambda = theta(idx_start+(1:D));
else
    lnlambda = log(vp.lambda(:));
end
if vp.optimize_weights
    eta = theta(end-K+1:end);
else
    eta = [];
end

lnscale = bsxfun(@plus,lnsigma(:)',lnlambda(:));
theta_ext = [];
if vp.optimize_mu; theta_ext = [theta_ext; mu(:)]; end
if vp.optimize_sigma || vp.optimize_lambda; theta_ext = [theta_ext; lnscale(:)]; end
if vp.optimize_weights; theta_ext = [theta_ext; eta(:)]; end

if compute_grad
    [L,dL] = softbndloss(theta_ext,thetabnd.lb(:),thetabnd.ub(:),TolCon);
    if vp.optimize_mu
        dmu = dL(1:D*K);
        idx_start = D*K;
    else
        dmu = [];
        idx_start = 0;
    end
    if vp.optimize_sigma || vp.optimize_lambda
        dlnscale = reshape(dL((1:D*K)+idx_start),[D,K]);
        if vp.optimize_sigma
            dsigma = sum(dlnscale,1);
        else
            dsigma = [];
        end
        if vp.optimize_lambda
            dlambda = sum(dlnscale,2);
        else
            dlambda = [];
        end
    else
        dsigma = []; dlambda = [];
    end
    if vp.optimize_weights
        deta = dL(end-K+1:end);
    else
        deta = [];
    end
    dL = [dmu(:); dsigma(:); dlambda(:); deta(:)];
else
    L = softbndloss(theta_ext,thetabnd.lb(:),thetabnd.ub(:),TolCon);
end

end

author: acerbilab

repo: vbmc

save_path: github-repos/MATLAB/acerbilab-vbmc


text:
function [pval x K P xx fx pperm] = fexact( varargin )
%FEXACT Fisher's exact test
%
% Fisher's exact test is a statistical test used to compare binary outcomes
% between two groups. For example, a laboratory test might be positive or negative
% and we may be interested in knowing whether there is a difference in frequency of positive
% results among people with a certain disease (cases) compared to those in healthy
% (controls). The test is applicable whenever there are independent tests conducted on two
% groups. A crosstable of results (pos/neg) and status (case/control) produces a 2x2 contingency
% table shown below.
%       cases  controls
%  pos   a        c       K
%  neg   b        d       -
% total  N        -       M
%
% a,b,c,d are the numbers of positive and negative results for the cases
% and controls. K, N, M are the row, column and grand totals.
%
% This function returns a pvalue against the hypothesis that the number of
% positive cases (and by extension other cells in the table) observed or a more extreme
% distribution might have occured by chance given K, N,and M.
%
% usage
%   p = fexact(X,y)
%       y is a vetor of status (1=case/0=control). X is a MxP matrix of
%       binary results (pos=1/neg=0). P can be very large for
%       genotyping assays, each of which can be considered a different way
%       to categorize the cases and controls.
%
%  p = fexact( a,M,K,N, options)
%      M and N are constants described above. a and K are P-vectors. No checks for valid
%      input are made.
%
%  [p a K C] = fexact(X,y);
%       also returns vector a. a(i) is the upper left square of the contingency
%       corresponding to the ith column of X crosstabulated by y. K is a
%       vector containing a+b for the ith column of X.
%       C is a lookup table for CDFs that can be used in subsequent calls
%       to fexact to further improve performance.
%       C( x(i)+1, K(i)+1) is the cdf for the tail specified in
%       options.
%       NB. this lookup table is only for the given M and N values
%
%  p = fexact( a,M,K,N,C]
%      C is a lookup table that improve performance. Intended to be
%      used with permutation testing where 1000s of calls are made.
%      the option "tail" is ignored if C is provided, as it implies
%      whatever tail was used to generate them.
%
%  p = fexact(..., 'options', values)
%      options:
%      tail  l(eft)','r(ight)', 'b(oth)'
%            left tail tests a negative association. p( x<=a, M, K, N)
%            right tail tests a positive association. p( x<=b, M, K, M-N)
%            both (default) is either a positive or negative association.
%            this is the sum from the left and right tails of the
%            distribution including all x where p(x|M,K,N) is less than
%            or equal p(a|M,K,N).
%      perm  Q, where Q is an integer. Permutation testing when X has multiple Columns.
%             to correct for multiple testing. The entire set of tests is repeated Q*repsz
%            times with permuted y variables. The reported p-values are
%            corrected for multiple tests by interpolating into the emprical
%            distribution of the minimum p-value obtained from each Q*repsz rounds .
%            To use this option the X and y calling convention must be used
%            (as opposed to (a,M,K,N,...) convention.
%            and Q must be greater than 1.
%      repsz S, where S is an integer. Used with the PERM option to specify
%            the number of replicates to compute at one time (default=100).
%            Larger S is faster but requires more memory.
%
%
% NOTES and LIMITATIONS:
%      This function is extremely fast when doing multiple tests and is
%      acceptable with a small number of tests. However, it uses large
%      tracks of memory. On my 2Gb home Intel Core2 Quad Core Vista machine
%      this function does 250,000 tests with 100 observations each in 0.10
%      seconds. I run out of memory when I do more
%
% example
%  x = unidrnd(2, 200, 1000)-1;
%  y = unidrnd(2,200,1)-1;
%  p = fexact( x, y );                  % generates p-values for 1000 tests in x
%  p1 = fexact( x, y, 'perm', 10); % reports p-value relative to empirial
%                                               % cdf, since these are randomly generated
%                                               % the values in p1 should be near  1
%

% $Id: fexact.m 642 2012-06-11 18:25:47Z mboedigh $
% Copyright 2012 Mike Boedigheimer
% Amgen Inc.
% Department of Computational Biology
% mboedigh@amgen.com

p = inputParser;
p.addRequired('A');
p.addOptional('y',[]);
p.addOptional('K',[]);
p.addOptional('N',[]);
p.addOptional('C',[]);
p.addOptional('tail', 'b', @(c)ismember( c, {'b','l','r'} ));
p.addOptional('perm',  0, @(x) isnumeric(x)&isscalar(x));
p.addOptional('repsz', 100, @(x) isnumeric(x)&isscalar(x));
p.parse(varargin{:});

P = p.Results.C;
if isempty(p.Results.N)     % using fexact(X,y,options)
    X = p.Results.A;
    y = p.Results.y;
    if (~islogical(y) && ~all( ismember( y, [0 1])) )
        error('linstats:fexact:InvalidArgument', 'y must be in (0,1)' );
    end
    if (~islogical(X) && ~all( ismember( X(:), [0 1])) )
        error('linstats:fexact:InvalidArgument', 'X must be in (0,1)' );
    end

    y = logical(y);
    N = sum(y);           % number of cases
    M = length(y);
    K = sum(X,1)';
    x = X'*sparse(y);  % in unofficial testing using timeit. This was faster than sum(X(y==1,:))
    % and faster than sum(bsxfun(@eq,X,y))
else % using fexact( a,M,K,N ...)
    x = p.Results.A;
    M = p.Results.y;
    N = p.Results.N;
    K = p.Results.K;
end

switch p.Results.tail
    case 'l'; tail = -1;
    case 'r'; tail = 1;
    case 'b'; tail = 2;
end

if N==0 || N==M
    pval = ones(1,length(x));
    return;
end;

if isscalar(M) && isscalar(N) && isscalar(K) && isscalar(x) && ~(p.Results.perm > 1)
    pval = doTest( x, M, K, N, tail );
    pval(pval>1) = 1;   % fix roundoff errors
    return;
end

if isempty(P)
    P = getLookup( M,N,tail);
end

pval = P( sub2ind(size(P), x+1, K+1));

if p.Results.perm > 1
    pperm   = doPerm(X,y,M,K,P,p.Results.perm,p.Results.repsz);
    [fx xx] = mecdf(pperm(:));
    xx(end) = max(pval);
    pval    = interp1( xx, fx, pval);
end

pval(pval>1) = 1;   % fix roundoff errors

function pperm = doPerm(X,y,M,K,P,nperms,repsz)
% returns pperm a matrix with each element representing the smallest p
% value from a complete study where the rows of y have been randomly
% permuted
pperm = ones(repsz,nperms); % each element is from a permutation test. It is the minimum pvalue among all variables (columns of X)

% sub2ind was previously 85% of the execution time.
% this method computes the index directly rather than
% computing x and then calliung sub2ind. It is much
% faster because it uses the fast multiple engine
% precompute K*nrows +1, then add x after it is calculated
X = [ K*size(P,1)+1 X' ];
B = ones(1,repsz);

for permi = 1:nperms
    % generate new random vector y, by shuffling the rows (do not change
    % the counts of M,K or M)
    [i i] = sort( randn(M,repsz) ); % note randn was faster than rand when I checked
    pperm(:,permi) = min( P( X*[B;y(i)] ),[], 1); % find the smallest p-value among all separate tests for each permutation
end

function P = doTest( x, M, N, K, tail)

minx = max(0,N+K-M);
maxx = min(N,K);
if x < minx
    error( 'x is can not be smaller than max(0,N+K-M)' );
end

if x > maxx;
    error( 'x cannot be larger than min(N,K)');
end

%log p(x|M,N,K) = Q(K,x)+Q(M-K,N-x)-Q(M,N) where Q(n,k) = lnchoosek(n,k)
if tail == -1
    a = minx:x;
    logP = lnchoosek( K, a) + lnchoosek(M-K,N-a) - lnchoosek(M,N);
elseif tail == 1
    a = x:min(K,N);
    logP = lnchoosek( K, a) + lnchoosek(M-K,N-a) - lnchoosek(M,N);
else
    a = 0:min(K,N);
    logP = lnchoosek( K, a) + lnchoosek(M-K,N-a) - lnchoosek(M,N);
    pcrit = logP(x+1);
    k = (logP - pcrit) <= M*eps(max(logP));
    logP  = logP(k);
end

P = sum(exp(logP));

function P = getLookup( M, N, tail)

F = gammaln( 1:M+1);  % used to compute factorials. For small problems
% It is overkill to generate all factorials up
% to M, but it takes on .001 seconds to compute
% all of them to 1..10,001, so we can afford
% it. Memory-wise may be more of an issue, and
% a reason to

% lookup value for lnchoosek(a,b) in matrix L for all 1<=b<=a<=M
% to lookup log nchoosek(n,k) use L(n+1,k+1);
[k n ind] = tri2sqind(M);
L = zeros(M);
L(ind) = lnchoosek( n, k, F);
L = blkdiag(0,L');      % pad first row and col with 0s for n=0, k=0

%log p(x|M,N,K) = Q(K,x)+Q(M-K,N-x)-Q(M,N) where Q(n,k) = lnchoosek(n,k)

% create a lookup table for Log P
% to lookup log p(x,|M,N,K) use P(x+1,K+1);
% P initially contains 1s for valid values of x
% x is in [min( 0, N+K-M), ..., min(N,K)]
% min(N,K) is on the diagonal. the first
% column that contains a 0 in the first row will be when the minimum legal
% value of x is 1, so N+K-M == 1. K will equal M-N+1, which is the column
% M-N+2
P = spdiags(ones(N+1,M+1),0:M-N,N+1,M+1);
[Xind Kind] = find(P);
MKind = M-Kind+2;
NXind = N-Xind+2;

q = sub2ind(size(P),Xind,Kind);
P = nan(size(P));
P(q) = exp( L(sub2ind(size(L),Kind,Xind))    + ... Q(K,x)
    L(sub2ind(size(L),MKind,NXind))  + ... Q(M-K, N-x)
    -L(M+1,N+1));                     % ... Q(M,N)
if tail == -1
    P = cumsum(tril(P, M-N));
elseif tail == 1
    P = flipud(triu(P));
    P = cumsum(P);
    R = N-(-1:size(P,1)-2)';
    P = P(R,:);
else
    % This function has been tested against two reference algorithm using all M
    % and N up to 30. The reference algorithms differ in how they find the
    % find the tail on the other side of the distribution.
    [q r] = size(P);
    [P idx] = sort(P,1);     % sort order, called rank, will be transformed to rank below
    ties = [abs(diff(P)) < M*numel(P)*eps(P(1:end-1,:));
        zeros(1,r)];
    R = repmat( (1:q)',1,r);
    R(col2ind(idx)) = R+ties;
    P = cumsum(P,1);          % reuse P to store cumulative P
    P = P(col2ind(R));
end

function q = lnchoosek( n, k, f)
% lnchoosek natural log of nchoosek
%
% usage
%       q = lnchoosek( n, k );
% usage
%       q = lnchoosek( n, k, f); where f is the gammaln(1:max(n,k)+1);
q    = n-k;       % allocate space same size as n
i    = q>0;       % only compute for valid n and k

if nargin < 3
    % q is 0 when n = k because log nchoosek(a,a) = 1 for any a >= 0
    % otherwise it is given by ...
    q(i) = -log(q(i)) - betaln(k(i)+1,q(i));     % otherwise this is it
else
    q(i) =  f(n(i)+1) - f(k(i)+1) - f( q(i)+1 );
end

q(~i) = 0;

function i = col2ind(order,siz)
%col2ind converts column specific 1-based indices to element based indices
%i = col2ind(A, siz)
%   return i, a set of integer indices into a m x n matrix
%   A is a p x n set of column specific indices as you'd get from the second
%   output of sort
%   size is optional.
%        if present, P is set to SIZ(1), otherwise
%        P is set to the size of the first dimension of ORDER,
%
% example
%   [xsort order] = sort(x);
%   wrong = x(order);    % this isn't what you want
%   xsort = x(col2ind(order));   % this is

[m p] = size(order);
if nargin < 2
    q = m;
else
    q = siz(1);
end
oset = 0:q:(q*p-1);
i = order + repmat( oset, m, 1 );

function [i,j,k] = tri2sqind( m, k )
%TRI2SQIND subscript and linear indices for upper tri portion of matrix
%
% get indices into a square matrix for a vector representing a the upper
% triangular portion of a matrix such as those returned by pdist.
%
% [i,j,k] = tri2sqind( m, k )
%  If V is a hypothetical vector representing the upper triangular portion
%  of a matrix (not including the diagonal) and
%  M is the size of a square matrix and
%  K is an optional vector of indices into V then tri2sqind returns
%  (i,j) the subscripted indices into the equivalent square matrix.
%  K is an integer index into the equivalent square matrix
%
% Example
%  X = randn(5, 20);
%  Y = pdist(X, 'euclidean');
%  [i,j,k] = tri2sqind( 5 );
%  S = squareform(Y);
%  isequal( Y(:), S(k) );
%  Z = zeros(5);
%  Z(k) = Y;

max_k = m*(m-1)/2;

if ( nargin < 2 )
    k = (1:max_k)';
end;

if any( k > max_k )
    error('linstats:tri2sqind:InvalidArgument', 'ind2subl:Out of range subscript');
end;

i = floor(m+1/2-sqrt(m^2-m+1/4-2.*(k-1)));
j = k - (i-1).*(m-i/2)+i;
k = sub2ind( [m m], i, j );

function [fx x] = mecdf( y )
% MECDF my ECDF function to avoid dependendencies on stats toolbox
% assumes 0 <= y <= 1
% x always includes 0 and fx(x>1) = 1;
n = length(y);
[x i] = unique(sort(y));
fx = linspace( 0, 1, n+1 )';
fx = fx(i+1);
if x(1)>0
    x = [0;x];
    fx = [0;fx];
end

if x(end) < 1    % set up the maximum x value as Inf, so interp1 never fails
    x  = [x;Inf];
    fx = [fx;1];
end

author: Sable

repo: mcbench-benchmarks


text:
function names_numbered = DisplayMaskNames(Msk_IDs,MASKs)
if Msk_IDs ~= 0,
    names = MASKs.MaskDatabaseNames(Msk_IDs);
    names_numbered = cell(size(names));
    disp('Regions:');
    for i = 1:length(names),
        names_numbered{i} = [num2str(Msk_IDs(i)), ': ', names{i}];
        disp(names_numbered{i});
    end
end
end

author: xiuyechen

repo: FishExplorer


text:
function events = in_events_bids(sFile, EventFile)
% IN_EVENTS_BIDS: Read a BIDS _events.tsv file (columns "onset", "duration", "trial_type").
%
% OUTPUT:
%    - events(i): array of structures with following fields (one structure per event type)
%        |- label   : Identifier of event #i
%        |- samples : Array of unique time indices for event #i in the corresponding raw file
%        |- times   : Array of unique time latencies (in seconds) for event #i in the corresponding raw file

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2019-2022

% Read tsv file
% https://bids-specification.readthedocs.io/en/stable/04-modality-specific-files/05-task-events.html
% trial_type would be a subset of possible events that identify trial category.
% value can be numbers or strings identifying the event.
% Only 'onset' and 'duration' are required.
Markers = in_tsv(EventFile, {'onset', 'duration', 'trial_type', 'channel', 'value'}, 0);
if isempty(Markers) || isempty(Markers{1,1})
    events = [];
    return;
end
% Standardize empty/missing values to empty char.
Markers(cellfun(@(c) (isequal(c,'n/a') || isequal(c,'N/A') || isempty(c)), Markers(:))) = {''};
% If there is no trial_type and no value information: use the filename as the event name
if all(cellfun(@isempty, Markers(:,3)) & cellfun(@isempty, Markers(:,5)))
    [fPath, fbase, fExt] = bst_fileparts(EventFile);
    Markers(:,3) = repmat({fbase}, size(Markers(:,3)));
elseif all(cellfun(@isempty, Markers(:,3)))
    Markers(:,3) = Markers(:,5);
elseif ~all(cellfun(@isempty, Markers(:,5)))
    % Both columns contain data
    for iM = 1:size(Markers, 1)
        if isempty(Markers{iM,3})
            Markers{iM,3} = Markers{iM,5};
        elseif ~isempty(Markers{iM,5})
            % Merge
            Markers{iM,3} = [Markers{iM,3} '-' Markers{iM,5}];
        end
    end
end
% List of events from trial_type and/or value
uniqueEvt = unique(Markers(:,3)');

% Initialize returned structure
events = repmat(db_template('event'), [1, length(uniqueEvt)]);
% Create events list
for iEvt = 1:length(uniqueEvt)
    % Find all the occurrences of event #iEvt
    iMrk = find(strcmpi(Markers(:,3)', uniqueEvt{iEvt}));
    % Get event onsets and durations
    onsets = cellfun(@(c)sscanf(c,'%f',1), Markers(iMrk,1), 'UniformOutput', 0);
    durations = cellfun(@(c)sscanf(c,'%f',1), Markers(iMrk,2), 'UniformOutput', 0);
    channels = Markers(iMrk,4)';
    % Find and reject events with no latency
    iEmpty = find(cellfun(@isempty, onsets));
    if ~isempty(iEmpty)
        iMrk(iEmpty) = [];
        onsets(iEmpty) = [];
        durations(iEmpty) = [];
        channels(iEmpty) = [];
    end
    % Channel names
    for iOcc = 1:length(channels)
        if (isempty(channels{iOcc}) || strcmpi(channels{iOcc}, 'n/a'))
            channels{iOcc} = [];
        else
            channels{iOcc} = {channels{iOcc}};
        end
    end
    % Add event structure
    events(iEvt).label  = uniqueEvt{iEvt};
    events(iEvt).epochs = ones(1, length(iMrk));
    events(iEvt).times  = [onsets{:}];
    % Convert latencies from the first sample in the file to the real timing
    % See specs: https://bids-specification.readthedocs.io/en/stable/04-modality-specific-files/05-task-events.html
    events(iEvt).times = events(iEvt).times + sFile.prop.times(1);
    % Extended events if durations are defined for all the markers
    if all(~cellfun(@isempty, durations)) && all(~cellfun(@(c)isequal(c,0), durations))
        events(iEvt).times(2,:) = events(iEvt).times + [durations{:}];
    end
    % Make the function independent of sFile
    if ~isempty(sFile)
        events(iEvt).times  = round(events(iEvt).times .* sFile.prop.sfreq) ./ sFile.prop.sfreq;
    end
    events(iEvt).reactTimes = [];
    events(iEvt).select     = 1;
    events(iEvt).channels   = channels;
    events(iEvt).notes      = [];
end

author: brainstorm-tools


text:
function b = isPointInCircle(point, circle, varargin)
%ISPOINTINCIRCLE Test if a point is located inside a given circle
%
%   B = isPointInCircle(POINT, CIRCLE)
%   Returns true if point is located inside the circle, i.e. if distance to
%   circle center is lower than the circle radius.
%
%   B = isPointInCircle(POINT, CIRCLE, TOL)
%   Specifies the tolerance value
%
%   Example:
%   isPointInCircle([1 0], [0 0 1])
%   isPointInCircle([0 0], [0 0 1])
%   returns true, whereas
%   isPointInCircle([1 1], [0 0 1])
%   return false
%
%   See also:
%   circles2d, isPointOnCircle
%
%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 07/04/2004.
%

%   HISTORY
%   22/05/2009 rename to isPointInCircle, add psb to specify tolerance

% extract computation tolerance
tol = 1e-14;
if ~isempty(varargin)
    tol = varargin{1};
end

d = sqrt(sum(power(point - circle(:,1:2), 2), 2));
b = d-circle(:,3)<=tol;

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function Dout = spm_eeg_merge(S)
% Concatenate epoched single trial files
% FORMAT Dout = spm_eeg_merge(S)
%
% S           - input structure (optional)
%  fields of S:
%   S.D       - character array containing filename of M/EEG mat-files
%               or cell array of D's
%   S.recode  - this field specifies how the condition labels will be
%               translated from the original files to the merged file.
%               Several options are possible:
%                 'same'        - leave the condition labels unchanged
%                 'addfilename' - add the original file name to condition
%                                 label
%                 old way specification - (for backward compatibility)
%                       a cell array where each cell contains a condition
%                       label. The ordering of these labels must be such
%                       that each row in the cell matrix specifies the
%                       conditionlabels for one of the selected files.
%                 specification via recoding rules - for this S.recode
%                       should be a structure array where each element
%                       specifies a rule using the following fields:
%                            file - can be a cell array of strings with
%                                   file names, a vector of file indices
%                                   or a string with regular expression
%                                   matching the files to which the rule
%                                   will apply.
%                            labelorg - can be a cell array of condition
%                                   labels or a string with regular
%                                   expression matching the condition
%                                   labels to which this rule will apply.
%                            labelnew - new label for the merged file. It
%                                   can contain special tokens #file# and
%                                   #labelorg# that will be replaced by
%                                   the original file name and original
%                                   condition label respectively.
%                       The rule will be applied one after the other so
%                       the last rule takes precedences. Trials not
%                       matched by any of the rules will keep their
%                       original labels.
%                       Example:
%                          S.recode(1).file     = '.*';
%                          S.recode(1).labelorg = '.*';
%                          S.recode(1).labelnew = '#labelorg# #file#';
%                       has the same effect as the 'addfilename' option.
%   S.prefix  - prefix for the output file (default - 'c')
%
%
% Dout        - MEEG object (also written to disk)
%__________________________________________________________________________
%
% This function can be used to merge M/EEG files to one file. This is
% useful whenever the data are distributed over multiple files, but one
% wants to use all information in one file. For example, when displaying
% data (SPM displays data from only one file at a time), or merging
% information that has been measured in multiple sessions.
%__________________________________________________________________________
% Copyright (C) 2008-2017 Wellcome Trust Centre for Neuroimaging

% Stefan Kiebel, Vladimir Litvak, Doris Eckstein, Rik Henson
% $Id: spm_eeg_merge.m 7125 2017-06-23 09:49:29Z guillaume $

SVNrev = '$Rev: 7125 $';

%-Startup
%--------------------------------------------------------------------------
spm('FnBanner', mfilename, SVNrev);
spm('FigName','M/EEG Merge'); spm('Pointer','Watch');

if ~isfield(S,'prefix'), S.prefix = 'c';    end
if ~isfield(S,'recode'), S.recode = 'same'; end

%-Load MEEG data
%--------------------------------------------------------------------------
D = S.D;

if ischar(D)
    F = cell(1,size(D,1));
    try
        for i = 1:size(D, 1)
            F{i} = spm_eeg_load(deblank(D(i, :)));
        end
        D = F;
    catch
        error('Trouble reading files.');
    end
end

Nfiles = length(D);

if Nfiles < 2
    %error('Need at least two files for merging.');
end

%-Check input and determine number of new number of trial types
%--------------------------------------------------------------------------
Ntrials = [];
megsens = [];
eegsens = [];
fid     = [];
isTF    =  strncmpi(D{1}.transformtype,'TF',2); % TF and TFphase

for i = 1:Nfiles
    if ~isequal(D{i}.transformtype, D{1}.transformtype)
        error(['The datasets do not contain the same kind of data.\n'...
               'There is a difference between files\n\t%s\nand\n\t%s.'], ...
               D{1}.fname, D{i}.fname);
    end

    if D{1}.nchannels ~= D{i}.nchannels
        error(['Data don''t have the same number of channels.\n' ...
               'There is a difference between files\n\t%s\nand\n\t%s.'], ...
               D{1}.fname, D{i}.fname);
    end

    if D{1}.nsamples ~= D{i}.nsamples
        error(['Data don''t have the same number of time points.\n' ...
               'There is a difference between files\n\t%s\nand\n\t%s.'], ...
               D{1}.fname, D{i}.fname);
    end

    if D{1}.fsample ~= D{i}.fsample
        error(['Data don''t have the same sampling rate.\n' ...
               'There is a difference between files\n\t%s\nand\n\t%s.'], ...
               D{1}.fname, D{i}.fname);
    end

    if isTF &&  ~isequal(D{1}.frequencies, D{i}.frequencies)
        error(['Data don''t have the same frequencies.\n' ...
               'There is a difference between files\n\t%s\nand\n\t%s.'], ...
               D{1}.fname, D{i}.fname);
    end

    if ~isempty(D{i}.sensors('MEG'))
        megsens = spm_cat_struct(megsens, D{i}.sensors('MEG'));
    end

    if ~isempty(D{i}.sensors('EEG'))
        eegsens = spm_cat_struct(eegsens, D{i}.sensors('EEG'));
    end

    if ~isempty(megsens) || ~isempty(eegsens)
        fid = spm_cat_struct(fid, D{i}.fiducials);
    end

    Ntrials = [Ntrials D{i}.ntrials];
end

%-Prepare some useful lists
%--------------------------------------------------------------------------
F        = {};
Find     = [];
clb      = {};
for i = 1:Nfiles
    F{i} = fname(D{i});
    clb  = [clb D{i}.conditions];
    Find = [Find i*ones(1, D{i}.ntrials)];
end
uclb     = unique(clb);

%-Specify condition labels recoding
%--------------------------------------------------------------------------
if ~isfield(S, 'recode')
    S.recode = spm_input('What to do with condition labels?', 1, 'm',...
        'Leave as they are|Add file name|Specify rules for recoding|Specify recoding the old way', strvcat('same', 'addfilename', 'rules', 'old'));
end

if isequal(S.recode, 'old')
    S.recode = {};
    for i = 1:Nfiles
        for j = 1:nconditions(D{i})
            S.recode{i}{j} = spm_input(sprintf('Labels: %s', spm_file(D{i}.fname, 'basename')),...
                '+1', 's', D{i}.condlist{j});
        end
    end
elseif isequal(S.recode, 'rules')
    S.recode = [];
    stop     = 0;
    ind      = 1;
    while ~stop
        spm_input(['Please define rule ' num2str(ind) ':'], 1, 'd');
        switch spm_input('To which files will this rule apply?', '+1', 'm',...
                'All the files|Specify indices|Select files|Wildcard expression (*,?)|Regular expression',...
                strvcat('all', 'indices', 'select', 'wildcard', 'regexp'))
            case 'all'
                S.recode(ind).file = '.*';
            case 'indices'
                S.recode(ind).file =  spm_input('Input file indices', '+1', 'n', num2str(ind), [1 Inf]);
            case 'select'
                [selection, ok]= listdlg('ListString', F, 'SelectionMode', 'multiple' ,'Name', 'Select files' , 'ListSize', [400 300]);
                if ok
                    S.recode(ind).file = F(selection);
                else
                    continue;
                end
            case 'wildcard'
                S.recode(ind).file = regexptranslate('wildcard' , spm_input('Input wildcard expresssion', '+1', 's',  '*'));
            case 'regexp'
                S.recode(ind).file = spm_input('Input regular expresssion', '+1', 's',  '.*');
        end

        switch spm_input('What conditions will be renamed?', '+1', 'm',...
                'All|Select|Specify by wildcard expression (*,?)|Specify by regular expression',...
                strvcat('all', 'select', 'wildcard', 'regexp'))
            case 'all'
                S.recode(ind).labelorg = '.*';
            case 'select'
                [selection, ok]= listdlg('ListString', uclb, 'SelectionMode', 'multiple' ,'Name', 'Select conditions' , 'ListSize', [400 300]);
                if ok
                    S.recode(ind).labelorg = uclb(selection);
                else
                    continue;
                end
            case 'wildcard'
                S.recode(ind).labelorg = regexptranslate('wildcard' , spm_input('Input wildcard expresssion', '+1', 's',  '*'));
            case 'regexp'
                S.recode(ind).labelorg = spm_input('Input regular expresssion', '+1', 's',  '.*');
        end

        S.recode(ind).labelnew = spm_input('Input the new name?', '+1', 's',  '');

        stop = spm_input('Define another rule?','+1','yes|stop', [0 1], 0);
        ind  = ind+1;
    end
end

%-Generate new meeg object with new filename
%--------------------------------------------------------------------------
Dout = D{1};
newfilename = spm_file(fnamedat(Dout), 'path',pwd, 'prefix',S.prefix);

if ~isTF
    Dout = clone(Dout, newfilename, [Dout.nchannels Dout.nsamples sum(Ntrials)]);
else
    Dout = clone(Dout, newfilename, [Dout.nchannels Dout.nfrequencies Dout.nsamples sum(Ntrials)]);
end

%-Perform condition labels recoding
%--------------------------------------------------------------------------
if isequal(S.recode, 'same')
    Dout = conditions(Dout, ':', clb);
elseif isequal(S.recode, 'addfilename')
    for i = 1:numel(clb)
        clb{i} = [clb{i} ' ' spm_file(F{Find(i)}, 'basename')];
    end
    Dout = conditions(Dout, ':', clb);
elseif iscell(S.recode)
    for i = 1:Nfiles
        ind = find(Find == i);

        for j = 1:D{i}.nconditions
            clb(ind(strmatch(D{i}.condlist{j}, clb(ind), 'exact'))) = S.recode{i}(j);
        end
    end
    Dout = conditions(Dout, ':', clb);
elseif isstruct(S.recode)
    clbnew = clb;

    for i = 1:numel(S.recode)
        if isnumeric(S.recode(i).file)
            ind = S.recode(i).file;
        elseif iscell(S.recode(i).file)
            ind = spm_match_str(F, S.recode(i).file);
        elseif ischar(S.recode(i).file)
            ind = find(~cellfun('isempty', regexp(F, S.recode(i).file)));
        else
            error('Invalid file specification in recoding rule.');
        end

        ind = find(ismember(Find, ind));

        if iscell(S.recode(i).labelorg)
            ind = ind(ismember(clb(ind), S.recode(i).labelorg));
        elseif ischar(S.recode(i).labelorg)
            ind = ind(~cellfun('isempty', regexp(clb(ind), S.recode(i).labelorg)));
        else
            error('Invalid original condition label specification in recoding rule.');
        end

        for j = 1:length(ind)
            labelnew = S.recode(i).labelnew;
            labelnew = strrep(labelnew, '#file#', spm_file(F{Find(ind(j))}, 'basename'));
            labelnew = strrep(labelnew, '#labelorg#', clb(ind(j)));

            clbnew{ind(j)} = labelnew;
        end
    end

    Dout = conditions(Dout, ':', clbnew);
end

%-Average sensor locations
%--------------------------------------------------------------------------
CmdLine = spm('CmdLine');
if CmdLine, h = []; end
if ~isempty(megsens)
    if ~CmdLine, spm_figure('GetWin','Graphics');clf; end
    if ~isempty(eegsens)
        if ~CmdLine, h = subplot(2, 1, 1); end
        aeegsens = ft_average_sens(eegsens, 'weights', Ntrials, 'feedback', h);
        Dout = sensors(Dout, 'EEG', aeegsens);

        if ~CmdLine, h = subplot(2, 1, 2); end
    else
        if ~CmdLine, h = axes; end
    end

    [amegsens,afid] = ft_average_sens(megsens, 'fiducials', fid, 'weights', Ntrials, 'feedback', h);
    Dout = sensors(Dout, 'MEG', amegsens);
    Dout = fiducials(Dout, afid);
elseif ~isempty(eegsens)
    if ~CmdLine, spm_figure('GetWin','Graphics');clf; end
    if ~CmdLine, h = axes; end
    [aeegsens,afid] = ft_average_sens(eegsens, 'fiducials', fid, 'weights', Ntrials, 'feedback', h);
    Dout = sensors(Dout, 'EEG', aeegsens);
    Dout = fiducials(Dout, afid);
end

%-Write files
%--------------------------------------------------------------------------
spm_progress_bar('Init', Nfiles, 'Files merged');

k = 0;

for i = 1:Nfiles

    ind = union(Dout.badchannels, D{i}.badchannels);
    if ~isempty(ind)
        Dout = badchannels(Dout, ind, 1);
    end

    % write trial-wise to avoid memory mapping error
    for j = 1:D{i}.ntrials
        k = k + 1;
        if ~isTF
            Dout(1:Dout.nchannels, 1:Dout.nsamples, k) =  D{i}(:,:,j);
        else
            Dout(1:Dout.nchannels, 1:Dout.nfrequencies, 1:Dout.nsamples, k) =  D{i}(:,:,:,j);
        end
        Dout = badtrials(Dout, k, badtrials(D{i}, j));
    end

    % Propagate some useful information from the original files to the
    % merged file
    Dout = repl(Dout, find(Find == i), D{i}.repl);
    Dout = trialonset(Dout, find(Find == i), D{i}.trialonset);
    Dout = trialtag(Dout, find(Find == i), D{i}.trialtag);
    Dout = events(Dout, find(Find == i), D{i}.events);

    spm_progress_bar('Set', i);

end

%-Save new M/EEG data
%--------------------------------------------------------------------------
Dout = Dout.history('spm_eeg_merge', S, 'reset');
save(Dout);

%-Cleanup
%--------------------------------------------------------------------------
spm_progress_bar('Clear');
fprintf('%-40s: %30s\n','Completed',spm('time'));                       %-#
spm('FigName','M/EEG merge: done'); spm('Pointer','Arrow');

author: spm


text:
classdef LaunchVehicleGroundObjectWayPt < matlab.mixin.SetGet
    %LaunchVehicleGroundObjectWayPt Summary of this class goes here
    %   Detailed explanation goes here

    properties
        elemSet(1,:) GeographicElementSet
        timeToNextWayPt(1,1) double = 0
    end

    methods
        function obj = LaunchVehicleGroundObjectWayPt(elemSet, timeToNextWayPt)
            obj.elemSet = elemSet;
            obj.timeToNextWayPt = timeToNextWayPt;
        end

        function str = getDisplayStr(obj)
            str = sprintf('Lat: %0.3f deg, Long: %0.3f deg, Alt: %0.3f km', ...
                          rad2deg(obj.getLatitude()), rad2deg(obj.getLongitude()), obj.getAltitude());
        end

        function timesToNextWayPt = getTimesToNextWaypt(obj)
            timesToNextWayPt = [obj.timeToNextWayPt];
        end

        function gElemSet = getElemSet(obj)
            gElemSet = obj.elemSet;
        end

        function lat = getLatitude(obj)
            lat = obj.elemSet.lat;
        end

        function long = getLongitude(obj)
            long = obj.elemSet.long;
        end

        function alt = getAltitude(obj)
            alt = obj.elemSet.alt;
        end

        function setLatitude(obj, newLat)
            obj.elemSet.lat = newLat;
        end

        function setLongitude(obj, newLong)
            obj.elemSet.long = newLong;
        end

        function setAltitude(obj, newAlt)
            obj.elemSet.alt = newAlt;
        end
    end

    methods(Static)
        function wayPt = getDefaultWayPt(initialTime, frameToUse)
            durToNextWayPt1 = 3600;
            elemSet = GeographicElementSet(initialTime + durToNextWayPt1, 0, 0, 0, 0, 0, 0, frameToUse);
            wayPt = LaunchVehicleGroundObjectWayPt(elemSet, durToNextWayPt1);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [khs, ghs, def] = HashinShtrikmanModulus(C,Ko,Go)
% Hashin Shtrikman moduli
%
% Description
% The equations are from Peselnick and Meister 1965 through Watt and
% Peselnick 1980 which are based on Hashin and Shtrikman 1963 JMB 8/2013
%
% Syntax:
%
%   [khs, ghs, def] = HashinShtrikmanModulus(C,Ko,Go)
%
% Input
%  C  - @stiffnessTensor
%  Ko - bulk modulus of the reference isotropic material
%  Go - shear modulus of the reference isotropic material
%
% Output
%  khs - Hashin Shtrikman bulk modulus
%  ghs - Hashin Shtrikman shear modulus
%  def - +1 positive definite, -1 for negative definite R
%

% ensure K0 and G0 have same size
if numel(Go) == 1, Go = repmat(Go,size(Ko)); end

% the isotropic compliances elements (eq. 10-11 Watt Peselnick 1980)
alpha = -3 ./ (3*Ko + 4*Go);
beta = -3 * (Ko + 2*Go) ./ ( 5*Go .* (3*Ko + 4*Go) );
gamma = (alpha - 3 * beta) ./ 9;

% define the isotropic tensor corresponding to Ko and Go
Co = 2 * Go .* stiffnessTensor.eye + (Ko - 2/3 * Go) .* dyad(tensor.eye,tensor.eye);

% take difference between anisotropic and isotropic matrices and take
% inverse to get the residual compliance matrix (eq 4 & 8 Watt Peselnick 1980)
R = C - Co;
H = inv(R);

% Subtract the isotropic compliances from the residual compliance matrix
% and invert back to moduli space.  The matrix B is the moduli matrix that
% when multiplied by the strain difference between the reference isotropic
% response and the isotropic response of the actual material gives the
% average stres <pij> (eq 15 - 16 Watt Peselnick 1980).
A = H - beta .* complianceTensor.eye - gamma .* dyad(tensor.eye,tensor.eye);
B = inv(A);

% now perform the averaging of the B matrix to get two numbers (B1 and B2)
% that are related to the two isotropic moduli - perturbations of the
% reference values (eq 21 and 22 Watt Peselnick 1980)
sB1 = EinsteinSum(B,[-1 -1 -2 -2]);
sB2 = EinsteinSum(B,[-1 -2 -1 -2]);

B1 = (2*sB1 - sB2) / 15;
B2 = (3*sB2 - sB1) / 30;

% The Hashin-Shtrikman moduli are then calculated as changes from the
% reference isotropic body. (eq 25 & 27 Watt and Peselnick 1980)
khs = Ko + ( 3 * B1 + 2 * B2 ) ./ ( 3 + alpha .* (3 * B1 + 2 * B2) );
ghs = Go + B2 ./ ( 1 + 2 * beta .* B2 );

% The moduli are valid in two limits - minimizing or maximizing the
% anisotropic difference elastic energy. Think of it as the maximum
% positive deviations from the reference state or the maximum negative
% deviations from the reference state. These are either in the positive
% definite or negative definite regime of the matrix R. Here is a test of
% the properties of R.  A +1 is returned if positive definite, a -1 is
% retunred if negative definite.  0 returned otherwise.

def = zeros(size(R));
D = eig(R);
def(all(D>0,1)) = 1;
def(all(D<0,1)) = -1;

% isf minmax == 0 this means ...
khs(def==0) = NaN;
ghs(def==0) = NaN;

author: mtex-toolbox

repo: mtex

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


text:
% This small example illustrates how to use the remote API
% synchronous mode. The synchronous mode needs to be
% pre-enabled on the server side. You would do this by
% starting the server (e.g. in a child script) with:
%
% simRemoteApi.start(19999,1300,false,true)
%
% But in this example we try to connect on port
% 19997 where there should be a continuous remote API
% server service already running and pre-enabled for
% synchronous mode.
%
% IMPORTANT: for each successful call to simxStart, there
% should be a corresponding call to simxFinish at the end!

function kuka_kart_robot()
    % init basic data. Robots and number of collision objects.
    coppelia = [];
    n_joints = 10;
    coppelia.n_collision_objects = 0;
    coppelia.robots{1}.n_joints = n_joints;

    coppelia = coppelia_init(coppelia);

    q = 0.1*ones(n_joints, 1);
    dq = 0.1*q;
    % Draw robot 1: R1, which possesses joints R1_q1, R1_q2...
    for i=1:40
       q = q + dq;
       drawrobotcoppelia(coppelia, 1, q)
       pause(0.01)
    end

    coppelia_end(coppelia);
end

repo: ARTE


text:
function [ekg_RSlinB_am_ELF_CtO, ekg_RSlinB_bw_ELF_CtO, ekg_RSlinB_fm_ELF_CtO] = estimate_rr_is(ekg_RSlinB_am_ELF, ekg_RSlinB_bw_ELF, ekg_RSlinB_fm_ELF, up)
% Estimation of respiration rate using multiple methods

subj = 1;
loaded_this_subj_respSigs = 0;
%% Make window timings if necessary
[wins] = identify_subj_wins_is(ekg_RSlinB_am_ELF, ekg_RSlinB_bw_ELF, ekg_RSlinB_fm_ELF, up);
%% Cycle through each resp signal
for respSig_no = 1:3    % 3 input signals
    rel_name = ['ekg_RSlinB_' up.al.options.FMe{respSig_no} '_ELF'];
    for option_no = 1 : length(up.al.options.estimate_rr)
        % Skip if this processing has been done previously
        if strcmp(up.al.options.estimate_rr{option_no}, 'GCE')
            %save_name = [ respSigs{respSig_no}(1:3) '_' up.al.options.estimate_rr{option_no} ];
            save_name = [ rel_name '_' up.al.options.estimate_rr{option_no} ];
        else
            %save_name = [ respSigs{respSig_no} '_' up.al.options.estimate_rr{option_no} ];
            save_name = [ rel_name '_' up.al.options.estimate_rr{option_no} ];
        end
%         savepath = [up.paths.data_save_folder, num2str(subj), up.paths.filenames.rrEsts, '.mat'];
%         exist_log = check_exists(savepath, save_name);
%         if exist_log
%             continue
%         end

%         % load data if it hasn't yet been loaded
%         if ~loaded_this_subj_respSigs
%             % Signals
%             load([up.paths.data_save_folder, num2str(subj), up.paths.filenames.respSigs]);
%             % Window timings
%             load([up.paths.data_save_folder, num2str(subj), up.paths.filenames.win_timings, '.mat']);
%             loaded_this_subj_respSigs = 1;
%         end
        % Identify the relevant respSig data
        eval(['rel_data = ' rel_name ';']);
        % add current subject and resp sig for any methods that do not use a resp sig.
        rel_data.subj = subj;
        rel_data.respSig = rel_name;
        %eval(['rel_data.respSig = ' rel_name ';']); %respSigs{respSig_no};
        %% Calculate RR from this resp sig using each option for estimating RR
        if (length(rel_data.t) == 1 && isnan(rel_data.t)) || sum(isnan(rel_data.v))==length(rel_data.v)
            temp_rr.t = mean([wins.t_start(:)' ; wins.t_end(:)']); temp_rr.t = temp_rr.t(:);
            temp_rr.v = nan(length(temp_rr.t),1);
            [temp_rr.f, temp_rr.p] = deal(cell(length(temp_rr.t),1));
        else
            temp_rr = feval(up.al.options.estimate_rr{option_no}, rel_data, wins, up);
        end
        % store this series of rrs:
        eval([rel_name '_' up.al.options.estimate_rr{option_no} ' = temp_rr;']);

        clear temp_rr
%         %% Save RRs to file
%         save_or_append_data
    end
end

end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function C = fdct_wrapping(x, is_real, finest, nbscales, nbangles_coarse)

% fdct_wrapping.m - Fast Discrete Curvelet Transform via wedge wrapping - Version 1.0
%
% Inputs
%   x           M-by-N matrix
%
% Optional Inputs
%   is_real     Type of the transform
%                   0: complex-valued curvelets
%                   1: real-valued curvelets
%               [default set to 0]
%   finest      Chooses one of two possibilities for the coefficients at the
%               finest level:
%                   1: curvelets
%                   2: wavelets
%               [default set to 2]
%   nbscales    number of scales including the coarsest wavelet level
%               [default set to ceil(log2(min(M,N)) - 3)]
%   nbangles_coarse
%               number of angles at the 2nd coarsest level, minimum 8,
%               must be a multiple of 4. [default set to 16]
%
% Outputs
%   C           Cell array of curvelet coefficients.
%               C{j}{l}(k1,k2) is the coefficient at
%                   - scale j: integer, from finest to coarsest scale,
%                   - angle l: integer, starts at the top-left corner and
%                   increases clockwise,
%                   - position k1,k2: both integers, size varies with j
%                   and l.
%               If is_real is 1, there are two types of curvelets,
%               'cosine' and 'sine'. For a given scale j, the 'cosine'
%               coefficients are stored in the first two quadrants (low
%               values of l), the 'sine' coefficients in the last two
%               quadrants (high values of l).
%
% See also ifdct_wrapping.m, fdct_wrapping_param.m
%
% By Laurent Demanet, 2004

X = fftshift(fft2(ifftshift(x)))/sqrt(prod(size(x)));
[N1,N2] = size(X);
if nargin < 2, is_real = 0; end;
if nargin < 3, finest = 2; end;
if nargin < 4, nbscales = ceil(log2(min(N1,N2)) - 3); end;
if nargin < 5, nbangles_coarse = 16; end;

% Initialization: data structure
nbangles = [1, nbangles_coarse .* 2.^(ceil((nbscales-(nbscales:-1:2))/2))];
if finest == 2, nbangles(nbscales) = 1; end;
C = cell(1,nbscales);
for j = 1:nbscales
    C{j} = cell(1,nbangles(j));
end;

% Loop: pyramidal scale decomposition
M1 = N1/3;
M2 = N2/3;
if finest == 1,

    % Initialization: smooth periodic extension of high frequencies
    bigN1 = 2*floor(2*M1)+1;
    bigN2 = 2*floor(2*M2)+1;
    equiv_index_1 = 1+mod(floor(N1/2)-floor(2*M1)+(1:bigN1)-1,N1);
    equiv_index_2 = 1+mod(floor(N2/2)-floor(2*M2)+(1:bigN2)-1,N2);
    X = X(equiv_index_1,equiv_index_2);
        % Invariant: equiv_index_1(floor(2*M1)+1) == (N1 + 2 - mod(N1,2))/2
        % is the center in frequency. Same for M2, N2.
    window_length_1 = floor(2*M1) - floor(M1) - 1 - (mod(N1,3)==0);
    window_length_2 = floor(2*M2) - floor(M2) - 1 - (mod(N2,3)==0);
        % Invariant: floor(M1) + floor(2*M1) == N1 - (mod(M1,3)~=0)
        % Same for M2, N2.
    coord_1 = 0:(1/window_length_1):1;
    coord_2 = 0:(1/window_length_2):1;
    [wl_1,wr_1] = fdct_wrapping_window(coord_1);
    [wl_2,wr_2] = fdct_wrapping_window(coord_2);
    lowpass_1 = [wl_1, ones(1,2*floor(M1)+1), wr_1];
    if mod(N1,3)==0, lowpass_1 = [0, lowpass_1, 0]; end;
    lowpass_2 = [wl_2, ones(1,2*floor(M2)+1), wr_2];
    if mod(N2,3)==0, lowpass_2 = [0, lowpass_2, 0]; end;
    lowpass = lowpass_1'*lowpass_2;
    Xlow = X .* lowpass;

    scales = nbscales:-1:2;

else

    M1 = M1/2;
    M2 = M2/2;
    window_length_1 = floor(2*M1) - floor(M1) - 1;
    window_length_2 = floor(2*M2) - floor(M2) - 1;
    coord_1 = 0:(1/window_length_1):1;
    coord_2 = 0:(1/window_length_2):1;
    [wl_1,wr_1] = fdct_wrapping_window(coord_1);
    [wl_2,wr_2] = fdct_wrapping_window(coord_2);
    lowpass_1 = [wl_1, ones(1,2*floor(M1)+1), wr_1];
    lowpass_2 = [wl_2, ones(1,2*floor(M2)+1), wr_2];
    lowpass = lowpass_1'*lowpass_2;
    hipass = sqrt(1 - lowpass.^2);
    Xlow_index_1 = ((-floor(2*M1)):floor(2*M1)) + ceil((N1+1)/2);
    Xlow_index_2 = ((-floor(2*M2)):floor(2*M2)) + ceil((N2+1)/2);
    Xlow = X(Xlow_index_1, Xlow_index_2) .* lowpass;
    Xhi = X;
    Xhi(Xlow_index_1, Xlow_index_2) = Xhi(Xlow_index_1, Xlow_index_2) .* hipass;
    C{nbscales}{1} = fftshift(ifft2(ifftshift(Xhi)))*sqrt(prod(size(Xhi)));
    if is_real, C{nbscales}{1} = real(C{nbscales}{1}); end;

    scales = (nbscales-1):-1:2;

end;
for j = scales,

    M1 = M1/2;
    M2 = M2/2;
    window_length_1 = floor(2*M1) - floor(M1) - 1;
    window_length_2 = floor(2*M2) - floor(M2) - 1;
    coord_1 = 0:(1/window_length_1):1;
    coord_2 = 0:(1/window_length_2):1;
    [wl_1,wr_1] = fdct_wrapping_window(coord_1);
    [wl_2,wr_2] = fdct_wrapping_window(coord_2);
    lowpass_1 = [wl_1, ones(1,2*floor(M1)+1), wr_1];
    lowpass_2 = [wl_2, ones(1,2*floor(M2)+1), wr_2];
    lowpass = lowpass_1'*lowpass_2;
    hipass = sqrt(1 - lowpass.^2);
    Xhi = Xlow;                 % size is 2*floor(4*M1)+1 - by - 2*floor(4*M2)+1
    Xlow_index_1 = ((-floor(2*M1)):floor(2*M1)) + floor(4*M1) + 1;
    Xlow_index_2 = ((-floor(2*M2)):floor(2*M2)) + floor(4*M2) + 1;
    Xlow = Xlow(Xlow_index_1, Xlow_index_2);
    Xhi(Xlow_index_1, Xlow_index_2) = Xlow .* hipass;
    Xlow = Xlow .* lowpass;     % size is 2*floor(2*M1)+1 - by - 2*floor(2*M2)+1

    % Loop: angular decomposition
    l = 0;
    nbquadrants = 2 + 2*(~is_real);
    nbangles_perquad = nbangles(j)/4;
    for quadrant = 1:nbquadrants
        M_horiz = M2 * (mod(quadrant,2)==1) + M1 * (mod(quadrant,2)==0);
        M_vert = M1 * (mod(quadrant,2)==1) + M2 * (mod(quadrant,2)==0);
        if mod(nbangles_perquad,2),
            wedge_ticks_left = round((0:(1/(2*nbangles_perquad)):.5)*2*floor(4*M_horiz) + 1);
            wedge_ticks_right = 2*floor(4*M_horiz) + 2 - wedge_ticks_left;
            wedge_ticks = [wedge_ticks_left, wedge_ticks_right(end:-1:1)];
        else
            wedge_ticks_left = round((0:(1/(2*nbangles_perquad)):.5)*2*floor(4*M_horiz) + 1);
            wedge_ticks_right = 2*floor(4*M_horiz) + 2 - wedge_ticks_left;
            wedge_ticks = [wedge_ticks_left, wedge_ticks_right((end-1):-1:1)];
        end;
        wedge_endpoints = wedge_ticks(2:2:(end-1));         % integers
        wedge_midpoints = (wedge_endpoints(1:(end-1)) + wedge_endpoints(2:end))/2;
                % integers or half-integers

        % Left corner wedge
        l = l+1;
        first_wedge_endpoint_vert = round(2*floor(4*M_vert)/(2*nbangles_perquad) + 1);
        length_corner_wedge = floor(4*M_vert) - floor(M_vert) + ceil(first_wedge_endpoint_vert/4);
        Y_corner = 1:length_corner_wedge;
        [XX,YY] = meshgrid(1:(2*floor(4*M_horiz)+1),Y_corner);
        width_wedge = wedge_endpoints(2) + wedge_endpoints(1) - 1;
        slope_wedge = (floor(4*M_horiz) + 1 - wedge_endpoints(1))/floor(4*M_vert);
        left_line = round(2 - wedge_endpoints(1) + slope_wedge*(Y_corner - 1));
                                                            % integers
        [wrapped_data, wrapped_XX, wrapped_YY] = deal(zeros(length_corner_wedge,width_wedge));
        first_row = floor(4*M_vert)+2-ceil((length_corner_wedge+1)/2)+...
            mod(length_corner_wedge+1,2)*(quadrant-2 == mod(quadrant-2,2));
        first_col = floor(4*M_horiz)+2-ceil((width_wedge+1)/2)+...
            mod(width_wedge+1,2)*(quadrant-3 == mod(quadrant-3,2));
                % Coordinates of the top-left corner of the wedge wrapped
                % around the origin. Some subtleties when the wedge is
                % even-sized because of the forthcoming 90 degrees rotation
        for row = Y_corner
            cols = left_line(row) + mod((0:(width_wedge-1))-(left_line(row)-first_col),width_wedge);
            admissible_cols = round(1/2*(cols+1+abs(cols-1)));
            new_row = 1 + mod(row - first_row, length_corner_wedge);
            wrapped_data(new_row,:) = Xhi(row,admissible_cols) .* (cols > 0);
            wrapped_XX(new_row,:) = XX(row,admissible_cols);
            wrapped_YY(new_row,:) = YY(row,admissible_cols);
        end;
        slope_wedge_right = (floor(4*M_horiz)+1 - wedge_midpoints(1))/floor(4*M_vert);
        mid_line_right = wedge_midpoints(1) + slope_wedge_right*(wrapped_YY - 1);
                % not integers in general
        coord_right = 1/2 + floor(4*M_vert)/(wedge_endpoints(2) - wedge_endpoints(1)) * ...
            (wrapped_XX - mid_line_right)./(floor(4*M_vert)+1 - wrapped_YY);
        C2 = 1/(1/(2*(floor(4*M_horiz))/(wedge_endpoints(1) - 1) - 1) + 1/(2*(floor(4*M_vert))/(first_wedge_endpoint_vert - 1) - 1));
        C1 = C2 / (2*(floor(4*M_vert))/(first_wedge_endpoint_vert - 1) - 1);
        wrapped_XX((wrapped_XX - 1)/floor(4*M_horiz) + (wrapped_YY-1)/floor(4*M_vert) == 2) = ...
            wrapped_XX((wrapped_XX - 1)/floor(4*M_horiz) + (wrapped_YY-1)/floor(4*M_vert) == 2) + 1;
        coord_corner = C1 + C2 * ((wrapped_XX - 1)/(floor(4*M_horiz)) - (wrapped_YY - 1)/(floor(4*M_vert))) ./ ...
            (2-((wrapped_XX - 1)/(floor(4*M_horiz)) + (wrapped_YY - 1)/(floor(4*M_vert))));
        wl_left = fdct_wrapping_window(coord_corner);
        [wl_right,wr_right] = fdct_wrapping_window(coord_right);
        wrapped_data = wrapped_data .* (wl_left .* wr_right);

        switch is_real
            case 0
                wrapped_data = rot90(wrapped_data,-(quadrant-1));
                C{j}{l} = fftshift(ifft2(ifftshift(wrapped_data)))*sqrt(prod(size(wrapped_data)));
            case 1
                wrapped_data = rot90(wrapped_data,-(quadrant-1));
                x = fftshift(ifft2(ifftshift(wrapped_data)))*sqrt(prod(size(wrapped_data)));
                C{j}{l} = sqrt(2)*real(x);
                C{j}{l+nbangles(j)/2} = sqrt(2)*imag(x);
        end;

        % Regular wedges
        length_wedge = floor(4*M_vert) - floor(M_vert);
        Y = 1:length_wedge;
        first_row = floor(4*M_vert)+2-ceil((length_wedge+1)/2)+...
            mod(length_wedge+1,2)*(quadrant-2 == mod(quadrant-2,2));
        for subl = 2:(nbangles_perquad-1);
            l = l+1;
            width_wedge = wedge_endpoints(subl+1) - wedge_endpoints(subl-1) + 1;
            slope_wedge = ((floor(4*M_horiz)+1) - wedge_endpoints(subl))/floor(4*M_vert);
            left_line = round(wedge_endpoints(subl-1) + slope_wedge*(Y - 1));
            [wrapped_data, wrapped_XX, wrapped_YY] = deal(zeros(length_wedge,width_wedge));
            first_col = floor(4*M_horiz)+2-ceil((width_wedge+1)/2)+...
                mod(width_wedge+1,2)*(quadrant-3 == mod(quadrant-3,2));
            for row = Y
                cols = left_line(row) + mod((0:(width_wedge-1))-(left_line(row)-first_col),width_wedge);
                new_row = 1 + mod(row - first_row, length_wedge);
                wrapped_data(new_row,:) = Xhi(row,cols);
                wrapped_XX(new_row,:) = XX(row,cols);
                wrapped_YY(new_row,:) = YY(row,cols);
            end;
            slope_wedge_left = ((floor(4*M_horiz)+1) - wedge_midpoints(subl-1))/floor(4*M_vert);
            mid_line_left = wedge_midpoints(subl-1) + slope_wedge_left*(wrapped_YY - 1);
            coord_left = 1/2 + floor(4*M_vert)/(wedge_endpoints(subl) - wedge_endpoints(subl-1)) * ...
                (wrapped_XX - mid_line_left)./(floor(4*M_vert)+1 - wrapped_YY);
            slope_wedge_right = ((floor(4*M_horiz)+1) - wedge_midpoints(subl))/floor(4*M_vert);
            mid_line_right = wedge_midpoints(subl) + slope_wedge_right*(wrapped_YY - 1);
            coord_right = 1/2 + floor(4*M_vert)/(wedge_endpoints(subl+1) - wedge_endpoints(subl)) * ...
                (wrapped_XX - mid_line_right)./(floor(4*M_vert)+1 - wrapped_YY);
            wl_left = fdct_wrapping_window(coord_left);
            [wl_right,wr_right] = fdct_wrapping_window(coord_right);
            wrapped_data = wrapped_data .* (wl_left .* wr_right);
            switch is_real
                case 0
                    wrapped_data = rot90(wrapped_data,-(quadrant-1));
                    C{j}{l} = fftshift(ifft2(ifftshift(wrapped_data)))*sqrt(prod(size(wrapped_data)));
                case 1
                    wrapped_data = rot90(wrapped_data,-(quadrant-1));
                    x = fftshift(ifft2(ifftshift(wrapped_data)))*sqrt(prod(size(wrapped_data)));
                    C{j}{l} = sqrt(2)*real(x);
                    C{j}{l+nbangles(j)/2} = sqrt(2)*imag(x);
            end;
        end;

        % Right corner wedge
        l = l+1;
        width_wedge = 4*floor(4*M_horiz) + 3 - wedge_endpoints(end) - wedge_endpoints(end-1);
        slope_wedge = ((floor(4*M_horiz)+1) - wedge_endpoints(end))/floor(4*M_vert);
        left_line = round(wedge_endpoints(end-1) + slope_wedge*(Y_corner - 1));
        [wrapped_data, wrapped_XX, wrapped_YY] = deal(zeros(length_corner_wedge,width_wedge));
        first_row = floor(4*M_vert)+2-ceil((length_corner_wedge+1)/2)+...
            mod(length_corner_wedge+1,2)*(quadrant-2 == mod(quadrant-2,2));
        first_col = floor(4*M_horiz)+2-ceil((width_wedge+1)/2)+...
            mod(width_wedge+1,2)*(quadrant-3 == mod(quadrant-3,2));
        for row = Y_corner
            cols = left_line(row) + mod((0:(width_wedge-1))-(left_line(row)-first_col),width_wedge);
            admissible_cols = round(1/2*(cols+2*floor(4*M_horiz)+1-abs(cols-(2*floor(4*M_horiz)+1))));
            new_row = 1 + mod(row - first_row, length_corner_wedge);
            wrapped_data(new_row,:) = Xhi(row,admissible_cols) .* (cols <= (2*floor(4*M_horiz)+1));
            wrapped_XX(new_row,:) = XX(row,admissible_cols);
            wrapped_YY(new_row,:) = YY(row,admissible_cols);
        end;
        slope_wedge_left = ((floor(4*M_horiz)+1) - wedge_midpoints(end))/floor(4*M_vert);
        mid_line_left = wedge_midpoints(end) + slope_wedge_left*(wrapped_YY - 1);
        coord_left = 1/2 + floor(4*M_vert)/(wedge_endpoints(end) - wedge_endpoints(end-1)) * ...
            (wrapped_XX - mid_line_left)./(floor(4*M_vert) + 1 - wrapped_YY);
        C2 = -1/(2*(floor(4*M_horiz))/(wedge_endpoints(end) - 1) - 1 + 1/(2*(floor(4*M_vert))/(first_wedge_endpoint_vert - 1) - 1));
        C1 = -C2 * (2*(floor(4*M_horiz))/(wedge_endpoints(end) - 1) - 1);
        wrapped_XX((wrapped_XX - 1)/floor(4*M_horiz) == (wrapped_YY - 1)/floor(4*M_vert)) = ...
            wrapped_XX((wrapped_XX - 1)/floor(4*M_horiz) == (wrapped_YY - 1)/floor(4*M_vert)) - 1;
        coord_corner = C1 + C2 * (2-((wrapped_XX - 1)/(floor(4*M_horiz)) + (wrapped_YY - 1)/(floor(4*M_vert)))) ./ ...
            ((wrapped_XX - 1)/(floor(4*M_horiz)) - (wrapped_YY - 1)/(floor(4*M_vert)));
        wl_left = fdct_wrapping_window(coord_left);
        [wl_right,wr_right] = fdct_wrapping_window(coord_corner);

        wrapped_data = wrapped_data .* (wl_left .* wr_right);
        switch is_real
            case 0
                wrapped_data = rot90(wrapped_data,-(quadrant-1));
                C{j}{l} = fftshift(ifft2(ifftshift(wrapped_data)))*sqrt(prod(size(wrapped_data)));
            case 1
                wrapped_data = rot90(wrapped_data,-(quadrant-1));
                x = fftshift(ifft2(ifftshift(wrapped_data)))*sqrt(prod(size(wrapped_data)));
                C{j}{l} = sqrt(2)*real(x);
                C{j}{l+nbangles(j)/2} = sqrt(2)*imag(x);
        end;

        if quadrant < nbquadrants, Xhi = rot90(Xhi); end;
    end;
end;

% Coarsest wavelet level
C{1}{1} = fftshift(ifft2(ifftshift(Xlow)))*sqrt(prod(size(Xlow)));
if is_real == 1,
    C{1}{1} = real(C{1}{1});
end;

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   PARAMETERS Returns a data structure containing the parameters of the
%   KUKA KR 30 3.
%
%   Author:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function robot = parameters()

robot.name= 'KR30_3';

robot.DH.theta= '[  q(1)  q(2)-pi/2     q(3)    q(4)    q(5)   q(6)]';
robot.DH.d='[       0.815     0           0       0.820   0    0.17]';
robot.DH.a='[       0.350   0.850         0.145    0       0      0]';
robot.DH.alpha= '[  -pi/2   0           -pi/2    pi/2    -pi/2   0]';
robot.J=[];

robot.inversekinematic_fn = 'inversekinematic_kuka_kr303(robot, T)';

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

%rotational: 0, translational: 1
robot.kind=['R' 'R' 'R' 'R' 'R' 'R'];

%minimum and maximum rotation angle in rad
robot.maxangle =[deg2rad(-185) deg2rad(185); %Axis 1, minimum, maximum
                deg2rad(-35) deg2rad(135); %Axis 2, minimum, maximum
                deg2rad(-120) deg2rad(158); %Axis 3
                deg2rad(-350) deg2rad(350); %Axis 4: Unlimited (400� default)
                deg2rad(-119) deg2rad(119); %Axis 5
                deg2rad(-350) deg2rad(350)]; %Axis 6: Unlimited (800� default)

%maximum absolute speed of each joint rad/s or m/s
robot.velmax = [deg2rad(140); %Axis 1, rad/s
                deg2rad(126); %Axis 2, rad/s
                deg2rad(140); %Axis 3, rad/s
                deg2rad(260); %Axis 4, rad/s
                deg2rad(245); %Axis 5, rad/s
                deg2rad(322)];%Axis 6, rad/s
% end effectors maximum velocity
robot.linear_velmax = 1.0; %m/s, not specified

robot.accelmax=robot.velmax/0.1; % 0.1 is here an acceleration time
%base reference system
robot.T0 = eye(4);

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

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

%DYNAMICS
robot.has_dynamics=1;

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

%link masses (kg)
robot.dynamics.masses=[266 133 99.75 66.5 66.5 33.25];

%COM of each link with respect to own reference system
robot.dynamics.r_com=[0       0          0; %(rx, ry, rz) link 0
                       -0.35	0.252	 0; %(rx, ry, rz) link 1
                       -0.675	   0	 -0.225;  %(rx, ry, rz) link 2
                        0       0       0.150;%(rx, ry, rz) link 4
                        0       0           0;%(rx, ry, rz) link 5
                        0       0         0.1075];%(rx, ry, rz) link 6

%Inertia matrices of each link with respect to its D-H reference system.
% Ixx   Iyy	Izz	Ixy	Iyz	Ixz, for each row
% % % % % robot.dynamics.Inertia=[0      0.35	0   	0	0	0;
% % % % %                         9.89     9.89	16	    0	0	0;
% % % % %                         38.78	38.78	2.56	0	0	0;
% % % % %                         34.2	34.2	0.9	    0	0	0;
% % % % %                         0.214	0.214	0.162	0	0	0;
% % % % %                         0.0356	0.0356	0.07	0	0	0];

robot.dynamics.Inertia=[0      0.35	0   	0	0	0;
    .13     .524	.539	0	0	0;
    .066	.086	.0125	0	0	0;
    1.8e-3	1.3e-3	1.8e-3	0	0	0;
    .3e-3	.4e-3	.3e-3	0	0	0;
    .15e-3	.15e-3	.04e-3	0	0	0];

robot.motors=load_motors([2 4 3 1 3 1]);
%Speed reductor at each joint
robot.motors.G=[600 600 600 100 100 50];


text:
function [ x, seed ] = truncated_normal_ab_sample ( mu, sigma, a, b, seed )

%*****************************************************************************80
%
%% TRUNCATED_NORMAL_AB_SAMPLE samples the truncated Normal PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 August 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real MU, SIGMA, the mean and standard deviation of the
%    parent Normal distribution.
%
%    Input, real A, B, the lower and upper truncation limits.
%
%    Input/output, integer SEED, a seed for the random number
%    generator.
%
%    Output, real X, a sample of the PDF.
%
  alpha = ( a - mu ) / sigma;
  beta = ( b - mu ) / sigma;

  alpha_cdf = normal_01_cdf ( alpha );
  beta_cdf = normal_01_cdf ( beta );

  [ u, seed ] = r8_uniform_01 ( seed );
  xi_cdf = alpha_cdf + u * ( beta_cdf - alpha_cdf );
  xi = normal_01_cdf_inv ( xi_cdf );

  x = mu + sigma * xi;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%TROTX SE(3) rotation about X axis
%
% T = TROTX(THETA) is a homogeneous transformation (4x4) representing a rotation
% of THETA radians about the x-axis.
%
% T = TROTX(THETA, 'deg') as above but THETA is in degrees.
%
% Notes::
% - Translational component is zero.
%
% See also rotx, troty, trotz, trot2, SE3.Rx.

%## 3d rotation

% Copyright (C) 1993-2019 Peter I. Corke
%
% This file is part of The Spatial Math Toolbox for MATLAB (SMTB).
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, including without limitation the rights
% to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
% of the Software, and to permit persons to whom the Software is furnished to do
% so, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in all
% copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
% FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
% COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
% IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
% CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
%
% https://github.com/petercorke/spatial-math

function T = trotx(t, varargin)
	T = [rotx(t, varargin{:}) [0 0 0]'; 0 0 0 1];

author: petercorke

repo: spatialmath-matlab


text:
function varargout = process_aec1n( varargin )
% PROCESS_PLV1N: Compute amplitude envelope correlation between all the pairs of signals, in one file.

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2012-2014
%          Peter Donhauser, 2017

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Amplitude Envelope Correlation NxN';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = 'Connectivity';
    sProcess.Index       = 682;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/Connectivity';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'data',     'results',  'matrix'};
    sProcess.OutputTypes = {'timefreq', 'timefreq', 'timefreq'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;

    % === CONNECT INPUT
    sProcess = process_corr1n('DefineConnectOptions', sProcess, 1);
    % === FREQ BANDS
    sProcess.options.freqbands.Comment = 'Frequency bands for the Hilbert transform:';
    sProcess.options.freqbands.Type    = 'groupbands';
    sProcess.options.freqbands.Value   = bst_get('DefaultFreqBands');
    % === KEEP TIME
    sProcess.options.isorth.Comment = 'Orthogonalize signal pairs before envelope computation';
    sProcess.options.isorth.Type    = 'checkbox';
    sProcess.options.isorth.Value   = 0;
    % === OUTPUT MODE
    sProcess.options.outputmode.Comment = {'Save individual results (one file per input file)', 'Concatenate input files before processing (one file)', 'Save average connectivity matrix (one file)'};
    sProcess.options.outputmode.Type    = 'radio';
    sProcess.options.outputmode.Value   = 1;
    sProcess.options.outputmode.Group   = 'output';
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    Comment = sProcess.Comment;
end

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputA) %#ok<DEFNU>
    % Input options
    OPTIONS = process_corr1n('GetConnectOptions', sProcess, sInputA);
    if isempty(OPTIONS)
        OutputFiles = {};
        return
    end
    OPTIONS.Method = 'aec';
    % Filtering bands options
    OPTIONS.Freqs = sProcess.options.freqbands.Value;
    OPTIONS.isOrth = sProcess.options.isorth.Value;

    % Compute metric
    OutputFiles = bst_connectivity({sInputA.FileName}, [], OPTIONS);
end


text:
% -------------------------------------------------------------------------
% Downloaded 10.08.2018 from:
% http://www.columbia.edu/~bsm2105/type2sdt/archive/index.html
% Author: Brian Maniscalco
% Contact: brian@psych.columbia.edu

% 10.08.2018: OKH added calculation for type 1 c to the results

% NOTE FROM WEBPAGE:
% If you use this analysis file, please reference the Consciousness &
% Cognition paper:
%       Maniscalco, B., & Lau, H. (2012). A signal detection
%       theoretic approach for estimating metacognitive sensitivity from
%       confidence ratings. Consciousness and Cognition, 21(1), 422–430.
%       doi:10.1016/j.concog.2011.09.021)
% And the website:
%       http://www.columbia.edu/~bsm2105/type2sdt/archive/index.html
% -------------------------------------------------------------------------

function fit = tapas_fit_meta_d_MLE(nR_S1, nR_S2, s, fncdf, fninv)

% fit = fit_meta_d_MLE(nR_S1, nR_S2, s, fncdf, fninv)
%
% Given data from an experiment where an observer discriminates between two
% stimulus alternatives on every trial and provides confidence ratings,
% provides a type 2 SDT analysis of the data.
%
% INPUTS
%
% * nR_S1, nR_S2
% these are vectors containing the total number of responses in
% each response category, conditional on presentation of S1 and S2.
%
% e.g. if nR_S1 = [100 50 20 10 5 1], then when stimulus S1 was
% presented, the subject had the following response counts:
% responded S1, rating=3 : 100 times
% responded S1, rating=2 : 50 times
% responded S1, rating=1 : 20 times
% responded S2, rating=1 : 10 times
% responded S2, rating=2 : 5 times
% responded S2, rating=3 : 1 time
%
% The ordering of response / rating counts for S2 should be the same as it
% is for S1. e.g. if nR_S2 = [3 7 8 12 27 89], then when stimulus S2 was
% presented, the subject had the following response counts:
% responded S1, rating=3 : 3 times
% responded S1, rating=2 : 7 times
% responded S1, rating=1 : 8 times
% responded S2, rating=1 : 12 times
% responded S2, rating=2 : 27 times
% responded S2, rating=3 : 89 times
%
% N.B. if nR_S1 or nR_S2 contain zeros, this may interfere with estimation of
% meta-d'.
%
% Some options for dealing with response cell counts containing zeros are:
%
% (1) Add a small adjustment factor, e.g. adj_f = 1/(length(nR_S1), to each
% input vector:
%
% adj_f = 1/length(nR_S1);
% nR_S1_adj = nR_S1 + adj_f;
% nR_S2_adj = nR_S2 + adj_f;
%
% This is a generalization of the correction for similar estimation issues of
% type 1 d' as recommended in
%
% Hautus, M. J. (1995). Corrections for extreme proportions and their biasing
%     effects on estimated values of d'. Behavior Research Methods, Instruments,
%     & Computers, 27, 46-51.
%
% When using this correction method, it is recommended to add the adjustment
% factor to ALL data for all subjects, even for those subjects whose data is
% not in need of such correction, in order to avoid biases in the analysis
% (cf Snodgrass & Corwin, 1988).
%
% (2) Collapse across rating categories.
%
% e.g. if your data set has 4 possible confidence ratings such that length(nR_S1)==8,
% defining new input vectors
%
% nR_S1_new = [sum(nR_S1(1:2)), sum(nR_S1(3:4)), sum(nR_S1(5:6)), sum(nR_S1(7:8))];
% nR_S2_new = [sum(nR_S2(1:2)), sum(nR_S2(3:4)), sum(nR_S2(5:6)), sum(nR_S2(7:8))];
%
% might be sufficient to eliminate zeros from the input without using an adjustment.
%
% * s
% this is the ratio of standard deviations for type 1 distributions, i.e.
%
% s = sd(S1) / sd(S2)
%
% if not specified, s is set to a default value of 1.
% For most purposes, we recommend setting s = 1.
% See http://www.columbia.edu/~bsm2105/type2sdt for further discussion.
%
% * fncdf
% a function handle for the CDF of the type 1 distribution.
% if not specified, fncdf defaults to @normcdf (i.e. CDF for normal
% distribution)
%
% * fninv
% a function handle for the inverse CDF of the type 1 distribution.
% if not specified, fninv defaults to @norminv
%
% OUTPUT
%
% Output is packaged in the struct "fit."
% In the following, let S1 and S2 represent the distributions of evidence
% generated by stimulus classes S1 and S2.
% Then the fields of "fit" are as follows:
%
% fit.d1        = mean(S2) - mean(S1), in room-mean-square(sd(S1),sd(S2)) units
% fit.c1        = type 1 criterion
% fit.s         = sd(S1) / sd(S2)
% fit.meta_d    = meta-d' in RMS units
% fit.M_diff    = meta_da - da
% fit.M_ratio   = meta_da / da
% fit.meta_c1   = type 1 criterion for meta-d' fit, RMS units
% fit.t2ca_rS1  = type 2 criteria of "S1" responses for meta-d' fit, RMS units
% fit.t2ca_rS2  = type 2 criteria of "S2" responses for meta-d' fit, RMS units
%
% fit.S1units   = contains same parameters in sd(S1) units.
%                 these may be of use since the data-fitting is conducted
%                 using parameters specified in sd(S1) units.
%
% fit.logL          = log likelihood of the data fit
%
% fit.est_HR2_rS1  = estimated (from meta-d' fit) type 2 hit rates for S1 responses
% fit.obs_HR2_rS1  = actual type 2 hit rates for S1 responses
% fit.est_FAR2_rS1 = estimated type 2 false alarm rates for S1 responses
% fit.obs_FAR2_rS1 = actual type 2 false alarm rates for S1 responses
%
% fit.est_HR2_rS2  = estimated type 2 hit rates for S2 responses
% fit.obs_HR2_rS2  = actual type 2 hit rates for S2 responses
% fit.est_FAR2_rS2 = estimated type 2 false alarm rates for S2 responses
% fit.obs_FAR2_rS2 = actual type 2 false alarm rates for S2 responses
%
% If there are N ratings, then there will be N-1 type 2 hit rates and false
% alarm rates.

% 2015/07/23 - fixed bug for output fit.meta_ca and fit.S1units.meta_c1.
%            - added comments to help section as well as a warning output
%              for nR_S1 or nR_S2 inputs containing zeros
% 2014/10/14 - updated discussion of "s" input in the help section above.
% 2010/09/07 - created

%% parse inputs

% check inputs
if ~mod(length(nR_S1),2)==0, error('input arrays must have an even number of elements'); end
if length(nR_S1)~=length(nR_S2), error('input arrays must have the same number of elements'); end

if any(nR_S1 == 0) || any(nR_S2 == 0)
    disp(' ')
    disp('WARNING!!')
    disp('---------')
    disp('Your inputs')
    disp(' ')
    disp(['nR_S1 = [' num2str(nR_S1) ']'])
    disp(['nR_S2 = [' num2str(nR_S2) ']'])
    disp(' ')
    disp('contain zeros! This may interfere with proper estimation of meta-d''.')
    disp('See ''help fit_meta_d_MLE'' for more information.')
    disp(' ')
    disp(' ')
end

% assign input default values
if ~exist('s','var') || isempty(s)
    s = 1;
end

if ~exist('fncdf','var') || isempty(fncdf)
    fncdf = @normcdf;
end

if ~exist('fninv','var') || isempty(fninv)
    fninv = @norminv;
end

nRatings = length(nR_S1) / 2;
nCriteria = 2*nRatings - 1;

%% set up constraints for MLE estimation

% parameters
% meta-d' - 1
% t2c     - nCriteria-1

A = [];
b = [];

% constrain type 2 criteria values,
% such that t2c(i) is always <= t2c(i+1)
% want t2c(i)   <= t2c(i+1)
% -->  t2c(i+1) >= c(i) + 1e-5 (i.e. very small deviation from equality)
% -->  t2c(i) - t2c(i+1) <= -1e-5
for i = 2 : nCriteria-1

    A(end+1,[i i+1]) = [1 -1];
    b(end+1) = -1e-5;

end

% lower bounds on parameters
LB = [];
LB = [LB -10];                         % meta-d'
LB = [LB -20*ones(1,(nCriteria-1)/2)]; % criteria lower than t1c
LB = [LB zeros(1,(nCriteria-1)/2)];    % criteria higher than t1c

% upper bounds on parameters
UB = [];
UB = [UB 10];                          % meta-d'
UB = [UB zeros(1,(nCriteria-1)/2)];    % criteria lower than t1c
UB = [UB 20*ones(1,(nCriteria-1)/2)];  % criteria higher than t1c

%% select constant criterion type

constant_criterion = 'meta_d1 * (t1c1 / d1)'; % relative criterion

%% set up initial guess at parameter values

ratingHR  = [];
ratingFAR = [];
for c = 2:nRatings*2
    ratingHR(end+1) = sum(nR_S2(c:end)) / sum(nR_S2);
    ratingFAR(end+1) = sum(nR_S1(c:end)) / sum(nR_S1);
end

t1_index = nRatings;
t2_index = setdiff(1:2*nRatings-1, t1_index);

d1 = (1/s) * fninv( ratingHR(t1_index) ) - fninv( ratingFAR(t1_index) );
meta_d1 = d1;
c = -0.5 .* (fninv(ratingHR(t1_index)) + fninv(ratingFAR(t1_index)));
c1 = (-1/(1+s)) * ( fninv( ratingHR ) + fninv( ratingFAR ) );
t1c1 = c1(t1_index);
t2c1 = c1(t2_index);

guess = [meta_d1 t2c1 - eval(constant_criterion)];

%% find the best fit for type 2 hits and FAs

% save fit_meta_d_MLE.mat nR_S1 nR_S2 t2FAR_rS2 t2HR_rS2 t2FAR_rS1 t2HR_rS1 nRatings t1c1 s d1 fncdf constant_criterion
save fit_meta_d_MLE.mat nR_S1 nR_S2 nRatings d1 t1c1 s constant_criterion fncdf fninv

op = optimset(@fmincon);
op = optimset(op,'MaxFunEvals',100000);

[x f] = fmincon(@fit_meta_d_logL,guess,A,b,[],[],LB,UB,[],op);

meta_d1  = x(1);
t2c1     = x(2:end) + eval(constant_criterion);
logL     = -f;

%% data is fit, now to package it...

%% find observed t2FAR and t2HR

% I_nR and C_nR are rating trial counts for incorrect and correct trials
% element i corresponds to # (in)correct w/ rating i
I_nR_rS2 = nR_S1(nRatings+1:end);
I_nR_rS1 = nR_S2(nRatings:-1:1);

C_nR_rS2 = nR_S2(nRatings+1:end);
C_nR_rS1 = nR_S1(nRatings:-1:1);

for i = 2:nRatings
    obs_FAR2_rS2(i-1) = sum( I_nR_rS2(i:end) ) / sum(I_nR_rS2);
    obs_HR2_rS2(i-1)  = sum( C_nR_rS2(i:end) ) / sum(C_nR_rS2);

    obs_FAR2_rS1(i-1) = sum( I_nR_rS1(i:end) ) / sum(I_nR_rS1);
    obs_HR2_rS1(i-1)  = sum( C_nR_rS1(i:end) ) / sum(C_nR_rS1);
end

%% find estimated t2FAR and t2HR

S1mu = -meta_d1/2; S1sd = 1;
S2mu =  meta_d1/2; S2sd = S1sd/s;

mt1c1 = eval(constant_criterion);

C_area_rS2 = 1-fncdf(mt1c1,S2mu,S2sd);
I_area_rS2 = 1-fncdf(mt1c1,S1mu,S1sd);

C_area_rS1 = fncdf(mt1c1,S1mu,S1sd);
I_area_rS1 = fncdf(mt1c1,S2mu,S2sd);

for i=1:nRatings-1

    t2c1_lower = t2c1(nRatings-i);
    t2c1_upper = t2c1(nRatings-1+i);

    I_FAR_area_rS2 = 1-fncdf(t2c1_upper,S1mu,S1sd);
    C_HR_area_rS2  = 1-fncdf(t2c1_upper,S2mu,S2sd);

    I_FAR_area_rS1 = fncdf(t2c1_lower,S2mu,S2sd);
    C_HR_area_rS1  = fncdf(t2c1_lower,S1mu,S1sd);

    est_FAR2_rS2(i) = I_FAR_area_rS2 / I_area_rS2;
    est_HR2_rS2(i)  = C_HR_area_rS2 / C_area_rS2;

    est_FAR2_rS1(i) = I_FAR_area_rS1 / I_area_rS1;
    est_HR2_rS1(i)  = C_HR_area_rS1 / C_area_rS1;

end

%% package output

fit.d1        = sqrt(2/(1+s^2)) * s * d1;
fit.c1        = c;
fit.s         = s;
fit.meta_d    = sqrt(2/(1+s^2)) * s * meta_d1;
fit.M_diff    = fit.meta_d - fit.d1;
fit.M_ratio   = fit.meta_d / fit.d1;

mt1c1         = eval(constant_criterion);
fit.meta_c1   = ( sqrt(2).*s ./ sqrt(1+s.^2) ) .* mt1c1;

t2ca          = ( sqrt(2).*s ./ sqrt(1+s.^2) ) .* t2c1;
fit.t2ca_rS1  = t2ca(1:nRatings-1);
fit.t2ca_rS2  = t2ca(nRatings:end);

fit.S1units.d1        = d1;
fit.S1units.meta_d1   = meta_d1;
fit.S1units.s         = s;
fit.S1units.meta_c1   = mt1c1;
fit.S1units.t2c1_rS1  = t2c1(1:nRatings-1);
fit.S1units.t2c1_rS2  = t2c1(nRatings:end);

fit.logL           = logL;

fit.est_HR2_rS1  = est_HR2_rS1;
fit.obs_HR2_rS1  = obs_HR2_rS1;

fit.est_FAR2_rS1 = est_FAR2_rS1;
fit.obs_FAR2_rS1 = obs_FAR2_rS1;

fit.est_HR2_rS2  = est_HR2_rS2;
fit.obs_HR2_rS2  = obs_HR2_rS2;

fit.est_FAR2_rS2 = est_FAR2_rS2;
fit.obs_FAR2_rS2 = obs_FAR2_rS2;

%% clean up
delete fit_meta_d_MLE.mat

end

%% function to find the likelihood of parameter values, given observed data
function logL = fit_meta_d_logL(parameters)

% set up parameters
meta_d1  = parameters(1);
t2c1     = parameters(2:end);

% loads:
% nR_S1 nR_S2 nRatings d1 t1c1 s constant_criterion fncdf fninv
load fit_meta_d_MLE.mat

% define mean and SD of S1 and S2 distributions
S1mu = -meta_d1/2; S1sd = 1;
S2mu =  meta_d1/2; S2sd = S1sd/s;

% adjust so that the type 1 criterion is set at 0
% (this is just to work with optimization toolbox constraints...
%  to simplify defining the upper and lower bounds of type 2 criteria)
S1mu = S1mu - eval(constant_criterion);
S2mu = S2mu - eval(constant_criterion);

t1c1 = 0;

%%% set up MLE analysis

% get type 2 response counts
for i = 1:nRatings

    % S1 responses
    nC_rS1(i) = nR_S1(i);
    nI_rS1(i) = nR_S2(i);

    % S2 responses
    nC_rS2(i) = nR_S2(nRatings+i);
    nI_rS2(i) = nR_S1(nRatings+i);

end

% get type 2 probabilities
C_area_rS1 = fncdf(t1c1,S1mu,S1sd);
I_area_rS1 = fncdf(t1c1,S2mu,S2sd);

C_area_rS2 = 1-fncdf(t1c1,S2mu,S2sd);
I_area_rS2 = 1-fncdf(t1c1,S1mu,S1sd);

t2c1x = [-Inf t2c1(1:nRatings-1) t1c1 t2c1(nRatings:end) Inf];

for i = 1:nRatings
    prC_rS1(i) = ( fncdf(t2c1x(i+1),S1mu,S1sd) - fncdf(t2c1x(i),S1mu,S1sd) ) / C_area_rS1;
    prI_rS1(i) = ( fncdf(t2c1x(i+1),S2mu,S2sd) - fncdf(t2c1x(i),S2mu,S2sd) ) / I_area_rS1;

    prC_rS2(i) = ( (1-fncdf(t2c1x(nRatings+i),S2mu,S2sd)) - (1-fncdf(t2c1x(nRatings+i+1),S2mu,S2sd)) ) / C_area_rS2;
    prI_rS2(i) = ( (1-fncdf(t2c1x(nRatings+i),S1mu,S1sd)) - (1-fncdf(t2c1x(nRatings+i+1),S1mu,S1sd)) ) / I_area_rS2;
end

% calculate logL
logL = 0;
for i = 1:nRatings

    logL = logL + nC_rS1(i)*log(prC_rS1(i)) + nI_rS1(i)*log(prI_rS1(i)) + ...
                  nC_rS2(i)*log(prC_rS2(i)) + nI_rS2(i)*log(prI_rS2(i));

end

if isnan(logL), logL=-Inf; end

logL = -logL;

end

author: translationalneuromodeling

repo: tapas


text:
function obj = Make_f11(obj,D_filename,file_suffixes,bathyfile)
% obj = Make_f11(obj,D_filename,file_suffixes)
% Input a msh class object get the values of the density over the depth
% based on D_filename, and computes the depth-averaged value which
% populates the f11 struct in the msh class object.
%
%  Author:      William Pringle
%  Created:     March 19 2018
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if isempty(obj.b)
    error('No bathymetry data in grid to calculate the depth-averaged value')
end
rho0 = 1000;
%g = 9.807;

% %% Do some projection for obj.p
% proj = 'Mercator';
% R = 6378206.4; %[m]
% m_proj(proj,'lon',[ min(obj.p(:,1)) max(obj.p(:,1)) ],...
%             'lat',[ min(obj.p(:,2)) max(obj.p(:,2))])
% [xx,yy] = m_ll2xy(obj.p(:,1),obj.p(:,2));
% xx = R*xx; yy= R*yy;
%
% %% Get the element areas and slopes, connectivity table
% A = polyarea(xx(obj.t(:,1:3))',yy(obj.t(:,1:3))')';
% % Get x differences
% a = [ xx(obj.t(:,3)) - xx(obj.t(:,2))
%       xx(obj.t(:,1)) - xx(obj.t(:,3))
%       xx(obj.t(:,2)) - xx(obj.t(:,1)) ] ;
% a = reshape(a,[],3);
%
% % Get y differences
% b = [ yy(obj.t(:,2)) - yy(obj.t(:,3))
%       yy(obj.t(:,3)) - yy(obj.t(:,1))
%       yy(obj.t(:,1)) - yy(obj.t(:,2)) ];
% b = reshape(b,[],3);
%
% % Get the vertex to element table
% vtoe = VertToEle(obj.t);
% vtoe(vtoe == 0) = length(obj.t) + 1;
% A(end+1) = 0;
% An = sum(A(vtoe))';

% nearest neighbour extrapolation of T and S (0 for no extrapolation)?
FillNaN = 0;

[direc , title , ext ] = fileparts( D_filename);

if isempty(title)
    D_filename = dir(direc);
    D_filename(1:2) = [];
    idx = zeros(length(D_filename),1);
    for jj = 1:size(file_suffixes,1)
        for ii = 1:length(D_filename)
            if idx(ii) == 0
                idx(ii) = contains(D_filename(ii).name,file_suffixes(jj,:));
            end
        end
    end
    D_filename = D_filename(find(idx));
else
   name = D_filename; D_filename = [];
   D_filename.folder = direc;
   D_filename.name = [title ext];
end

Bx = NaN(length(obj.b),size(D_filename,1));
By = NaN(length(obj.b),size(D_filename,1));
tic
for ii = 1:size(D_filename,1)
    Dn = [D_filename(ii).folder '/' D_filename(ii).name];
    [~ , ~ , ext ] = fileparts( Dn);
    %% Read the grid and density data
    if strcmp(ext,'.nc')
        try
            title = ncreadatt(Dn,'/','title'); % Reading some
            lon = ncread(Dn,'lon');
            lat = ncread(Dn,'lat');
            z   = ncread(Dn,'depth');
            sigma_t = ncread(Dn,'I_an');
            time_all = 0;
        catch
            title = ncreadatt(Dn,'/','History'); % Reading some
            [time_all(ii),lon,lat,z,~,~,dpx,dpy] = Make_Gridded_rho(Dn,...
         {'time','depth','lon','lat','salinity','water_temp','surf_el'}); %,...
                                                                %bathyfile);
            disp(['Read time ' datestr(time_all(ii))])
            %if ii == 1
            %    sigma_t = zeros([size(rho) size(D_filename,1)]);
            %end
            %sigma_t(:,:,:,ii) = rho - 1000; clear rho
        end
    elseif strcmp(ext,'.mat')
        load(Dn);
        %sigma_t = rho - 1000; clear rho
    else
        error('does not understand file extension')
    end
%     % zeta on the finite-element mesh
%     [Lon,Lat] = ndgrid(lon,lat);
%     F = griddedInterpolant(Lon,Lat,zeta,'linear','none');
%     zeta = F(obj.p);
%     % Get rho on the finite-element mesh
%     [~,~,~,rho] = Gridded_to_Mesh_SeaBed_DepthAveraged(...
%             obj.p(:,1),obj.p(:,2),obj.b,z,rho,lon,lat,FillNaN);
%
%     % Compute bx, by for each element (loop over to reduce memory issues)
%     bxe = zeros(length(obj.t),length(z));
%     bye = zeros(length(obj.t),length(z));
%     for k = 1:length(z)
%         if k == 1
%             BCPe = zeta.*rho{k};
%         else
%             BCPe = rho{k} - rho0;
%         end
%         BCPe = BCPe(obj.t);
%         bxe(:,k) = 0.5 * sum(BCPe.*b,2);
%         bye(:,k) = 0.5 * sum(BCPe.*a,2);
%     end
%     bxe(end+1,:) = 0; bye(end+1,:) = 0;
%     bx = zeros(length(obj.b),length(z));
%     by = zeros(length(obj.b),length(z));
%     for k = length(z):-1:1
%         if k > 1
%             bxe(:,k) = trapz(z(1:k),bxe(:,1:k),2);
%             bye(:,k) = trapz(z(1:k),bye(:,1:k),2);
%         end
%         temp = sum(reshape(bxe(vtoe,k),size(vtoe,1),[]),'omitnan')';
%         bx(:,k) = temp./An;
%         by(:,k) = temp./An;
%     end

    % Interpolate the gradient at each depth
    dpx(isnan(dpy)) = NaN; dpy(isnan(dpx)) = NaN;
    [lonN,latN] = ndgrid(lon,lat);
    bx = NaN(length(obj.p),length(z));
    by = NaN(length(obj.p),length(z));
    for kk = 1:length(z)
        if kk == 1
            Fx = griddedInterpolant(lonN,latN,dpx(:,:,kk),'linear','none');
            Fy = griddedInterpolant(lonN,latN,dpy(:,:,kk),'linear','none');
        else
            Fx.Values = dpx(:,:,kk);
            Fy.Values = dpy(:,:,kk);
        end
        bb = obj.b >= z(kk);
        pt = obj.p(bb,:);
        bxt = Fx(obj.p(bb,:));
        byt = Fy(obj.p(bb,:));
        if ~isempty(find(isnan(bxt), 1))
            idn = find(isnan(bxt));
            idg = find(~isnan(bxt));
            idx = knnsearch(pt(idg,:),pt(idn,:));
            bxt(idn) = bxt(idg(idx));
        end
        if ~isempty(find(isnan(byt), 1))
            idn1 = find(isnan(byt));
            idg = find(~isnan(byt));
            if sum(idn - idn1) ~= 0
                idx = knnsearch(pt(idg,:),pt(idn1,:));
            end
            byt(idn1) = byt(idg(idx));
        end
        bx(bb,kk) = bxt;
        by(bb,kk) = byt;
    end
    % Now do the integration over the depth taking into account actual
    % depth on the computational mesh
    for kk = length(z):-1:1
        bb = isnan(Bx(:,ii)) & obj.b > z(kk) & ...
            ~isnan(bx(:,kk)) & ~isnan(by(:,kk));
        Bx(bb,ii) = trapz(z(1:kk),bx(bb,1:kk),2) + ...
                          bx(bb,kk).*(obj.b(bb) - z(kk));
        By(bb,ii) = trapz(z(1:kk),by(bb,1:kk),2) + ...
                          by(bb,kk).*(obj.b(bb) - z(kk));
    end
    % Get the depth averaged
    Bx(:,ii) = Bx(:,ii)./obj.b/rho0;
    By(:,ii) = By(:,ii)./obj.b/rho0;

%     figure;
%     fastscatter(obj.p(:,1),obj.p(:,2),9.807*hypot(Bx,By))
%     colormap(cmocean('speed'))
%     caxis([0 1e-4])
%     colorbar;
end
toc
disp(time_all)
DTIMINC = round(seconds(median(diff(time_all))));

%
% tic
% Bx = zeros(length(obj.b),size(sigma_t,4));
% By = zeros(length(obj.b),size(sigma_t,4));
% for t = 1:size(sigma_t,4)
%     disp(['Computing ' datestr(time_all(t))])
%     %% Calculate the depth-averaged density
%     [~,Sigma_tm,~,Sigma_3D] = Gridded_to_Mesh_SeaBed_DepthAveraged(...
%             obj.p(:,1),obj.p(:,2),obj.b,z,sigma_t(:,:,:,t),lon,lat,FillNaN);
%
%     %% Calculate the integrated barcolinic pressure's
%     Sigma_3D = reshape(cell2mat(Sigma_3D),[],length(Sigma_3D));
%     % Make sure derivative with "ground" returns NaN;
%     BCP = NaN(length(obj.b),length(z));
%     for k = 1:length(z)
%         bb = obj.b >= z(k);
%         BCP(bb,k) = trapz(z(1:k),Sigma_3D(bb,1:k),2);
%     end
%     BCP = BCP/rho0;
%
%     %% Calculate the bx/by's for each element
%     % Compute bx, by for each element (loop over to reduce memory issues)
%     bxe = zeros(length(obj.t),length(z));
%     bye = zeros(length(obj.t),length(z));
%     for k = 1:length(z)
%         BCPe = BCP(:,k);
%         BCPe = BCPe(obj.t);
%         bxe(:,k) = 0.5 * sum(BCPe.*b,2);
%         bye(:,k) = 0.5 * sum(BCPe.*a,2);
%     end
%     %% Sum the contributions of each element connected to a node
%     % Add in a zero ghost element and refer to the zero indices of vtoe to this
%     bxe(end+1,:) = 0; bye(end+1,:) = 0;
%
%     % Sum up all element contributions to each node and take the integral
%     bx = zeros(length(obj.b),length(z));
%     by = zeros(length(obj.b),length(z));
%     for k = 1:length(z)
%         bb = obj.b >= z(k);
%         temp = sum(reshape(bxe(vtoe,k),size(vtoe,1),[]),'omitnan')';
%         bx(bb,k) = temp(bb)./An(bb);
%         temp = sum(reshape(bye(vtoe,k),size(vtoe,1),[]),'omitnan')';
%         by(bb,k) = temp(bb)./An(bb);
%     end
%     Bx(:,t) = trapz(z,bx,2)./obj.b;
%     By(:,t) = trapz(z,by,2)./obj.b;
%     Bx(obj.b == 0,t) = NaN; By(obj.b == 0,t) = NaN;
% %     % On the structured grid
%     %% Calculate the integrated barcolinic pressure's
%     Sigma_3D = sigma_t(:,:,:,t);
%     % Make sure derivative with "ground" returns NaN;
%     BCP = NaN(size(Sigma_3D));
%     for k = 1:length(z)
%         BCP(:,:,k) = trapz(z(1:k),Sigma_3D(:,:,1:k),3);
%     end
%     BCP = BCP/rho0;
%
%     bx = zeros(size(Sigma_3D));
%     by = zeros(size(Sigma_3D));
%     for k = 1:length(z)
%         [bx(:,:,k),by(:,:,k)] = gradient(BCP(:,:,k),(lon(2)-lon(1))*111e3,(lat(2)-lat(1))*111e3);
%     end
%     bxv = reshape(bx,[],length(z));
%     byv = reshape(by,[],length(z));
%     Bx_struc = NaN(size(bxv,1),1);
%     By_struc = NaN(size(bxv,1),1);
%     for k = 1:length(z)
%         bb = ~isnan(bxv(:,k));
%         Bx_struc(bb) = trapz(z(1:k),bxv(bb,1:k),2)./z(k);
%         bb = ~isnan(byv(:,k));
%         By_struc(bb) = trapz(z(1:k),byv(bb,1:k),2)./z(k);
%     end
%     Bx_struc(Bx_struc == 0) = NaN;
%     By_struc(By_struc == 0) = NaN;
%     Bx_struc = reshape(Bx_struc,size(Sigma_3D,1),[])';
%     By_struc = reshape(By_struc,size(Sigma_3D,1),[])';
%     figure;
%     pcolor(hypot(Bx_struc,By_struc))
%     shading interp
%     colormap(cmocean('ice'))
%     caxis([0 2e-6])
%     colorbar;
%end
%toc
% Calculate the depth-averaged
% rhoxe = sum(Sigma_tm(obj.t).*b,2)/rho0;
% rhoye = sum(Sigma_tm(obj.t).*a,2)/rho0;
% rhoxe(end+1,:) = 0; rhoye(end+1,:) = 0;
% rhox = sum(rhoxe(vtoe))'./An;
% rhoy = sum(rhoye(vtoe))'./An;
% Bx2D = 0.5*obj.b.^2.*rhox;
% By2D = 0.5*obj.b.^2.*rhoy;

%% Make into f11 struct
%summary = ncreadatt(D_filename,'/','summary'); % global attributes
obj.f11.DataTitle = title;
%[~,name,~] = fileparts(D_filename);
%obj.f11.DataSubTitle = [summary ': ' name];
obj.f11.DTIMINC = DTIMINC;
obj.f11.TimeVec = time_all;
obj.f11.NumOfNodes = length(obj.b);
obj.f11.Val = cell(size(Bx,2),1); %zeros(3,length(obj.b),size(sigma_t,4));
for t = 1:size(Bx,2)
    idx = find(~isnan(Bx(:,t)) & ~isnan(By(:,t)));
    obj.f11.Val{t} = [idx'; Bx(idx,t)'; By(idx,t)']; % Sigma_tm';
end
%EOF
end

author: CHLNDDEV


text:
%%% [3dmatrix] = convert_video_to_3d(struct)
%
function [V] = convert_video_to_3d(video)
  for i = 1:video.nrFramesTotal
    frame = video.frames(i).cdata;
    if(size(frame,3) == 3)
      frame = rgb2gray(frame);
    end
    V(:,:,i) = frame;
  end
end

author: andrewssobral

repo: lrslibrary