text:
function err = getH1errorbd(node,Neumann,pde,u,A,quadorder)

if ~exist('quadorder','var'), quadorder = 4; end

el = sqrt(sum((node(Neumann(:,1),:) - node(Neumann(:,2),:)).^2,2));
[lambdagN,weightgN] = quadpts1(quadorder);
phigN = lambdagN;                 % linear bases
nQuadgN = size(lambdagN,1);
err = zeros(size(Neumann,1),1);
for pp = 1:nQuadgN
    % quadrature points in the x-y coordinate
    ppxy = lambdagN(pp,1)*node(Neumann(:,1),:) ...
         + lambdagN(pp,2)*node(Neumann(:,2),:);
    uhp = u(Neumann(:,1))*phigN(pp,1) + u(Neumann(:,2))*phigN(pp,2);
    up = pde.exactu(ppxy);
    gNp = pde.g_N(ppxy);
    err = err + weightgN(pp)*gNp.*(up - 2*uhp);
end
err = sum(err.*el) + u'*A*u;
err = sqrt(abs(err));

repo: ifem


text:
classdef SmoothingExponentComputer < handle

    properties (Access = protected)
        value
    end

    methods (Access = public, Static)

        function obj = create(cParams)
            f = SmoothingExponentComputerFactory();
            obj = f.create(cParams);
        end

    end

    methods (Access = public)

        function q = compute(obj)
           obj.computeExponent();
           q = obj.value;
        end

    end

    methods (Access = protected, Abstract)
       computeExponent(obj);
    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function r=ksrlin(x,y,h,N)
% KSRLIN   Local linear kernel smoothing regression
%
% r=ksrlin(x,y) returns the local linear Gaussian kernel regression in
% structure r such that r.f(r.x) = y(x) + e. The bandwidth and number of
% samples are also stored in r.h and r.n respectively.
%
% r=ksrlin(x,y,h) performs the regression using the specified bandwidth, h.
%
% r=ksrlin(x,y,h,n) calculates the regression in n points (default n=100).
%
% Without output, ksrlin(x,y) or ksrlin(x,y,h) will display the regression
% plot.
%
% Algorithm
% The kernel regression is a non-parametric approach to
% estimate the conditional expectation of a random variable:
%
% E(Y|X) = f(X)
%
% where f is a non-parametric function. The normal kernel regression is a
% local constant estimator. The extension of local linear estimator is
% obtained by solving the least squares problem:
%
% min sum (y-alpha-beta(x-X))^2 kerf((x-X)/h)
%
% The local linear estimator can be given an explicit formula
%
% f(x) = 1/n sum((s2-s1*(x-X)).*kerf((x-X)/h).*Y)/(s2*s0-s1^2)
%
%where si = sum((x-X)^i*kerf((x-X)/h))/n. Compare with local constant
%estimator, the local linear estimator can improve the estimation near the
%edge of the region over which the data have been collected.
%
% See also gkde, ksdensity, ksr

% Reference:
% Bowman, A.W. and Azzalini, A., Applied Smoothing Techniques for Data
% Analysis, Clarendon Press, 1997 (p.50)

% Example 1: smooth curve with noise
%{
x = 1:100;
y = sin(x/10)+(x/50).^2;
yn = y + 0.2*randn(1,100);
r=ksrlin(x,yn);
r1=ksr(x,yn); % downloadable from FEX ID 19195
plot(x,y,'b-',x,yn,'co',r.x,r.f,'r--',r1.x,r1.f,'m-.','linewidth',2)
legend('true','data','local linear','local constant','location','northwest');
title('Gaussian kernel regression')
%}
% Example 2: with missing data
%{
x = sort(rand(1,100)*99)+1;
y = sin(x/10)+(x/50).^2;
y(round(rand(1,20)*100)) = NaN;
yn = y + 0.2*randn(1,100);
r=ksrlin(x,yn);
r1=ksr(x,yn); % downloadable from FEX ID 19195
plot(x,y,'b-',x,yn,'co',r.x,r.f,'r--',r1.x,r1.f,'m-.','linewidth',2)
legend('true','data','local linear','local constant','location','northwest');
title('Gaussian kernel regression with 20% missing data')
%}

% By Yi Cao at Cranfield University on 12 April 2008.
%

% Check input and output
error(nargchk(2,4,nargin));
error(nargoutchk(0,1,nargout));
if numel(x)~=numel(y)
    error('x and y are in different sizes.');
end

x=x(:);
y=y(:);
% clean missing or invalid data points
inv=(x~=x)|(y~=y)|(abs(x)==Inf)|(abs(y)==Inf);
x(inv)=[];
y(inv)=[];

% Default parameters
if nargin<4
    N=100;
elseif ~isscalar(N)
    error('N must be a scalar.')
end
r.n=length(x);
if nargin<3
    % optimal bandwidth suggested by Bowman and Azzalini (1997) p.31
    hx=median(abs(x-median(x)))/0.6745*(4/3/r.n)^0.2;
    hy=median(abs(y-median(y)))/0.6745*(4/3/r.n)^0.2;
    h=sqrt(hy*hx);
elseif ~isscalar(h)
    error('h must be a scalar.')
end
r.h=h;

% Gaussian kernel function
kerf=@(z)exp(-z.*z/2)/sqrt(2*pi);

r.x=linspace(min(x),max(x),N);
r.f=zeros(1,N);
for k=1:N
    d=r.x(k)-x;
    z=kerf(d/h);
    s1=d.*z;
    s2=sum(d.*s1);
    s1=sum(s1);
    r.f(k)=sum((s2-s1*d).*z.*y)/(s2*sum(z)-s1*s1);
end

% Plot
if ~nargout
    plot(r.x,r.f,'r',x,y,'bo')
    ylabel('f(x)')
    xlabel('x')
    title('Kernel Smoothing Regression');
end

author: Sable

repo: mcbench-benchmarks


text:
function experiments = stack_vottir2016()

set_global_variable('bundle', 'http://data.votchallenge.net/vot2016/vot-tir2016.zip');

baseline.name = 'baseline';
baseline.converter = [];
baseline.type = 'supervised';
baseline.tags = {'camera_motion', 'dynamics_change', 'occlusion', 'size_change', 'motion_change', 'empty'};

baseline.parameters.repetitions = 15;
baseline.parameters.burnin = 10;
baseline.parameters.skip_initialize = 5;
baseline.parameters.failure_overlap = 0;

experiments = {baseline};

end

author: votchallenge

repo: toolkit-legacy


text:
function xhat = kalman(z, A, C, R, Q)
% Simple Kalman Filter (linear) with optimal gain, no control signal
%
% z Measurement signal              m observations X # of observations
% A State transition model          n X n, n = # of state values
% C Observation model               m X n
% R Covariance of observation noise m X m
% Q Covariance of process noise     n X n
%
% Based on http://en.wikipedia.org/wiki/Kalman_filter, but matrices named
% A, C, G instead of F, H, K.
%
% See http://home.wlu.edu/~levys/kalman_tutorial for background
%
% Copyright (C) 2014 Simon D. Levy
%
% This code is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as
% published by the Free Software Foundation, either version 3 of the
% License, or (at your option) any later version.
%
% This code is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with this code. If not, see <http:#www.gnu.org/licenses/>.

% Initializtion =====================================================

% Number of sensors
m = size(C, 1);

% Number of state values
n = size(C, 2);

% Number of observations
numobs = size(z, 2);

% Use linear least squares to estimate initial state from initial
% observation
xhat = zeros(n, numobs);
xhat(:,1) = C \ z(:,1);

% Initialize P, I
P = ones(size(A));
I = eye(size(A));

% Kalman Filter =====================================================

for k = 2:numobs

    % Predict
    xhat(:,k) = A * xhat(:,k-1);
    P         = A * P * A' + Q;

    % Update
    G         = P  * C' / (C * P * C' + R);
    P         = (I - G * C) * P;
    xhat(:,k) = xhat(:,k) + G * (z(:,k) - C * xhat(:,k));

end

author: simondlevy

repo: SensorFusion


text:
%%******************************************************************
%% geometric_mean: an example with rotated cone variables
%%
%%  max {prod(d+B*x) : d+B*x > 0, x <= 10}
%%
%%  where B = 4xn matrix,
%%        d = 4x1 vector
%%
%% [blk,At,C,b,xx] = geometric_mean(B,d,solve);
%%
%% E.g. p = 6; m = 4; B = [rand(2,p); -rand(2,p)]; d = rand(m,1);
%%
%% SDPT3: version 3.0
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last modified: 2 Feb 01
%%******************************************************************

   function [blk,At,C,b,xx] = geometric_mean(B,d,solve);

   if (nargin == 2); solve = 0; end

   n = size(B,2);
   zz = zeros(1,n);
   r2 = sqrt(2);
   blk{1,1} = 'r'; blk{1,2} = [3,3,3];
   At{1,1} = -[B(1,:),0,0,0; B(2,:),0,0,0; zz,r2,0,0; ...
               B(3,:),0,0,0; B(4,:),0,0,0; zz,0,r2,0; ...
	       zz, 1,0,0; zz, 0,1,0; zz, 0,0,r2];
   C{1,1} = [d(1:2); 0; d(3:4); 0; 0;0;0];
   b = [zeros(n,1); 0;0;1];
   blk{2,1} = 'l'; blk{2,2} = n;
   At{2,1} = [eye(n), zeros(n,3)]; C{2,1} = 10*ones(n,1);
   if (solve)
      [bblk,AAt,CC,bb,T] = convertRcone(blk,At,C,b);
      [obj,X,y,Z] = sqlp(bblk,AAt,CC,bb);
      xx = y(1:n);
   end
%%******************************************************************

repo: REST


text:
%------------------------------ PolyMesher -------------------------------%
% Ref: C Talischi, GH Paulino, A Pereira, IFM Menezes, "PolyMesher: A     %
%      general-purpose mesh generator for polygonal elements written in   %
%      Matlab," Struct Multidisc Optim, DOI 10.1007/s00158-011-0706-z     %
%-------------------------------------------------------------------------%
function [x] = MichellDomain(Demand,Arg)
  BdBox = [0 5 -2 2];
  switch(Demand)
    case('Dist');  x = DistFnc(Arg,BdBox);
    case('BC');    x = BndryCnds(Arg{:},BdBox);
    case('BdBox'); x = BdBox;
    case('PFix');  x = FixedPoints(BdBox);
  end
%----------------------------------------------- COMPUTE DISTANCE FUNCTIONS
function Dist = DistFnc(P,BdBox)
  d1 = dRectangle(P,BdBox(1),BdBox(2),BdBox(3),BdBox(4));
  d2 = dCircle(P,0,0,BdBox(4)/2);
  Dist = dDiff(d1,d2);
%---------------------------------------------- SPECIFY BOUNDARY CONDITIONS
function [x] = BndryCnds(Node,Element,BdBox)
  eps = 0.1*sqrt((BdBox(2)-BdBox(1))*(BdBox(4)-BdBox(3))/size(Node,1));
  CircleNodes = find(abs(sqrt(Node(:,1).^2+Node(:,2).^2)-1.0)<eps);
  Supp = ones(size(CircleNodes,1),3);
  Supp(:,1) = CircleNodes;
  MidRightFace = sqrt((Node(:,1)-BdBox(2)).^2+...
                      (Node(:,2)-(BdBox(3)+BdBox(4))/2).^2);
  [foo,MidRightFace] = sort(MidRightFace);
  Load = [MidRightFace(1),0,-1];
  x = {Supp,Load};
%----------------------------------------------------- SPECIFY FIXED POINTS
function [PFix] = FixedPoints(BdBox)
  PFix = [5 0];
%-------------------------------------------------------------------------%


text:
function [intra, inter, Qnodes, Fnodes, Onode] = mk_hhmm_topo_F1(D, all_Q_to_Qs, Ops)
% MK_HHMM_TOPO Make Hierarchical HMM topology assuming level 1 can finish
% function [intra, inter, Qnodes, Fnodes, Onode] = mk_hhmm_topo(D, all_Q_to_Qs, Ops, F1)
%
% D is the depth of the hierarchy
% If all_Q_to_Qs = 1, level i connects to all levels below, else just to i+1 [0]
% Ops are the Q parents of the observed node [Qnodes(end)]
% If F1=1, level 1 can finish (restart), else there is no F1->Q1 arc [0]

Qnodes = 1:D;

if nargin < 2, all_Q_to_Qs = 1; end
if nargin < 3, Ops = Qnodes(D); end
if nargin < 4, F1 = 0; end

if F1
  Fnodes = 2*D:-1:D+1; % must number from bottom to top
  Onode = 2*D+1;
  ss = 2*D+1;
else
  Fnodes = (2*D)-1:-1:D+1;
  Onode = 2*D;
  ss = 2*D;
end

intra = zeros(ss);
intra(Ops, Onode) = 1;
for d=1:D-1
  if all_Q_to_Qs
    intra(Qnodes(d), Qnodes(d+1:end)) = 1;
  else
    intra(Qnodes(d), Qnodes(d+1)) = 1;
  end
end
for d=D:-1:3
  intra(Fnodes(d), Fnodes(d-1)) = 1;
end
if F1
  intra(Fnodes(2), Fnodes(1)) = 1;
end
if all_Q_to_Qs
  for d=1:D
    intra(Qnodes(d), Fnodes(d:end)) = 1;
  end
else
  for d=1:D
    if d < D
      intra(Qnodes(d), Fnodes([d d+1])) = 1;
    else
      intra(Qnodes(d), Fnodes(d)) = 1;
    end
  end
end

inter = zeros(ss);
for d=1:D
  inter(Qnodes(d), Qnodes(d)) = 1;
end
for d=1:D
  if d==1
    inter(Fnodes(d), Qnodes(d)) = 1;
  else
    inter(Fnodes(d), Qnodes([d-1 d])) = 1;
  end
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [grad,isBdNode] = gradmatrix(edge,isBdEdge)
%% GRADMATRIX matrix for the gradient of a nodal linear element
%
% grad = gradmatrix(edge) returns the sparse matrix grad which is an NE by 3N
% matrix mapping linear nodal element (a vector with length N) to lowest
% order linear edge element (a vector with length NE). The trace is zero at
% certain boundary edges given by isBdEdge.
%
% grad = gradmatrix(edge) returns the corresponding matrix without boundary
% conditions.
%
% This function is used in mgMaxwell to construct transfer operators for
% the HX preconditioner.
%
% See also mgMaxwell, mgMaxwell1, mgMaxwell2, node2edgematrix1
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

if nargin<=1, isBdEdge = []; end
NE = size(edge,1); N  = double(max(edge(:)));
i = repmat((1:NE)',2,1);
j = double(edge(:));
s = [-ones(NE,1),ones(NE,1)];
if any(isBdEdge) % no grad on boundary edges
    bdEdge = edge(isBdEdge,:);
    isBdNode = false(N,1);
    isBdNode(bdEdge(:)) = true;
    idx = ~(isBdEdge(i) | isBdNode(j));
else
    idx = 1:length(i);
end
grad = sparse(i(idx),j(idx),s(idx),NE,N);


text:
%% ArUco Grid Board Image Demo
% This example shows how to create an ArUco grid board image.
%
% Sources:
%
% * <https://github.com/opencv/opencv_contrib/blob/3.1.0/modules/aruco/samples/create_board.cpp>
%

%% Board

% params
markersX = 5;                % Number of markers in X direction
markersY = 7;                % Number of markers in Y direction
markerLength = 60;           % Marker side length (in pixels)
markerSeparation = 15;       % Separation between two consecutive markers in the grid (in pixels)
dictionaryId = '6x6_250';    % dictionary id
margins = markerSeparation;  % Margins size (in pixels)
borderBits = 1;              % Number of bits in marker borders

imageSize = [markersY, markersX] .* ...
    (markerLength + markerSeparation) - markerSeparation + 2 * margins;

% create board
dictionary = {'Predefined', dictionaryId};
board = {'GridBoard', ...
    markersX, markersY, markerLength, markerSeparation, dictionary};

%% Display

% show created board
boardImage = cv.drawPlanarBoard(board, fliplr(imageSize), ...
    'MarginSize',margins, 'BorderBits',borderBits);
imshow(boardImage), title('GridBoard')

% save image
imwrite(boardImage, fullfile(tempdir(),'GridBoard.png'));

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
% Author: Ricardo Baptista and Matthias Poloczek
% Date:   June 2018
%
% See LICENSE.md for copyright information
%

function Q = rand_ising_grid(n_vars)
% RAND_ISING: Function generates a random interaction matrix
% for an Ising Model on a grid graph with n_vars total variables

% Check that n_side is an integer
n_side = sqrt(n_vars);
if floor(n_side) ~= n_side
	error('Number of nodes is not square')
end

% Connect nodes horizontally
for i=1:n_side
	for j=1:n_side-1

		% Determine node idx
		node = (i-1)*n_side + j;

		Q(node,node+1) = 4.95*rand() + 0.05;%0.95*rand() + 0.05;
		Q(node+1,node) = Q(node,node+1);

	end
end

% Connect nodes vertically
for i=1:n_side-1
	for j=1:n_side

		% Determine node idx
		node = (i-1)*n_side + j;

		Q(node,node+n_side) = 4.95*rand() + 0.05;%0.95*rand() + 0.05;
		Q(node+n_side,node) = Q(node,node+n_side);

	end
end

% Apply random sign flips to Q
rand_sign = tril((rand(n_vars,n_vars) > 0.5)*2 - 1,-1);
rand_sign = rand_sign + rand_sign';
Q = rand_sign.*Q;

end

author: baptistar

repo: BOCS

save_path: github-repos/MATLAB/baptistar-BOCS


text:
function img = readFlowFile(filename)

% readFlowFile read a flow file FILENAME into 2-band image IMG

%   According to the c++ source code of Daniel Scharstein
%   Contact: schar@middlebury.edu

%   Author: Deqing Sun, Department of Computer Science, Brown University
%   Contact: dqsun@cs.brown.edu
%   $Date: 2007-10-31 16:45:40 (Wed, 31 Oct 2006) $

% Copyright 2007-2010, Brown University, Providence, RI. USA
%
%
%                         All Rights Reserved
%
% Permission to use, copy, modify, and distribute this software and its
% documentation for any purpose other than its incorporation into a
% commercial product is hereby granted without fee, provided that the
% above copyright notice appear in all copies and that both that
% copyright notice and this permission notice appear in supporting
% documentation, and that the name of the author and Brown University not be used in
% advertising or publicity pertaining to distribution of the software
% without specific, written prior permission.
%
% THE AUTHOR AND BROWN UNIVERSITY DISCLAIM ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
% INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR ANY
% PARTICULAR PURPOSE.  IN NO EVENT SHALL THE AUTHOR OR BROWN UNIVERSITY BE LIABLE FOR
% ANY SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
% WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
% ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
% OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.

TAG_FLOAT = 202021.25;  % check for this when READING the file

% sanity check
if isempty(filename) == 1
    error('readFlowFile: empty filename');
end;

idx = findstr(filename, '.');
idx = idx(end);

if length(filename(idx:end)) == 1
    error('readFlowFile: extension required in filename %s', filename);
end;

if strcmp(filename(idx:end), '.flo') ~= 1
    error('readFlowFile: filename %s should have extension ''.flo''', filename);
end;

fid = fopen(filename, 'r');
if (fid < 0)
    error('readFlowFile: could not open %s', filename);
end;

tag     = fread(fid, 1, 'float32');
width   = fread(fid, 1, 'int32');
height  = fread(fid, 1, 'int32');

% sanity check

if (tag ~= TAG_FLOAT)
    error('readFlowFile(%s): wrong tag (possibly due to big-endian machine?)', filename);
end;

if (width < 1 || width > 99999)
    error('readFlowFile(%s): illegal width %d', filename, width);
end;

if (height < 1 || height > 99999)
    error('readFlowFile(%s): illegal height %d', filename, height);
end;

nBands = 2;

% arrange into matrix form
tmp = fread(fid, inf, 'float32');
tmp = reshape(tmp, [width*nBands, height]);
tmp = tmp';
img(:,:,1) = tmp(:, (1:width)*nBands-1);
img(:,:,2) = tmp(:, (1:width)*nBands);

fclose(fid);

author: kristinbranson

repo: JAABA


text:
function X = skew(x)

X=[0 -x(3) x(2) ; x(3) 0 -x(1) ; -x(2) x(1) 0 ];

end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function file_exist=my_check_file(cache_file, be_quite)

if nargin<2
    be_quite=false;
end

if exist(cache_file, 'file')

    file_exist=true;

    if ~be_quite
        fprintf('file exist:%s\n',cache_file);
    end
else

    file_exist=false;

    if ~be_quite
        fprintf('file not exist:%s\n',cache_file);
    end
end

end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function [h, varargout] = nethess_weighted(w, net, x, t, eso_w, varargin)
%NETHESS Evaluate network Hessian
%
%	Description
%
%	H = NETHESS(W, NET, X, T) takes a weight vector W and a network data
%	structure NET, together with the matrix X of input vectors and the
%	matrix T of target vectors, and returns the value of the Hessian
%	evaluated at W.
%
%	[E, VARARGOUT] = NETHESS(W, NET, X, T, VARARGIN) also returns any
%	additional return values from the network Hessian function, and
%	passes additional arguments to that function.
%
%	See also
%	NETERR, NETGRAD, NETOPT
%

%	Copyright (c) Ian T Nabney (1996-9)

hess_str = [net.type, 'hess_weighted'];

net = netunpak(net, w);

[s{1:nargout}] = feval(hess_str, net, x, t, eso_w, varargin{:});
h = s{1};
for i = 2:nargout
  varargout{i-1} = s{i};
end


text:
% FIGURE Create figure window.
%    FIGURE, by itself, creates a new figure window, and returns
%    its handle.
%
%    FIGURE(H) makes H the current figure, forces it to become visible,
%    and raises it above all other figures on the screen.  If Figure H
%    does not exist, and H is an integer, a new figure is created with
%    handle H.
%
%    GCF returns the handle to the current figure.
%
%    Execute GET(H) to see a list of figure properties and
%    their current values. Execute SET(H) to see a list of figure
%    properties and their possible values.
%
%    See also SUBPLOT, AXES, GCF, CLF.
%
%    Reference page in Doc Center
%       doc figure
%
%

author: jmaih


text:
function [nll,g,H,T] = LogisticLoss(w,X,y)
% w(feature,1)
% X(instance,feature)
% y(instance,1)

[n,p] = size(X);

Xw = X*w;
yXw = y.*Xw;

nll = sum(mylogsumexp([zeros(n,1) -yXw]));

if nargout > 1
    if nargout > 2
        sig = 1./(1+exp(-yXw));
        g = -X.'*(y.*(1-sig));
    else
        g = -X.'*(y./(1+exp(yXw)));
    end
end

if nargout > 2
    H = X.'*diag(sparse(sig.*(1-sig)))*X;
end

if nargout > 3
    T = zeros(p,p,p);
    for j1 = 1:p
        for j2 = 1:p
            for j3 = 1:p
                T(j1,j2,j3) = sum(y(:).^3.*X(:,j1).*X(:,j2).*X(:,j3).*sig.*(1-sig).*(1-2*sig));
            end
        end
    end
end

author: huashiyiqike

repo: LSTM-MATLAB


text:
function D = driving_function_mono_wfs(x0,xs,src,f,conf)
%DRIVING_FUNCTION_MONO_WFS driving signal for WFS
%
%   Usage: D = driving_function_mono_wfs(x0,xs,src,f,conf)
%
%   Input parameters:
%       x0          - position and direction of the secondary source / m [nx7]
%       xs          - position of virtual source or direction of plane
%                     wave / m [1x3] or [1x6]
%       src         - source type of the virtual source
%                         'pw' - plane wave (xs is the direction of the
%                                plane wave in this case)
%                         'ps' - point source
%                         'ls' - line source
%                         'fs' - focused source
%       f           - frequency of the monochromatic source / Hz
%       conf        - configuration struct (see SFS_config)
%
%   Output parameters:
%       D           - driving function signal [nx1]
%
%   See also: plot_sound_field, sound_field_mono_wfs_25d,
%             driving_function_imp_wfs_25d

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input  parameters ==================================
nargmin = 5;
nargmax = 5;
narginchk(nargmin,nargmax);
isargsecondarysource(x0);
isargxs(xs);
isargpositivescalar(f);
isargchar(src);
isargstruct(conf);

%% ===== Computation ====================================================
% Calculate the driving function in time-frequency domain

% Secondary source positions and directions
nx0 = x0(:,4:6);
x0 = x0(:,1:3);

% Source position/direction/orientation
xs = repmat(xs,[size(x0,1) 1]);

% Get driving signals
if strcmp('pw',src)
    % === Plane wave ===
    % Direction of plane wave
    nk = bsxfun(@rdivide,xs,vector_norm(xs(:,1:3),2));
    % Driving signal
    D = driving_function_mono_wfs_pw(x0,nx0,nk,f,conf);

elseif strcmp('ps',src)
    % === Point source ===
    D = driving_function_mono_wfs_ps(x0,nx0,xs(:,1:3),f,conf);

elseif strcmp('ls',src)
    % === Line source ===
    D = driving_function_mono_wfs_ls(x0,nx0,xs,f,conf);

elseif strcmp('fs',src)
    % === Focused source ===
    D = driving_function_mono_wfs_fs(x0,nx0,xs(:,1:3),f,conf);

else
    error('%s: %s is not a known source type.',upper(mfilename),src);
end

author: sfstoolbox

repo: sfs-matlab


text:
function h = drawCameras(R, C, scale)
    h = hggroup;
    for i = 1:size(C,2)
        Ci = C(:,i);
        Ri = R(:,:,i);
        if ~isequal(Ri, zeros(3))
            drawCamera(Ci, Ri, scale, h);
            text(Ci(1), Ci(2), Ci(3), num2str(i-1), 'Parent', h);
        end
    end
end

function h = drawCamera(C, R, windowScale, handle)
h = hggroup;
% windowScale - scalar, specify the distance between optical center to
% camera frontal face; equal to half of frontal face side length

assert(isequal(size(C), [3,1]), 'size of C must be 3x1!');

window11 = windowScale*[1; 1; 1];
window12 = windowScale*[-1; 1; 1];
window21 = windowScale*[-1;-1;1];
window22 = windowScale*[1; -1; 1];
windowPrime11 = R'*window11+C;
windowPrime12 = R'*window12+C;
windowPrime21 = R'*window21+C;
windowPrime22 = R'*window22+C;
AxisEndPoints = C + windowScale * R'; % 3x3, each column is an endpoint coordinates of each axis
windowVertices = [windowPrime11 windowPrime12 windowPrime22 windowPrime21]; % 3x4, each column is a vertex coornidates of the camera window

% draw three camera axis
line([C(1) AxisEndPoints(1,1)], [C(2) AxisEndPoints(2,1)], [C(3) AxisEndPoints(3,1)], 'Color', 'r', 'LineWidth', 1, 'Parent', h);
line([C(1) AxisEndPoints(1,2)], [C(2) AxisEndPoints(2,2)], [C(3) AxisEndPoints(3,2)], 'Color', 'g', 'LineWidth', 1, 'Parent', h);
line([C(1) AxisEndPoints(1,3)], [C(2) AxisEndPoints(2,3)], [C(3) AxisEndPoints(3,3)], 'Color', 'b', 'LineWidth', 1, 'Parent', h);

% draw square window
line([windowPrime11(1), windowPrime12(1), windowPrime21(1), windowPrime22(1), windowPrime11(1)],...
    [windowPrime11(2), windowPrime12(2), windowPrime21(2), windowPrime22(2), windowPrime11(2)], ...
    [windowPrime11(3), windowPrime12(3), windowPrime21(3), windowPrime22(3), windowPrime11(3)], 'Color', [0.5 0.5 0.5], 'Parent', h);
% draw lines connecting camera center to four vertices of the square window
line([windowPrime11(1) C(1)], [windowPrime11(2) C(2)], [windowPrime11(3) C(3)], 'Color', [0.5 0.5 0.5], 'Parent', h);
line([windowPrime12(1) C(1)], [windowPrime12(2) C(2)], [windowPrime12(3) C(3)], 'Color', [0.5 0.5 0.5], 'Parent', h);
line([windowPrime22(1) C(1)], [windowPrime22(2) C(2)], [windowPrime22(3) C(3)], 'Color', [0.5 0.5 0.5], 'Parent', h);
line([windowPrime21(1) C(1)], [windowPrime21(2) C(2)], [windowPrime21(3) C(3)], 'Color', [0.5 0.5 0.5], 'Parent', h);

if exist('handle', 'var')
    h.Parent = handle;
end

end

author: zhixuany

repo: HUMBI

save_path: github-repos/MATLAB/zhixuany-HUMBI


text:
function img = pgmread(filename)
% function img = pgmread(filename)
%   this is my version of pgmread for the pgm file created by XV.
%
%  this program also corrects for the shifts in the image from pm file.

fid = fopen(filename,'r');
fscanf(fid, 'P5\n');
cmt = '#';
while findstr(cmt, '#'),
  cmt = fgets(fid);
   if length(findstr(cmt, '#')) ~= 1,
      YX = sscanf(cmt, '%d %d');
      y = YX(1); x = YX(2);
   end
end

%fgets(fid);

%img = fscanf(fid,'%d',size);
%img = img';

img = fread(fid,[y,x],'uint8');
img = img';
fclose(fid);

repo: ekfmonoslam


text:
function [ fea, out ] = ex_planestrain1( varargin )
%EX_PLANESTRAIN1 Plane strain analysis of a pressure vessel.
%
%   [ FEA, OUT ] = EX_PLANESTRAIN1( VARARGIN ) Benchmark example for plane strain
%   approximation of a pressure vessel (annular cross section with symmetry).
%
%   Reference.  B. J. Mac Donald, Practical Stress Analysis with Finite Elements (2nd Ed),
%   case study E on page 327, 2007.
%
%   Accepts the following property/value pairs.
%
%       Input       Value/{Default}        Description
%       -----------------------------------------------------------------------------------
%       E           scalar {207e9}         Modulus of elasticity
%       nu          scalar {0.27}          Poissons ratio
%       sfun        string {sflag1}        Shape function for displacements
%       iplot       scalar 0/{1}           Plot solution (=1)
%                                                                                .
%       Output      Value/(Size)           Description
%       -----------------------------------------------------------------------------------
%       fea         struct                 Problem definition struct
%       out         struct                 Output struct

% Copyright 2013-2022 Precise Simulation, Ltd.

cOptDef = { ...
  'E',        207e9; ...
  'nu',       0.27; ...
  'sfun',     'sflag1'; ...
  'iplot',    1; ...
  'igrid',    1; ...
  'tol',      0.1; ...
  'fid',      1 };
[got,opt] = parseopt( cOptDef, varargin{:} );
fid       = opt.fid;

% Geometry and grid.
fea.sdim = { 'x' 'y' };   % Coordinate names.
fea.grid = ringgrid( 12, 216, 100e-3, 120e-3 );
fea.grid = delcells( fea.grid, selcells( fea.grid, '(x<=eps) | (y<=eps)') );
if( opt.igrid~=1 )
  fea.grid = quad2tri( fea.grid );
end
n_bdr = max(fea.grid.b(3,:));   % Number of boundaries.

% Problem definition.
fea = addphys( fea, @planestrain );
fea.phys.psn.eqn.coef{1,end} = { opt.nu };
fea.phys.psn.eqn.coef{2,end} = { opt.E  };
fea.phys.psn.sfun            = { opt.sfun opt.sfun };

% Boundary conditions.
bctype = mat2cell( zeros(2,n_bdr), [1 1], ones(1,n_bdr) );
bctype{1,4} = 1;
bctype{2,3} = 1;
fea.phys.psn.bdr.coef{1,5} = bctype;

bccoef = mat2cell( zeros(2,n_bdr), [1 1], ones(1,n_bdr) );
bccoef{1,1} = '-nx*1e4';
bccoef{2,1} = '-ny*1e4';
fea.phys.psn.bdr.coef{1,end} = bccoef;

% Parse and solve problem.
fea       = parsephys( fea );
fea       = parseprob( fea );                          % Check and parse problem struct.
fea.sol.u = solvestat( fea, 'fid', fid, 'icub', 1+str2num(strrep(opt.sfun,'sflag','')) );   % Call to stationary solver.

% Postprocessing.
s_disp = fea.phys.psn.eqn.vars{2,end};
if( opt.iplot>0 )
  figure
  postplot( fea, 'surfexpr', s_disp )
  title( 'Total displacement' )
end

% Error checking.
s_sx     = fea.phys.psn.eqn.vars{5,end};
s_sy     = fea.phys.psn.eqn.vars{6,end};
s_sxy    = fea.phys.psn.eqn.vars{8,end};
s_sp1    = fea.phys.psn.eqn.vars{9,end};
s_sp2    = fea.phys.psn.eqn.vars{10,end};
s_sp3    = fea.phys.psn.eqn.vars{11,end};
v_disp   = evalexpr( s_disp, [100e-3 120e-3-2*sqrt(eps);0 0]+sqrt(eps), fea )';
v_dref   = [2.64e-8 2.41e-8];
[v_sx(1),v_sx(2)] = minmaxsubd( s_sx, fea );
v_sxref  = [-10000 55454];
[v_sy(1),v_sy(2)] = minmaxsubd( s_sy, fea );
v_syref  = [-10000 55454];
[v_sxy(1),v_sxy(2)] = minmaxsubd( s_sxy, fea );
v_sxyref = [-32730 0];
[v_sp1(1),v_sp1(2)] = minmaxsubd( s_sp1, fea );
v_sp1ref = [4.5e4 55454];
[v_sp2(1),v_sp2(2)] = minmaxsubd( s_sp2, fea );
v_sp2ref = [1.227e4 1.227e4];
[v_sp3(1),v_sp3(2)] = minmaxsubd( s_sp3, fea );
v_sp3ref = [-1e4 0];
out.err  = [ abs([v_dref-v_disp])./v_dref ;
             abs([v_sxref-v_sx])./v_sxref ;
             abs([v_syref-v_sy])./v_syref ;
             abs([v_sxyref(1)-v_sxy(1)])./v_sxyref(1) 0 ;
             abs([v_sp1ref(2)-v_sp1(2)])./v_sp1ref(2) 0 ;
             abs([v_sp2ref-v_sp2])./v_sp2ref      ;
             abs([v_sp3ref(1)-v_sp3(1)])./v_sp3ref(1) 0 ];
out.pass = all( out.err(:) <= opt.tol );

if( nargout==0 )
  clear fea out
end

author: precise-simulation

repo: featool-multiphysics


text:
%DETECTMARKERS  Basic ArUco marker detection
%
%     [corners, ids] = cv.detectMarkers(img, dictionary)
%     [corners, ids, rejectedImgPoints] = cv.detectMarkers(img, dictionary)
%     [...] = cv.detectMarkers(..., 'OptionName',optionValue, ...)
%
% ## Input
% * __img__ input image (8-bit grayscale or color).
% * __dictionary__ indicates the type of markers that will be searched. You
%   can specify the dictionary as a cell-array that starts with the type name
%   followed by option arguments `{Type, ...}`. There are three types of
%   dictionaries available:
%   * __Predefined__ `{'Predefined', name}` or simply as a string `name`.
%     Returns one of the predefined dictionaries.
%   * __Custom__ `{'Custom', nMarkers, markerSize, 'BaseDictionary',baseDict}`.
%     Generates a new customizable marker dictionary. This creates a new
%     dictionary composed by `nMarkers` markers and each markers composed by
%     `markerSize*markerSize` bits. If `BaseDictionary` is provided, its
%     markers are directly included and the rest are generated based on them.
%     If the size of `baseDict` is higher than `nMarkers`, only the first
%     `nMarkers` in `baseDict` are taken and no new marker is added.
%   * __Manual__ `{'Manual', bytesList, markerSize, maxCorrectionBits}`.
%     Creates a dictionary/set of markers manually. It contains the inner
%     codification.
%
% ## Output
% * __corners__ cell array of detected marker corners. For each marker, its
%   four corners are provided `{{[x1,y1],[x2,y2],[x3,y3],[x4,y4]}, ..}`. The
%   order of the corners is clockwise.
% * __ids__ vector of identifiers of the detected markers. The identifier is
%   of integer type (0-based). For N detected markers, the size of ids is also
%   N. The identifiers have the same order than the markers in the `corners`
%   array.
% * __rejectedImgPoints__ contains the `corners` of those squares whose inner
%   code has not a correct codification. Useful for debugging purposes.
%
% ## Options
% * __DetectorParameters__ marker detection parameters. A struct of parameters
%   for the detection process, with the following fields:
%   * __adaptiveThreshWinSizeMin__ minimum window size for adaptive
%     thresholding before finding contours (default 3).
%   * __adaptiveThreshWinSizeMax__ maximum window size for adaptive
%     thresholding before finding contours (default 23).
%   * __adaptiveThreshWinSizeStep__ increments from `AdaptiveThreshWinSizeMin`
%     to `AdaptiveThreshWinSizeMax` during the thresholding (default 10).
%   * __adaptiveThreshConstant__ constant for adaptive thresholding before
%     finding contours (default 7)
%   * __minMarkerPerimeterRate__ determine minimum perimeter for marker
%     contour to be detected. This is defined as a rate respect to the maximum
%     dimension of the input image (default 0.03).
%   * __maxMarkerPerimeterRate__ determine maximum perimeter for marker
%     contour to be detected. This is defined as a rate respect to the maximum
%     dimension of the input image (default 4.0).
%   * __polygonalApproxAccuracyRate__ minimum accuracy during the polygonal
%     approximation process to determine which contours are squares
%     (default 0.03).
%   * __minCornerDistanceRate__ minimum distance between corners for detected
%     markers relative to its perimeter (default 0.05)
%   * __minDistanceToBorder__ minimum distance of any corner to the image
%     border for detected markers (in pixels) (default 3)
%   * __minMarkerDistanceRate__ minimum mean distance beetween two marker
%     corners to be considered similar, so that the smaller one is removed.
%     The rate is relative to the smaller perimeter of the two markers
%     (default 0.05).
%   * __cornerRefinementMethod__ corner refinement method, one of:
%     * __None__ (default) no refinement.
%     * __Subpix__ do subpixel refinement (cv.cornerSubPix).
%     * __Contour__ refine the corners using the contour-points.
%   * __cornerRefinementWinSize__ window size for the corner refinement
%     process (in pixels) (default 5).
%   * __cornerRefinementMaxIterations__ maximum number of iterations for stop
%     criteria of the corner refinement process (default 30).
%   * __cornerRefinementMinAccuracy__ minimum error for the stop cristeria of
%     the corner refinement process (default 0.1)
%   * __markerBorderBits__ number of bits of the marker border, i.e. marker
%     border width (default 1).
%   * __perpectiveRemovePixelPerCell__ number of bits (per dimension) for each
%     cell of the marker when removing the perspective (default 4).
%   * __perspectiveRemoveIgnoredMarginPerCell__ width of the margin of pixels
%     on each cell not considered for the determination of the cell bit.
%     Represents the rate respect to the total size of the cell, i.e.
%     `PerpectiveRemovePixelPerCell` (default 0.13)
%   * __maxErroneousBitsInBorderRate__ maximum number of accepted erroneous
%     bits in the border (i.e. number of allowed white bits in the border).
%     Represented as a rate respect to the total number of bits per marker
%     (default 0.35).
%   * __minOtsuStdDev__ minimun standard deviation in pixels values during the
%     decodification step to apply Otsu thresholding (otherwise, all the bits
%     are set to 0 or 1 depending on mean higher than 128 or not) (default 5.0)
%   * __errorCorrectionRate__ error correction rate respect to the maximum
%     error correction capability for each dictionary. (default 0.6).
% * __CameraMatrix__ Optional 3x3 camera calibration matrix
%   `A = [fx 0 cx; 0 fy cy; 0 0 1]`.
% * __DistCoeffs__ Optional vector of camera distortion coefficients
%   `[k1,k2,p1,p2,k3,k4,k5,k6,s1,s2,s3,s4]` of 4, 5, 8 or 12 elements.
%
% ## Inputs for Predefined Dictionary
% * __name__ name of predefined markers dictionaries/sets. Each dictionary
%   indicates the number of bits and the number of markers contained:
%   * `4x4_50`, `4x4_100`, `4x4_250`, `4x4_1000`: 4x4 bits with
%     (50|100|250|1000) markers.
%   * `5x5_50`, `5x5_100`, `5x5_250`, `5x5_1000`: 5x5 bits with
%     (50|100|250|1000) markers.
%   * `6x6_50`, `6x6_100`, `6x6_250`, `6x6_1000`: 6x6 bits with
%     (50|100|250|1000) markers.
%   * `7x7_50`, `7x7_100`, `7x7_250`, `7x7_1000`: 7x7 bits with
%     (50|100|250|1000) markers.
%   * `ArucoOriginal`: standard ArUco Library Markers. 1024 markers, 5x5 bits,
%     0 minimum distance.
%
% ## Inputs for Custom Dictionary
% * __nMarkers__ number of markers in the dictionary.
% * __markerSize__ number of bits per dimension of each markers.
% * __baseDict__ (Optional) include the markers in this dictionary at the
%   beginning. It is specified in the same format as the parent dictionary
%   (i.e. a cell-array `{Type, ...}`).
%
% ## Inputs for Manual Dictionary
% * __bytesList__ marker code information. Note that the `bytesList` is given
%   in the form of a matrix of bits which is then converted to list of bytes
%   in the 4 rotations using `Dictionary::getByteListFromBits`. It contains
%   the marker codewords where:
%   * the number of rows is the dictionary size
%   * each marker is encoded using `nbytes = ceil(markerSize*markerSize/8)`
%   * each row contains all 4 rotations of the marker, so its length is
%     `4*nbytes`
%   * `bytesList(i,k*nbytes+j)` is then the j-th byte of i-th marker, in its
%     k-th rotation.
% * __markerSize__ number of bits per dimension.
% * __maxCorrectionBits__ maximum number of bits that can be corrected.
%
% The cv.detectMarkers function performs marker detection in the input image.
% Only markers included in the specific dictionary are searched. For each
% detected marker, it returns the 2D position of its corner in the image and
% its corresponding identifier. Note that this function does not perform pose
% estimation.
%
% # ArUco Marker Detection
%
% This module is dedicated to square fiducial markers (also known as Augmented
% Reality Markers). These markers are useful for easy, fast and robust camera
% pose estimation.
%
% The main functionalities are:
%
% - Detection of markers in a image
% - Pose estimation from a single marker or from a board/set of markers
% - Detection of ChArUco board for high subpixel accuracy
% - Camera calibration from both, ArUco boards and ChArUco boards.
% - Detection of ChArUco diamond markers
%
% The implementation is based on the ArUco Library [2] by [1].
%
% This module has been originally developed by Sergio Garrido-Jurado as a
% project for Google Summer of Code 2015 (GSoC 15).
%
% ### ArUco
%
% ArUco markers are easy to detect pattern grids that yield up to 1024
% different patterns. They were built for augmented reality and later used for
% camera calibration. Since the grid uniquely orients the square, the
% detection algorithm can determing the pose of the grid.
%
% ### ChArUco
%
% ArUco markers were improved by interspersing them inside a checkerboard
% called ChArUco. Checkerboard corner intersectionsa provide more stable
% corners because the edge location bias on one square is countered by the
% opposite edge orientation in the connecting square. By interspersing ArUco
% markers inside the checkerboard, each checkerboard corner gets a label which
% enables it to be used in complex calibration or pose scenarios where you
% cannot see all the corners of the checkerboard.
%
% The smallest ChArUco board is 5 checkers and 4 markers called a
% "Diamond Marker".
%
% ## References
% [1]:
% > S. Garrido-Jurado, R. Munoz-Salinas, F. J. Madrid-Cuevas, and
% > M. J. Marin-Jimenez. 2014. "Automatic generation and detection of highly
% > reliable fiducial markers under occlusion". Pattern Recogn. 47, 6
% > (June 2014), 2280-2292. DOI=10.1016/j.patcog.2014.01.005
%
% [2]:
% > ArUco: a minimal library for Augmented Reality applications
% > based on OpenCV, [ArUco](http://www.uco.es/investiga/grupos/ava/node/26)
%
% See also: cv.estimatePoseSingleMarkers, cv.estimatePoseBoard,
%  cv.refineDetectedMarkers, cv.drawDetectedMarkers, cv.adaptiveThreshold,
%  cv.threshold, cv.findContours, cv.cornerSubPix
%


text:
filename='CantileverBeam_Quadrilateral_Bilinear_Structured';
ptype = 'MACRO';
method = 'SIMPALL'; % !! Instead of proportional to material density !!
materialType = 'ISOTROPIC';
initial_case = 'holes';
cost = {'compliance','perimeter'};
weights = [1 0.1];
constraint = {'volumeConstraint'};
optimizer = 'HAMILTON-JACOBI';
incrementFactor = 1;
designVariable = 'LevelSet';
filterType = 'P1';
constraint_case = 'INEQUALITY';
line_search_initiator = 'INCREASING LAST STEP';

nsteps = 10;
Vfrac_final = 0.1;
Perimeter_target=3.5;
optimality_final =1e-3;
constr_final =1e-3;

Vfrac_initial = 0.4;
optimality_initial = 1e-3;
constr_initial = 1e-3;

TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 0.3;
TOL.nu_minus = 0.3;

plotting = 1;
printing = 0;
monitoring = 1;
monitoring_interval = 1;


text:
function hbg = bgCircleGraph(bg,NodeSizes,varargin)
% (optional) NodeSizes should be an array specifying the radius of each node
% as a fraction (0 to 1) of the radius of the circlegraph

if nargin<2
    NodeSizes = [];
elseif numel(NodeSizes)==1
    NodeSizes = NodeSizes*ones(1,length(bg.Nodes));
end

% I create a random graph with 28 nodes:
% bg = biograph(28);

% Display it in the GUI and get a handle
% back to work with it (this will use by default
% the hierarchical layout):
hbg = view(bg);
% Get an idea of what are the actual extents of
% the page in the GUI by looking at all the current
% node positions:
page_Size = max(cell2mat(arrayfun(@(x) get(x,'Position'),...
            get(hbg,'Nodes'),'Uniform',false)));
% Place the nodes in a circular layout, in my graph
% page_size was around [1000,1000], so I will select a
% center at [500 500] and a radius of 300, this will keep
% the layout approximately in the same scale and I will
% not have to manually change the size of the nodes or
% the fonts:
radius = 0.6*(min(page_Size)/2);
center = [page_Size(1)/2 page_Size(2)/2];
numNodes = length(hbg.Nodes);
for i = 1:numNodes
set(hbg.Nodes(i),'Position',...
   [center(1)+radius.*sin((i*2*pi/numNodes)),...
    center(2)+radius.*cos((i*2*pi/numNodes))],'shape','circle')

    if ~isempty(NodeSizes)
        sz = 2*radius*NodeSizes(i);
        set(hbg.Nodes(i),'Size',[sz sz]);
    end
end

% Now I can use dolayout:
set(bg,'NodeAutoSize','off','EdgeType','curved',varargin{:});
dolayout(hbg,'pathsOnly',true)


function feature_map = get_cnn_layers(im, fparams, gparams)

% Get layers from a cnn.

if size(im,3) == 1
    im = repmat(im, [1 1 3]);
end

im_sample_size = size(im);

%preprocess the image
if ~isequal(im_sample_size(1:2), fparams.net.meta.normalization.imageSize(1:2))
    im = imresize(single(im), fparams.net.meta.normalization.imageSize(1:2));
else
    im = single(im);
end

% Normalize with average image
im = bsxfun(@minus, im, fparams.net.meta.normalization.averageImage);

if fparams.use_gpu
    cnn_feat = vl_simplenn(fparams.net,im,[],[],'CuDNN',true, 'Mode', 'test');
else
    cnn_feat = vl_simplenn(fparams.net, im, [], [], 'Mode', 'test');
end

feature_map = cell(1,1,length(fparams.output_layer));

for k = 1:length(fparams.output_layer)
    % Move data to CPU if using GPU
%     if fparams.use_gpu
%         cnn_feat(fparams.output_layer(k) + 1).x = gather(cnn_feat(fparams.output_layer(k) + 1).x);
%     end

    if fparams.use_gpu
        if fparams.downsample_factor(k) == 1
            temp_feat = cnn_feat(fparams.output_layer(k) + 1).x(fparams.start_ind(k,1):fparams.end_ind(k,1), fparams.start_ind(k,2):fparams.end_ind(k,2), :, :);
            if fparams.return_gpu
                feature_map{k} = temp_feat;
            else
                feature_map{k} = gather(temp_feat);
            end
        else
            temp_feat = average_feature_region(cnn_feat(fparams.output_layer(k) + 1).x(fparams.start_ind(k,1):fparams.end_ind(k,1), fparams.start_ind(k,2):fparams.end_ind(k,2), :, :), fparams.downsample_factor(k));
            if fparams.return_gpu
                feature_map{k} = temp_feat;
            else
                feature_map{k} = gather(temp_feat);
            end
        end
    else
        if fparams.downsample_factor(k) == 1
            feature_map{k} = cnn_feat(fparams.output_layer(k) + 1).x(fparams.start_ind(k,1):fparams.end_ind(k,1), fparams.start_ind(k,2):fparams.end_ind(k,2), :, :);
        else
            feature_map{k} = average_feature_region(cnn_feat(fparams.output_layer(k) + 1).x(fparams.start_ind(k,1):fparams.end_ind(k,1), fparams.start_ind(k,2):fparams.end_ind(k,2), :, :), fparams.downsample_factor(k));
        end
    end
end


text:
function [Ka, Kb] = getKernels(ops, nup, sig)
% this function makes upsampling kernels for the temporal components.
% those are used for interpolating the biggest negative peak,
% and aligning the template to that peak with sub-sample resolution
% needs nup, the interpolation factor (default = 10)
% also needs sig, the interpolation smoothness (default = 1)

nt0min = getOr(ops, 'nt0min', 20);
nt0    = getOr(ops, 'nt0',    61);

xs = 1:nt0;
ys = linspace(.5, nt0+.5, nt0*nup+1);
ys(end) = [];

% these kernels are just standard kriging interpolators

% first compute distances between the sample coordinates
% for some reason, this seems to be circular, although the waveforms are not circular
% I think the reason had to do with some constant offsets in some channels?
d = rem(xs(:) - xs(:)' + nt0, nt0);
d = min(d, nt0-d);
Kxx = exp(-d.^2 / (sig^2)); % the kernel covariance uses a squared exponential of spatial scale sig

% do the same for the kernel similarities between upsampled "test" timepoints and the original coordinates
d = rem(ys(:) - xs(:)' + nt0, nt0);
d = min(d, nt0-d);
Kyx = exp(-d.^2 / (sig^2));

% the upsampling matrix is given by the following formula,
% with some light diagonal regularization of the matrix inversion
B = Kyx/(Kxx + .01 * eye(nt0));
B = reshape(B, nup, nt0, nt0);

% A is just a slice through this upsampling matrix corresponding to the most negative point
% this is used to compute the biggest negative deflection (after upsampling)
A = squeeze(B(:, nt0min, :));
B = permute(B, [2 3 1]);

% move to the GPU and make it a double
Ka = gpuArray(double(A));
Kb = gpuArray(double(B));

author: MouseLand

repo: Kilosort

save_path: github-repos/MATLAB/MouseLand-Kilosort


text:
function D = spm_eeg_ft_datareg_manual(varargin)
% Data registration user-interface routine
% commands the EEG/MEG data co-registration within original sMRI space
%
% FORMAT D = spm_eeg_inv_datareg_ui(D,[val], modality)
% Input:
% Output:
% D         - same data struct including the new required files and variables
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: spm_eeg_ft_datareg_manual.m 5397 2013-04-11 15:13:47Z vladimir $

% initialise
%--------------------------------------------------------------------------
[Finter, Fgraph] = spm('FnUIsetup','Fieldtrip MEEG/MRI manual coregistration', 0);

[D,val] = spm_eeg_inv_check(varargin{:});

[D, ok] = check(D, 'sensfid');

if ~ok
    if check(D, 'basic')
        error(['The requested file is not ready for source reconstruction.'...
            'Use prep to specify sensors and fiducials.']);
    else
        error('The meeg file is corrupt or incomplete');
    end
end

try
    D.inv{val}.mesh.template;
catch
    error('Please add a head model to the file before coregistering');
end

usepolhemus = spm_input('Use polhemus?', 1, 'yes|no', [1, 0]);

if usepolhemus
    meegfid = ft_read_headshape(spm_select(1, '\.*', 'Select polhemus file'));
    meegfid = ft_convert_units(meegfid, 'mm');
else
    meegfid = D.fiducials;
end

mrifid = D.inv{val}.mesh.fid;

if spm_input('Use individual head surface?', 1, 'yes|no', [1, 0])
    isurf = export(gifti(spm_select(1, 'mesh', 'Select head surface file')), 'ft');
    mrifid.pnt = isurf.pnt;
    mrifid.tri = isurf.tri;
end

meeglbl = meegfid.fid.label;
mrilbl = mrifid.fid.label;

newmrifid = mrifid;
newmrifid.fid.pnt = [];
newmrifid.fid.label = {};

if numel(meeglbl)> 3
    [selection ok]= listdlg('ListString', meeglbl, 'SelectionMode', 'multiple',...
        'InitialValue', spm_match_str(upper(meeglbl), upper(mrilbl)), ...
        'Name', 'Select at least 3 fiducials', 'ListSize', [400 300]);

    if ~ok || length(selection) < 3
        error('At least 3 M/EEG fiducials are required for coregistration');
    end

    meegfid.fid.pnt   = meegfid.fid.pnt(selection, :);
    meegfid.fid.label = meegfid.fid.label(selection);
    meeglbl = meeglbl(selection);
end

if numel(meeglbl)>=3
    for i = 1:length(meeglbl)
        switch spm_input(['How to specify ' meeglbl{i} ' position?'] , 1, 'select|type|click|skip')
            case 'select'
                [selection ok]= listdlg('ListString', mrilbl, 'SelectionMode', 'single',...
                    'InitialValue', strmatch(upper(meeglbl{i}), upper(mrilbl)), ...
                    'Name', ['Select matching MRI fiducial for ' meeglbl{i}], 'ListSize', [400 300]);
                if ~ok
                    continue
                end

                newmrifid.fid.pnt   = [newmrifid.fid.pnt; mrifid.fid.pnt(selection, :)];
            case 'type'
                pnt = spm_input('Input MNI coordinates', '+1', 'r', '', 3);
                newmrifid.fid.pnt   = [newmrifid.fid.pnt; pnt(:)'];
            case 'click'
                while 1
                    figure(Fgraph); clf;
                    mri = spm_vol(D.inv{val}.mesh.sMRI);
                    spm_orthviews('Reset');
                    spm_orthviews('Image', mri);
                    rotate3d off;
                    if spm_input(['Select ' meeglbl{i} ' position and click'] , 1,'OK|Retry', [1,0], 1)
                        newmrifid.fid.pnt   = [newmrifid.fid.pnt; spm_orthviews('Pos')'];
                        spm_orthviews('Reset');
                        break;
                    end
                end
            case 'skip'
                meegfid.fid.pnt(i, :) = [];
                meegfid.fid.label(i)  = [];
                continue;
        end
        newmrifid.fid.label = [newmrifid.fid.label  meeglbl{i}];
    end

    if size(newmrifid.fid.label) < 3
        error('At least 3 M/EEG fiducials are required for coregistration');
    end

    switch spm_input('Choose initial coregistration', 1, 'rigid|align1|align2|spm');
        case 'align1'
            M1 = spm_eeg_inv_headcoordinates(meegfid.fid.pnt(1, :), meegfid.fid.pnt(2, :), meegfid.fid.pnt(3, :));
            M =  spm_eeg_inv_headcoordinates(newmrifid.fid.pnt(1, :), newmrifid.fid.pnt(2, :), newmrifid.fid.pnt(3, :));
            M1 = inv(M) * M1;
        case 'align2'
            M1 = spm_eeg_inv_rigidreg(newmrifid.fid.pnt', meegfid.fid.pnt');
            tempfid = ft_transform_headshape(M1, meegfid);
            tempfid.fid.pnt(:, 2) = tempfid.fid.pnt(:, 2)- tempfid.fid.pnt(1, 2)+ newmrifid.fid.pnt(1, 2);
            tempfid.fid.pnt(:, 3) = tempfid.fid.pnt(:, 3)- mean(tempfid.fid.pnt(2:3, 3))+ mean(newmrifid.fid.pnt(2:3, 3));
            M1 = spm_eeg_inv_rigidreg(tempfid.fid.pnt', meegfid.fid.pnt');
        case 'rigid'
            M1 = spm_eeg_inv_rigidreg(newmrifid.fid.pnt', meegfid.fid.pnt');
        case 'spm'
            S =[];
            S.sourcefid = meegfid;
            S.targetfid = newmrifid;
            S.template  = 0;
            S.useheadshape = ~isempty(S.sourcefid.pnt);
            M1 = spm_eeg_inv_datareg(S);
    end
    meegfid = ft_transform_headshape(M1, meegfid);
end
%%
cfg = [];
cfg.individual.headshape = meegfid;
cfg.template.headshape = newmrifid;
cfg = ft_interactiverealign(cfg);

meegfid = ft_transform_headshape(cfg.m, meegfid);

M1 = cfg.m * M1;

if usepolhemus
    origfid = D.fiducials;
    sel1 = spm_match_str(origfid.fid.label, {'nas', 'lpa', 'rpa'});
    sel2 = spm_match_str(meegfid.fid.label, {'nas', 'lpa', 'rpa'});
    M2 = spm_eeg_inv_headcoordinates(meegfid.fid.pnt(sel2(1), :), meegfid.fid.pnt(sel2(2), :), meegfid.fid.pnt(sel2(3), :));
    M3 = spm_eeg_inv_headcoordinates(origfid.fid.pnt(sel1(1), :), origfid.fid.pnt(sel1(2), :), origfid.fid.pnt(sel1(3), :));
    M1 = inv(M2) * M3;
end

ind = 1;
D.inv{val}.datareg = struct([]);

if ~isempty(D.sensors('EEG'))
    D.inv{val}.datareg(ind).sensors = ft_transform_sens(M1, D.sensors('EEG'));
    D.inv{val}.datareg(ind).fid_eeg = ft_transform_headshape(M1, D.fiducials);
    D.inv{val}.datareg(ind).fid_mri = newmrifid;
    D.inv{val}.datareg(ind).toMNI = D.inv{val}.mesh.Affine;
    D.inv{val}.datareg(ind).fromMNI = inv(D.inv{val}.datareg(ind).toMNI);
    D.inv{val}.datareg(ind).modality = 'EEG';
    ind = ind+1;
end

if ~isempty(D.sensors('MEG'))
    D.inv{val}.datareg(ind).sensors = D.sensors('MEG');
    D.inv{val}.datareg(ind).fid_eeg = ft_transform_headshape(inv(M1), meegfid);
    D.inv{val}.datareg(ind).fid_mri = ft_transform_headshape(inv(M1), newmrifid);
    D.inv{val}.datareg(ind).toMNI = D.inv{val}.mesh.Affine*M1;
    D.inv{val}.datareg(ind).fromMNI = inv(D.inv{val}.datareg(ind).toMNI);
    D.inv{val}.datareg(ind).modality = 'MEG';
end

% check and display registration
%--------------------------------------------------------------------------
spm_eeg_inv_checkdatareg(D);

save(D);

author: spm


text:
filename='Cantilever_tetrahedra_coarse';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
optimizer = 'SLERP';
incrementFactor = 1;
designVariable = 'LevelSet';
filterType = 'P1';

nsteps = 5;
Vfrac_final = 0.05;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-2;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan


text:
function ans = eq_str (s1,s2);
%
%		ans = eq_str (s1,s2);
%
%	returns 1 if strings are equal
%          0 if they're not
%
%		Ziad Saad July 20 97

ls1 = length (s1);
ls2 = length (s2);

if (ls1 ~= ls2),
	ans = 0;
	return;
else
	ans = 1;
	for i=1:1:ls1,
		if (s1(i) ~= s2(i)),
			ans=0;
		end;
	end;
	return;
end;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [S,C] = conncomp(G)
  % CONNCOMP Drop in replacement for graphconncomp.m from the bioinformatics
  % toobox. G is an n by n adjacency matrix, then this identifies the S
  % connected components C. This is also an order of magnitude faster.
  %
  % [S,C] = conncomp(G)
  %
  % Inputs:
  %   G  n by n adjacency matrix, G(i,j) = G(j,i) ≠ 0 implies ij are connected.
  % Outputs:
  %   S  scalar number of connected components
  %   C

  % Transpose to match graphconncomp
  G = G';

  A = G+speye(size(G));
  [p,~,r] = dmperm(A);
  S = numel(r)-1;
  C = cumsum(full(sparse(1,r(1:end-1),1,1,size(G,1))));
  C(p) = C;
end

author: alecjacobson

repo: gptoolbox


text:
function [varargout]=subTriLoop(varargin)

% function [Fs,Vs,Cs,CV]=subTriLoop(F,V,n,fixBoundaryOpt)
% ------------------------------------------------------------------------
%
%
%
%
% ------------------------------------------------------------------------

%%

switch nargin
    case 2
        F=varargin{1};
        V=varargin{2};
        n=1;
        fixBoundaryOpt=0;
        logicFixedFaces=[];
    case 3
        F=varargin{1};
        V=varargin{2};
        n=varargin{3};
        fixBoundaryOpt=0;
        logicFixedFaces=[];
    case 4
        F=varargin{1};
        V=varargin{2};
        n=varargin{3};
        fixBoundaryOpt=varargin{4};
        logicFixedFaces=[];
    case 5
        F=varargin{1};
        V=varargin{2};
        n=varargin{3};
        fixBoundaryOpt=varargin{4};
        logicFixedFaces=varargin{5};
    case 6
        F=varargin{1};
        V=varargin{2};
        n=varargin{3};
        fixBoundaryOpt=varargin{4};
        logicFixedFaces=varargin{5};
end

%%

if nargout>2 || ~isempty(logicFixedFaces)
    C=(1:1:size(F,1))'; %Face colors or indices
    if nargout==4
        CV=zeros(size(V,1),1); %Vertex labels/colors
    end
end

if n>0
    for qIter=1:1:n

        M=patchConnectivity(F,V,{'ev','ef','vv'});

        edgeVertexMat=M.edge.vertex;
        edgeFaceMat=M.edge.face;
        vertexVertexMat=M.vertex.vertex;

        if ~isempty(logicFixedFaces)
            indFixFaces=find(logicFixedFaces);
            indFixEdges=find(any(ismember(edgeFaceMat,indFixFaces),2));
        else
            indFixEdges=[];
        end

        numPoints = size(V,1);
        numEdges = size(edgeVertexMat,1);

        % Get indices of the three edges associated with each face
        A = sparse(edgeVertexMat(:,1),edgeVertexMat(:,2),(1:numEdges)+numPoints,numPoints,numPoints,numEdges);
        A = max(A,A'); %Copy symmetric

        %Indices for A matrix
        indA_12=F(:,1)+(F(:,2)-1)*numPoints;
        indA_23=F(:,2)+(F(:,3)-1)*numPoints;
        indA_31=F(:,3)+(F(:,1)-1)*numPoints;

        %Get indices for vertex array
        indV_12=full(A(indA_12));
        indV_23=full(A(indA_23));
        indV_31=full(A(indA_31));

        %Create faces array
        Fs=[[F(:,1)  indV_12 indV_31];...
            [F(:,2)  indV_23 indV_12];...
            [F(:,3)  indV_31 indV_23];...
            [indV_12 indV_23 indV_31]];

        indF1=F(edgeFaceMat(:,1),:);
        if size(edgeFaceMat,2)==2
            logicValid=edgeFaceMat(:,2)>0;
            indF2=zeros(numEdges,3);
            indF2(logicValid,:)=F(edgeFaceMat(logicValid,2),:);
        else
            indF2=zeros(numEdges,3);
        end

        indF=[indF1 indF2];

        logicPick1=indF~=edgeVertexMat(:,ones(1,6)) & indF~=edgeVertexMat(:,2*ones(1,6));
        indF(~logicPick1)=0;
        indF=sort(indF,2,'descend');
        indF=indF(:,1:2);

        logicBoundaryEdges=indF(:,2)==0;

        N=sum(vertexVertexMat>0,2);
        beta = (1./N) .* (5/8 - ( 3/8 + (1/4)*(cos(2*pi./N) )).^2);

        %Computed sum of neighbours
        V_sum=zeros(size(V));
        logicValid=vertexVertexMat>0;
        X=nan(size(vertexVertexMat));
        for q=1:1:size(V,2)
            X(logicValid)=V(vertexVertexMat(logicValid),q);
            V_sum(:,q)=sum(X,2,'omitnan');
        end

        %Compute replacement input vertices
        Vv=(1-N.*beta).*V+beta.*V_sum;

        if any(logicBoundaryEdges)
            eb=edgeVertexMat(logicBoundaryEdges,:);
            indBoundaryVertices=unique(eb(:));

            EV=[eb;fliplr(eb)]; %Non-unique edges
            vertexVertexConnectivityBoundary=sparse(EV(:,1),EV(:,2),EV(:,2),size(V,1),size(V,1));
            vertexVertexConnectivityBoundary=sort(vertexVertexConnectivityBoundary,2,'descend');
            [~,J,~] = find(vertexVertexConnectivityBoundary);
            vertexVertexConnectivityBoundary=full(vertexVertexConnectivityBoundary(:,1:max(J)));

            Vv(indBoundaryVertices,:)= 6/8*V(indBoundaryVertices,:) + 1/8*(V(vertexVertexConnectivityBoundary(indBoundaryVertices,1),:)+V(vertexVertexConnectivityBoundary(indBoundaryVertices,2),:));
        end

        %Compute the new edge vertices
        Vn=zeros(numEdges,size(V,2));
        if fixBoundaryOpt==1
            %Use normal mid-edge nodes for boundary edges
            Vn(logicBoundaryEdges,:)=1/2*(V(edgeVertexMat(logicBoundaryEdges,1),:) + V(edgeVertexMat(logicBoundaryEdges,2),:));
            %Replace boundary nodes with original
            indBoundaryVertices=unique(edgeVertexMat(logicBoundaryEdges,:));
            Vv(indBoundaryVertices,:)=V(indBoundaryVertices,:);
        else
            Vn(logicBoundaryEdges,:)=1/2*V(edgeVertexMat(logicBoundaryEdges,1),:) + 1/2*V(edgeVertexMat(logicBoundaryEdges,2),:);
        end
        Vn(~logicBoundaryEdges,:)=3/8*V(edgeVertexMat(~logicBoundaryEdges,1),:) + 3/8*V(edgeVertexMat(~logicBoundaryEdges,2),:) + 1/8*(V(indF(~logicBoundaryEdges,1),:)+V(indF(~logicBoundaryEdges,2),:));

        if ~isempty(indFixEdges)
            %Use normal mid-edge nodes for constrained edges
            Vn(indFixEdges,:)=1/2*(V(edgeVertexMat(indFixEdges,1),:) + V(edgeVertexMat(indFixEdges,2),:));

            %Replace constrained nodes with original
            indConstrainedVertices=unique(edgeVertexMat(indFixEdges,:));
            Vv(indConstrainedVertices,:)=V(indConstrainedVertices,:);
        end

        %Join point sets
        Vs = [Vv; Vn];

        %Color handling
        if nargout>2 || ~isempty(logicFixedFaces)
            %Override face color data
            C=repmat(C,[size(Fs,1)/size(F,1),1]);

            %Update vertex labels
            if nargout==4
                if qIter==1
                    CV=[zeros(size(Vv,1),1); 1*ones(size(Vn,1),1); ];
                else
                    CV=[CV; qIter.*ones(size(Vn,1),1); ];
                end
            end
        end

        %Override face/vertices
        F=Fs;
        V=Vs;
        if ~isempty(logicFixedFaces)
            logicFixedFaces=logicFixedFaces(C);
        end
    end
end

%% Collect output

varargout{1}=F;
varargout{2}=V;
if nargout>2
    varargout{3}=C;
    if nargout==4
        varargout{4}=CV;
    end
end

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function f = f12_f0 ( n, x, y )

%*****************************************************************************80
%
%% F12_F0 returns the value of function 12.
%
%  Discussion:
%
%    This is an example from Vicente Romero.
%
%      R = sqrt ( X^2 + Y^2 )
%      T = atan ( Y / X )
%      F(X,Y) = ( 0.8 * R + 0.35 * sin ( 2.4 * pi * R / sqrt ( 2 )  ) )
%               * 1.5 * sin ( 1.3 * T )
%
%    The mean and standard deviation of the function over the interval
%    are approximately:
%
%      mu    = 0.581608
%      sigma = 0.343208
%
%    Since the interpolation interval is the unit square, this means the
%    integral of the function over the interval can also be estimated as
%
%      I = 0.581608
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 August 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Vicente Romero, John Burkardt, Max Gunzburger, Janet Peterson,
%    Initial Evaluation of Centroidal Voronoi Tessellation for
%    Statistical Sampling and Function Integration,
%    Fourth International Symposium on Uncertainty Modeling and Analysis,
%    (ISUMA) 2003.
%
%  Parameters:
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real X(N,1), Y(N,1), the evalution points.
%
%    Output, real F(N,1), the function values.
%
  r(1:n,1) = sqrt ( x(1:n,1).^2 + y(1:n,1).^2 );
  t(1:n,1) = atan2 ( y(1:n,1), x(1:n,1) );

  f(1:n,1) = 1.5 * ( 0.8 * r(1:n,1) ...
    + 0.35 * sin ( 2.4 * pi * r(1:n,1) / sqrt ( 2.0 ) ) ) ...
    .* sin ( 1.3 * t(1:n,1) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function G =  getGroupNonOverlapColor(row, col)

   % N = (row-2)*(col-2);
 N =  row * col ;
    g = sparse(zeros(N,1));
    g = diag(g);
    groupCount = 0;

    %% build non-overlapping group
    for i=2:3:col-1
        for j=2:3:row-1
            groupCount = groupCount+1;

            g(groupCount, (i-2)*row+j-1 ) = 1;
            g(groupCount, (i-2)*row+j   ) = 1;
            g(groupCount, (i-2)*row+j+1 ) = 1;
            g(groupCount, (i-1)*row+j-1 ) = 1;
            g(groupCount, (i-1)*row+j   ) = 1;
            g(groupCount, (i-1)*row+j+1 ) = 1;
            g(groupCount, i*row+j-1 ) = 1;
            g(groupCount, i*row+j   ) = 1;
            g(groupCount, i*row+j+1 ) = 1;

        end
    end

    g= g';

    g =g(:, 1:groupCount);

    G = sparse(3*N, 3*groupCount);

    G(1:N, 1:groupCount) = g;
    G(N+1:2*N, groupCount+1:2*groupCount) =g;
    G(2*N+1:3*N, 2*groupCount+1:3*groupCount) =g;

end

author: andrewssobral

repo: lrslibrary


text:
function value = year_length_gregorian ( y )

%*****************************************************************************80
%
%% YEAR_LENGTH_GREGORIAN returns the number of days in a Gregorian year.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 September 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, the year to be checked.
%
%    Output, integer VALUE, the number of
%    days in the year.
%
  if ( y == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'YEAR_LENGTH_GREGORIAN - Fatal error!\n' );
    fprintf ( 1, '  Illegal Y = 0.\n' );
    error ( 'YEAR_LENGTH_GREGORIAN - Fatal error!' );
  end

  if ( year_is_leap_gregorian ( y ) )
    value = 366;
  else
    value = 365;
  end

  return
end


function h = maxhorsign(signrestriction)

for ii = 1 : size(signrestriction,2)
    tmp    = (signrestriction{ii});
    index  = strfind(tmp,',');
    eval(['tmp1 = ' tmp( index(1)+1 : index(2)-1 ) ';']);
    hh(ii) = max(tmp1);
end

h = max(hh);


text:
%% AI Clinician Building MIMIC-III dataset

% (c) Matthieu Komorowski, Imperial College London 2015-2019
% as seen in publication: https://www.nature.com/articles/s41591-018-0213-5

% version 16 Feb 19
% uses the sepsis-3 cohort previously defined
% builds the MIMIC-III dataset

% This code is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE

%% ########################################################################
%           INITIAL REFORMAT WITH CHARTEVENTS, LABS AND MECHVENT
% ########################################################################

% gives an array with all unique charttime (1 per row) and all items in columns.
% ################## IMPORTANT !!!!!!!!!!!!!!!!!!
% Here i use -24 -> +48 because that's for the MDP

reformat=NaN(2000000,68);  %final table
qstime=zeros(100000,4);
winb4=25;   %lower limit for inclusion of data (48h before time flag)
winaft=49;  % upper limit (24h after)
irow=1;  %recording row for summary table
h = waitbar(0,'Initializing waitbar...');

tic
for icustayidrow=1:size(sepsis,1)

qst=sepsis.sepsis_time(icustayidrow);%,3); %flag for presumed infection
icustayid=sepsis.icustayid(icustayidrow)-200000;
waitbar(icustayidrow/size(sepsis,1),h,icustayidrow/size(sepsis,1)*100) %moved here to save some time

% CHARTEVENTS
    if icustayid<10000
    temp=ce010(ce010(:,1)==icustayid+200000,:);
    elseif icustayid>=10000 & icustayid<20000
    temp=ce1020(ce1020(:,1)==icustayid+200000,:);
    elseif icustayid>=20000 & icustayid<30000
    temp=ce2030(ce2030(:,1)==icustayid+200000,:);
    elseif icustayid>=30000 && icustayid<40000
    temp=ce3040(ce3040(:,1)==icustayid+200000,:);
    elseif icustayid>=40000 & icustayid<50000
    temp=ce4050(ce4050(:,1)==icustayid+200000,:);
    elseif icustayid>=50000 & icustayid<60000
    temp=ce5060(ce5060(:,1)==icustayid+200000,:);
    elseif icustayid>=60000 & icustayid<70000
    temp=ce6070(ce6070(:,1)==icustayid+200000,:);
    elseif icustayid>=70000 & icustayid<80000
    temp=ce7080(ce7080(:,1)==icustayid+200000,:);
    elseif icustayid>=80000 & icustayid<90000
    temp=ce8090(ce8090(:,1)==icustayid+200000,:);
    elseif icustayid>=90000
    temp=ce90100(ce90100(:,1)==icustayid+200000,:);
    end

ii=temp(:,2)>= qst-(winb4+4)*3600 & temp(:,2)<=qst+(winaft+4)*3600; %time period of interest -4h and +4h
temp=temp(ii,:);   %only time period of interest

%LABEVENTS
ii=labU(:,1)==icustayid+200000;
temp2=labU(ii,:);
ii=temp2(:,2)>= qst-(winb4+4)*3600 & temp2(:,2)<=qst+(winaft+4)*3600; %time period of interest -4h and +4h
temp2=temp2(ii,:);   %only time period of interest

%Mech Vent + ?extubated
ii=MV(:,1)==icustayid+200000;
temp3=MV(ii,:);
ii=temp3(:,2)>= qst-(winb4+4)*3600 & temp3(:,2)<=qst+(winaft+4)*3600; %time period of interest -4h and +4h
temp3=temp3(ii,:);   %only time period of interest

t=unique([temp(:,2);temp2(:,2); temp3(:,2)]);   %list of unique timestamps from all 3 sources / sorted in ascending order

if t
for i=1:numel(t)

    %CHARTEVENTS
    ii=temp(:,2)==t(i);
    col=temp(ii,3);
    value=temp(ii,4);
    reformat(irow,1)=i; %timestep
    reformat(irow,2)=icustayid;
    reformat(irow,3)=t(i); %charttime
    reformat(irow,3+col)=value; %store available values

    %LAB VALUES
    ii=temp2(:,2)==t(i);
    col=temp2(ii,3);
    value=temp2(ii,4);
    reformat(irow,31+col)=value; %store available values

    %MV
    ii=temp3(:,2)==t(i);
    if nansum(ii)>0
    value=temp3(ii,3:4);
      reformat(irow,67:68)=value; %store available values
    else
      reformat(irow,67:68)=NaN;
    end

    irow=irow+1;

end

qstime(icustayid,1)=qst; %time of sepsis
%HERE I SAVE FIRST and LAST TIMESTAMPS, in QSTIME, for each ICUSTAYID
qstime(icustayid,2)=t(1);  %first timestamp
qstime(icustayid,3)=t(end);  %last timestamp
qstime(icustayid,4)=table2array(demog(demog.icustay_id==icustayid+200000,5)); % discharge time

end

% end
end
toc

close(h);
reformat(irow:end,:)=[];  %delete extra unused rows

%% ########################################################################
%                                   OUTLIERS
% ########################################################################

%weight
reformat=deloutabove(reformat,5,300);  %delete outlier above a threshold (300 kg), for variable # 5

%HR
reformat=deloutabove(reformat,8,250);

%BP
reformat=deloutabove(reformat,9,300);
reformat=deloutbelow(reformat,10,0);
reformat=deloutabove(reformat,10,200);
reformat=deloutbelow(reformat,11,0);
reformat=deloutabove(reformat,11,200);

%RR
reformat=deloutabove(reformat,12,80);

%SpO2
reformat=deloutabove(reformat,13,150);
ii=reformat(:,13)>100;reformat(ii,13)=100;
reformat=deloutbelow(reformat,13,50);

%temp
ii=reformat(:,14)>90 & isnan(reformat(:,15));reformat(ii,15)=reformat(ii,14);
reformat=deloutabove(reformat,14,90);
reformat=deloutbelow(reformat,14,25);

%interface / is in col 22

% FiO2
reformat=deloutabove(reformat,23,100);
ii=reformat(:,23)<1;reformat(ii,23)=reformat(ii,23)*100;
reformat=deloutbelow(reformat,23,20);
reformat=deloutabove(reformat,24,1.5);

% O2 FLOW
reformat=deloutabove(reformat,25,70);

%PEEP
reformat=deloutbelow(reformat,26,0);
reformat=deloutabove(reformat,26,40);

%TV
reformat=deloutabove(reformat,27,1800);

%MV
reformat=deloutabove(reformat,28,50);

%K+
reformat=deloutbelow(reformat,32,1);
reformat=deloutabove(reformat,32,15);

%Na
reformat=deloutbelow(reformat,33,95);
reformat=deloutabove(reformat,33,178);

%Cl
reformat=deloutbelow(reformat,34,70);
reformat=deloutabove(reformat,34,150);

%Glc
reformat=deloutbelow(reformat,35,1);
reformat=deloutabove(reformat,35,1000);

%Creat
reformat=deloutabove(reformat,37,150);

%Mg
reformat=deloutabove(reformat,38,10);

%Ca
reformat=deloutabove(reformat,39,20);

%ionized Ca
reformat=deloutabove(reformat,40,5);

%CO2
reformat=deloutabove(reformat,41,120);

%SGPT/SGOT
reformat=deloutabove(reformat,42,10000);
reformat=deloutabove(reformat,43,10000);

%Hb/Ht
reformat=deloutabove(reformat,50,20);
reformat=deloutabove(reformat,51,65);

%WBC
reformat=deloutabove(reformat,53,500);

%plt
reformat=deloutabove(reformat,54,2000);

%INR
reformat=deloutabove(reformat,58,20);

%pH
reformat=deloutbelow(reformat,59,6.7);
reformat=deloutabove(reformat,59,8);

%po2
reformat=deloutabove(reformat,60,700);

%pco2
reformat=deloutabove(reformat,61,200);

%BE
reformat=deloutbelow(reformat,62,-50);

%lactate
reformat=deloutabove(reformat,63,30);

% ####################################################################
% some more data manip / imputation from existing values

% estimate GCS from RASS - data from Wesley JAMA 2003
ii=isnan(reformat(:,6))&reformat(:,7)>=0;
reformat(ii,6)=15;
ii=isnan(reformat(:,6))&reformat(:,7)==-1;
reformat(ii,6)=14;
ii=isnan(reformat(:,6))&reformat(:,7)==-2;
reformat(ii,6)=12;
ii=isnan(reformat(:,6))&reformat(:,7)==-3;
reformat(ii,6)=11;
ii=isnan(reformat(:,6))&reformat(:,7)==-4;
reformat(ii,6)=6;
ii=isnan(reformat(:,6))&reformat(:,7)==-5;
reformat(ii,6)=3;

% FiO2
ii=~isnan(reformat(:,23)) & isnan(reformat(:,24));
reformat(ii,24)=reformat(ii,23)./100;
ii=~isnan(reformat(:,24)) & isnan(reformat(:,23));
reformat(ii,23)=reformat(ii,24).*100;

%ESTIMATE FiO2 /// with use of interface / device (cannula, mask, ventilator....)

reformatsah=SAH(reformat,sample_and_hold);  % do SAH first to handle this task

%NO FiO2, YES O2 flow, no interface OR cannula
ii=find(isnan(reformatsah(:,23))&~isnan(reformatsah(:,25))&(reformatsah(:,22)==0|reformatsah(:,22)==2));
reformat(ii(reformatsah(ii,25)<=15),23)=70;
reformat(ii(reformatsah(ii,25)<=12),23)=62;
reformat(ii(reformatsah(ii,25)<=10),23)=55;
reformat(ii(reformatsah(ii,25)<=8),23)=50;
reformat(ii(reformatsah(ii,25)<=6),23)=44;
reformat(ii(reformatsah(ii,25)<=5),23)=40;
reformat(ii(reformatsah(ii,25)<=4),23)=36;
reformat(ii(reformatsah(ii,25)<=3),23)=32;
reformat(ii(reformatsah(ii,25)<=2),23)=28;
reformat(ii(reformatsah(ii,25)<=1),23)=24;

%NO FiO2, NO O2 flow, no interface OR cannula
ii=find(isnan(reformatsah(:,23))&isnan(reformatsah(:,25))&(reformatsah(:,22)==0|reformatsah(:,22)==2));  %no fio2 given and o2flow given, no interface OR cannula
reformat(ii,23)=21;

%NO FiO2, YES O2 flow, face mask OR.... OR ventilator (assume it's face mask)
ii=find(isnan(reformatsah(:,23))&~isnan(reformatsah(:,25))&(reformatsah(:,22)==NaN|reformatsah(:,22)==1|reformatsah(:,22)==3|reformatsah(:,22)==4|reformatsah(:,22)==5|reformatsah(:,22)==6|reformatsah(:,22)==9|reformatsah(:,22)==10));
reformat(ii(reformatsah(ii,25)<=15),23)=75;
reformat(ii(reformatsah(ii,25)<=12),23)=69;
reformat(ii(reformatsah(ii,25)<=10),23)=66;
reformat(ii(reformatsah(ii,25)<=8),23)=58;
reformat(ii(reformatsah(ii,25)<=6),23)=40;
reformat(ii(reformatsah(ii,25)<=4),23)=36;

%NO FiO2, NO O2 flow, face mask OR ....OR ventilator
ii=find(isnan(reformatsah(:,23))&isnan(reformatsah(:,25))&(reformatsah(:,22)==NaN|reformatsah(:,22)==1|reformatsah(:,22)==3|reformatsah(:,22)==4|reformatsah(:,22)==5|reformatsah(:,22)==6|reformatsah(:,22)==9|reformatsah(:,22)==10));  %no fio2 given and o2flow given, no interface OR cannula
reformat(ii,23)=NaN;

%NO FiO2, YES O2 flow, Non rebreather mask
ii=find(isnan(reformatsah(:,23))&~isnan(reformatsah(:,25))&reformatsah(:,22)==7);
reformat(ii(reformatsah(ii,25)>=10),23)=90;
reformat(ii(reformatsah(ii,25)>=15),23)=100;
reformat(ii(reformatsah(ii,25)<10),23)=80;
reformat(ii(reformatsah(ii,25)<=8),23)=70;
reformat(ii(reformatsah(ii,25)<=6),23)=60;

%NO FiO2, NO O2 flow, NRM
ii=find(isnan(reformatsah(:,23))&isnan(reformatsah(:,25))&reformatsah(:,22)==7);  %no fio2 given and o2flow given, no interface OR cannula
reformat(ii,23)=NaN;

% update again FiO2 columns
ii=~isnan(reformat(:,23)) & isnan(reformat(:,24));
reformat(ii,24)=reformat(ii,23)./100;
ii=~isnan(reformat(:,24)) & isnan(reformat(:,23));
reformat(ii,23)=reformat(ii,24).*100;

%BP
ii=~isnan(reformat(:,9))&~isnan(reformat(:,10)) & isnan(reformat(:,11));
reformat(ii,11)=(3*reformat(ii,10)-reformat(ii,9))./2;
ii=~isnan(reformat(:,09))&~isnan(reformat(:,11)) & isnan(reformat(:,10));
reformat(ii,10)=(reformat(ii,9)+2*reformat(ii,11))./3;
ii=~isnan(reformat(:,10))&~isnan(reformat(:,11)) & isnan(reformat(:,9));
reformat(ii,9)=3*reformat(ii,10)-2*reformat(ii,11);

%TEMP
%some values recorded in the wrong column
ii=reformat(:,15)>25&reformat(:,15)<45; %tempF close to 37deg??!
reformat(ii,14)=reformat(ii,15);
reformat(ii,15)=NaN;
ii=reformat(:,14)>70;  %tempC > 70?!!! probably degF
reformat(ii,15)=reformat(ii,14);
reformat(ii,14)=NaN;
ii=~isnan(reformat(:,14)) & isnan(reformat(:,15));
reformat(ii,15)=reformat(ii,14)*1.8+32;
ii=~isnan(reformat(:,15)) & isnan(reformat(:,14));
reformat(ii,14)=(reformat(ii,15)-32)./1.8;

% Hb/Ht
ii=~isnan(reformat(:,50)) & isnan(reformat(:,51));
reformat(ii,51)=(reformat(ii,50)*2.862)+1.216;
ii=~isnan(reformat(:,51)) & isnan(reformat(:,50));
reformat(ii,50)=(reformat(ii,51)-1.216)./2.862;

%BILI
ii=~isnan(reformat(:,44)) & isnan(reformat(:,45));
reformat(ii,45)=(reformat(ii,44)*0.6934)-0.1752;
ii=~isnan(reformat(:,45)) & isnan(reformat(:,44));
reformat(ii,44)=(reformat(ii,45)+0.1752)./0.6934;

%% ########################################################################
%                      SAMPLE AND HOLD on RAW DATA
% ########################################################################

reformat=SAH(reformat(:,1:68),sample_and_hold);

%% ########################################################################
%                             DATA COMBINATION
% ########################################################################

tic
     save('D:\BACKUP MIT PC\Data_100219.mat', '-v7.3');
toc

% WARNING: the time window of interest has been defined above (here -48 -> +24)!

timestep=4;  %resolution of timesteps, in hours
irow=1;
icustayidlist=unique(reformat(:,2));
reformat2=nan(size(reformat,1),84);  %output array
h = waitbar(0,'Initializing waitbar...');
npt=numel(icustayidlist);  %number of patients
% Adding 2 empty cols for future shock index=HR/SBP and P/F
reformat(:,69:70)=NaN(size(reformat,1),2);

tic
for i=1:npt

    icustayid=icustayidlist(i);  %1 to 100000, NOT 200 to 300K!

        %CHARTEVENTS AND LAB VALUES
        temp=reformat(reformat(:,2)==icustayid,:);   %subtable of interest
        beg=temp(1,3);   %timestamp of first record

        % IV FLUID STUFF
        iv=find(inputMV(:,1)==icustayid+200000);   %rows of interest in inputMV
        input=inputMV(iv,:);    %subset of interest
        iv=find(inputCV(:,1)==icustayid+200000);   %rows of interest in inputCV
        input2=inputCV(iv,:);    %subset of interest
        startt=input(:,2); %start of all infusions and boluses
        endt=input(:,3); %end of all infusions and boluses
        rate=input(:,8);  %rate of infusion (is NaN for boluses)

        pread=inputpreadm(inputpreadm(:,1)==icustayid+200000,2) ;%preadmission volume
            if ~isempty(pread)             %store the value, if available
                totvol=nansum(pread);
                waitbar(i/npt,h,i/npt*100) %moved here to save some time
            else
                totvol=0;   %if not documented: it's zero
            end

        % compute volume of fluid given before start of record!!!
        t0=0;
        t1=beg;
        %input from MV (4 ways to compute)
        infu=  nansum(rate.*(endt-startt).*(endt<=t1&startt>=t0)/3600   +    rate.*(endt-t0).*(startt<=t0&endt<=t1&endt>=t0)/3600 +     rate.*(t1-startt).*(startt>=t0&endt>=t1&startt<=t1)/3600 +      rate.*(t1-t0).*(endt>=t1&startt<=t0)   /3600);
        %all boluses received during this timestep, from inputMV (need to check rate is NaN) and inputCV (simpler):
        bolus=nansum(input(isnan(input(:,6))& input(:,2)>=t0&input(:,2)<=t1,7)) + nansum(input2(input2(:,2)>=t0&input2(:,2)<=t1,5));
        totvol=nansum([totvol,infu,bolus]);

        %VASOPRESSORS
        iv=find(vasoMV(:,1)==icustayid+200000);   %rows of interest in vasoMV
        vaso1=vasoMV(iv,:);    %subset of interest
        iv=find(vasoCV(:,1)==icustayid+200000);   %rows of interest in vasoCV
        vaso2=vasoCV(iv,:);    %subset of interest
        startv=vaso1(:,3); %start of VP infusion
        endv=vaso1(:,4); %end of VP infusions
        ratev=vaso1(:,5);  %rate of VP infusion

        %DEMOGRAPHICS / gender, age, elixhauser, re-admit, died in hosp?, died within
        %48h of out_time (likely in ICU or soon after), died within 90d after admission?
        demogi=find(demog.icustay_id==icustayid+200000);
        dem=[  demog.gender(demogi) ; demog.age(demogi) ;demog.elixhauser(demogi) ; demog.adm_order(demogi)>1 ;  demog.morta_hosp(demogi); abs(demog.dod(demogi)-demog.outtime(demogi))<(24*3600*2); demog.morta_90(demogi) ; (qstime(icustayid,4)-qstime(icustayid,3))/3600];

        % URINE OUTPUT
        iu=find(UO(:,1)==icustayid+200000);   %rows of interest in inputMV
        output=UO(iu,:);    %subset of interest
        pread=UOpreadm(UOpreadm(:,1)==icustayid,4) ;%preadmission UO
            if ~isempty(pread)     %store the value, if available
                UOtot=nansum(pread);
            else
                UOtot=0;
            end
        % adding the volume of urine produced before start of recording!
        UOnow=nansum(output(output(:,2)>=t0&output(:,2)<=t1,4));  %t0 and t1 defined above
        UOtot=nansum([UOtot UOnow]);

    for j=0:timestep:79 % -28 until +52 = 80 hours in total
        t0=3600*j+ beg;   %left limit of time window
        t1=3600*(j+timestep)+beg;   %right limit of time window
        ii=temp(:,3)>=t0 & temp(:,3)<=t1;  %index of items in this time period
        if sum(ii)>0

        %ICUSTAY_ID, OUTCOMES, DEMOGRAPHICS
        reformat2(irow,1)=(j/timestep)+1;   %'bloc' = timestep (1,2,3...)
        reformat2(irow,2)=icustayid;        %icustay_ID
        reformat2(irow,3)=3600*j+ beg;      %t0 = lower limit of time window
        reformat2(irow,4:11)=dem;           %demographics and outcomes

        %CHARTEVENTS and LAB VALUES (+ includes empty cols for shock index and P/F)
        value=temp(ii,:);%records all values in this timestep

        if sum(ii)==1   %if only 1 row of values at this timestep
          reformat2(irow,12:78)=value(:,4:end);
        else
          reformat2(irow,12:78)=nanmean(value(:,4:end)); %mean of all available values

        end

        %VASOPRESSORS
            % for CV: dose at timestamps.
            % for MV: 4 possibles cases, each one needing a different way to compute the dose of VP actually administered:
            %----t0---start----end-----t1----
            %----start---t0----end----t1----
            %-----t0---start---t1---end
            %----start---t0----t1---end----

        %MV
        v=(endv>=t0&endv<=t1)|(startv>=t0&endv<=t1)|(startv>=t0&startv<=t1)|(startv<=t0&endv>=t1);
        %CV
        v2=vaso2(vaso2(:,3)>=t0&vaso2(:,3)<=t1,4);
        v1=nanmedian([ratev(v); v2]);
        v2=nanmax([ratev(v); v2]);
        if ~isempty(v1)&~isnan(v1)&~isempty(v2)&~isnan(v2)
        reformat2(irow,79)=v1;    %median of dose of VP
        reformat2(irow,80)=v2;    %max dose of VP
        end

        %INPUT FLUID
        %input from MV (4 ways to compute)
        infu=  nansum(rate.*(endt-startt).*(endt<=t1&startt>=t0)/3600   +    rate.*(endt-t0).*(startt<=t0&endt<=t1&endt>=t0)/3600 +     rate.*(t1-startt).*(startt>=t0&endt>=t1&startt<=t1)/3600 +      rate.*(t1-t0).*(endt>=t1&startt<=t0)   /3600);
        %all boluses received during this timestep, from inputMV (need to check rate is NaN) and inputCV (simpler):
        bolus=nansum(input(isnan(input(:,6))& input(:,2)>=t0&input(:,2)<=t1,7)) + nansum(input2(input2(:,2)>=t0&input2(:,2)<=t1,5));

        %sum fluid given
        totvol=nansum([totvol,infu,bolus]);
        reformat2(irow,81)=totvol;    %total fluid given
        reformat2(irow,82)=nansum([infu,bolus]);   %fluid given at this step

        %UO
        UOnow=nansum(output(output(:,2)>=t0&output(:,2)<=t1,4));
        UOtot=nansum([UOtot UOnow]);
        reformat2(irow,83)=UOtot;    %total UO
        reformat2(irow,84)=nansum(UOnow);   %UO at this step

        %CUMULATED BALANCE
        reformat2(irow,85)=totvol-UOtot;    %cumulated balance

        irow=irow+1;
        end
    end
end
toc

reformat2(irow:end,:)=[];
close(h);

tic
     save('D:\BACKUP MIT PC\Data_110219.mat', '-v7.3');
toc

%% ########################################################################
%             CONVERT TO TABLE AND KEEP ONLY WANTED VARIABLE
% ########################################################################

dataheaders=[sample_and_hold(1,:) {'Shock_Index' 'PaO2_FiO2'}];
dataheaders=regexprep(dataheaders,'['']','');
dataheaders = ['bloc','icustayid','charttime','gender','age','elixhauser','re_admission', 'died_in_hosp', 'died_within_48h_of_out_time','mortality_90d','delay_end_of_record_and_discharge_or_death',...
    dataheaders,  'median_dose_vaso','max_dose_vaso','input_total','input_4hourly','output_total','output_4hourly','cumulated_balance'];

reformat2t=array2table(reformat2);
reformat2t.Properties.VariableNames=dataheaders;

% headers I want to keep
dataheaders5 = {'bloc','icustayid','charttime','gender','age','elixhauser','re_admission', 'died_in_hosp', 'died_within_48h_of_out_time','mortality_90d','delay_end_of_record_and_discharge_or_death','SOFA','SIRS',...
    'Weight_kg','GCS','HR','SysBP','MeanBP','DiaBP','RR','SpO2','Temp_C','FiO2_1','Potassium','Sodium','Chloride','Glucose',...
    'BUN','Creatinine','Magnesium','Calcium','Ionised_Ca','CO2_mEqL','SGOT','SGPT','Total_bili','Albumin','Hb','WBC_count','Platelets_count','PTT','PT','INR',...
    'Arterial_pH','paO2','paCO2','Arterial_BE','HCO3','Arterial_lactate','mechvent','Shock_Index','PaO2_FiO2',...
    'median_dose_vaso','max_dose_vaso','input_total','input_4hourly','output_total','output_4hourly','cumulated_balance'};

ii=find(ismember(reformat2t.Properties.VariableNames,dataheaders5));
reformat3t=reformat2t(:,ii);

%% SOME DATA MANIP BEFORE IMPUTATION

% CORRECT GENDER
reformat3t.gender=reformat3t.gender-1;

%CORRECT AGE > 200 yo
ii=reformat3t.age>150*365.25;
reformat3t.age(ii)=91.4*365.25;

% FIX MECHVENT
reformat3t.mechvent(isnan(reformat3t.mechvent))=0;
reformat3t.mechvent(reformat3t.mechvent>0)=1;

% FIX Elixhauser missing values
reformat3t.elixhauser(isnan(reformat3t.elixhauser))=nanmedian(reformat3t.elixhauser);  %use the median value / only a few missing data points

%vasopressors / no NAN
a=find(ismember(reformat3t.Properties.VariableNames,{'median_dose_vaso'}));
ii=isnan(table2array(reformat3t(:,a)));
reformat3t(ii,a)=array2table(zeros(sum(ii),1));
a=find(ismember(reformat3t.Properties.VariableNames,{'max_dose_vaso'}));
ii=isnan(table2array(reformat3t(:,a)));
reformat3t(ii,a)=array2table(zeros(sum(ii),1));

% check prop of missingness here
miss=array2table(sum(isnan(table2array(reformat3t)))./size(reformat3t,1));
miss.Properties.VariableNames=reformat3t.Properties.VariableNames;
figure; bar(sum(isnan(table2array(reformat3t)))./size(reformat3t,1))

% I fill the values temporarily with zeros, otherwise kNN imp doesnt work
reformat3t.Shock_Index=zeros(size(reformat3t,1),1);
reformat3t.PaO2_FiO2=zeros(size(reformat3t,1),1);

%% ########################################################################
%        HANDLING OF MISSING VALUES & CREATE REFORMAT4T
% ########################################################################

% Do linear interpol where missingness is low (kNN imputation doesnt work if all rows have missing values)
reformat3=table2array(reformat3t);
miss=sum(isnan((reformat3)))./size(reformat3,1);
ii=miss>0&miss<0.05;  %less than 5% missingness
mechventcol=find(ismember(reformat3t.Properties.VariableNames,{'mechvent'}));

for i=11:mechventcol-1 % correct col by col, otherwise it does it wrongly
  if ii(i)==1
    reformat3(:,i)=fixgaps(reformat3(:,i));
  end
end

reformat3t(:,11:mechventcol-1)=array2table(reformat3(:,11:mechventcol-1));

% KNN IMPUTATION -  Done on chunks of 10K records.

mechventcol=find(ismember(reformat3t.Properties.VariableNames,{'mechvent'}));
ref=reformat3(:,11:mechventcol-1);  %columns of interest

tic
for i=1:10000:size(reformat3,1)-9999   %dataset divided in 5K rows chunks (otherwise too large)
    i
    ref(i:i+9999,:)=knnimpute(ref(i:i+9999,:)',1, 'distance','seuclidean')';
end

ref(end-9999:end,:)=knnimpute(ref(end-9999:end,:)',1, 'distance','seuclidean')';  %the last bit is imputed from the last 10K rows

toc

% I paste the data interpolated, but not the demographics and the treatments
reformat3t(:,11:mechventcol-1)=array2table(ref);

reformat4t=reformat3t;
reformat4=table2array(reformat4t);

%% ########################################################################
%        COMPUTE SOME DERIVED VARIABLES: P/F, Shock Index, SOFA, SIRS...
% ########################################################################

% re-compute P/F with no missing values...
p=find(ismember(reformat4t.Properties.VariableNames,{'paO2'}));
f=find(ismember(reformat4t.Properties.VariableNames,{'FiO2_1'}));
a=find(ismember(reformat4t.Properties.VariableNames,{'PaO2_FiO2'}));
reformat4t(:,a)=array2table(reformat4(:,p)./reformat4(:,f));

%recompute SHOCK INDEX without NAN and INF
p=find(ismember(reformat4t.Properties.VariableNames,{'HR'}));
f=find(ismember(reformat4t.Properties.VariableNames,{'SysBP'}));
a=find(ismember(reformat4t.Properties.VariableNames,{'Shock_Index'}));
reformat4(:,a)=reformat4(:,p)./reformat4(:,f);
reformat4(isinf(reformat4(:,a)),a)=NaN;
d=nanmean(reformat4(:,a));
reformat4(isnan(reformat4(:,a)),a)=d;  %replace NaN with average value ~ 0.8
reformat4t(:,a)=array2table(reformat4(:,a));
ii=reformat4t.Shock_Index>=quantile(reformat4t.Shock_Index,0.999); %replace outliers with 99.9th percentile
reformat4t.Shock_Index(ii)=quantile(reformat4t.Shock_Index,0.999);

% SOFA - at each timepoint
% need (in this order):  P/F  MV  PLT  TOT_BILI  MAP  NORAD(max)  GCS  CR  UO
a=zeros(8,1);                              % indices of vars used in SOFA
a(1)=find(ismember(reformat4t.Properties.VariableNames,{'PaO2_FiO2'}));
a(2)=find(ismember(reformat4t.Properties.VariableNames,{'Platelets_count'}));
a(3)=find(ismember(reformat4t.Properties.VariableNames,{'Total_bili'}));
a(4)=find(ismember(reformat4t.Properties.VariableNames,{'MeanBP'}));
a(5)=find(ismember(reformat4t.Properties.VariableNames,{'max_dose_vaso'}));
a(6)=find(ismember(reformat4t.Properties.VariableNames,{'GCS'}));
a(7)=find(ismember(reformat4t.Properties.VariableNames,{'Creatinine'}));
a(8)=find(ismember(reformat4t.Properties.VariableNames,{'output_4hourly'}));
s=table2array(reformat4t(:,a));

p=[0 1 2 3 4];

s1=[s(:,1)>400 s(:,1)>=300 &s(:,1)<400 s(:,1)>=200 &s(:,1)<300 s(:,1)>=100 &s(:,1)<200 s(:,1)<100 ];   %count of points for all 6 criteria of sofa
s2=[s(:,2)>150 s(:,2)>=100 &s(:,2)<150 s(:,2)>=50 &s(:,2)<100 s(:,2)>=20 &s(:,2)<50 s(:,2)<20 ];
s3=[s(:,3)<1.2 s(:,3)>=1.2 &s(:,3)<2 s(:,3)>=2 &s(:,3)<6 s(:,3)>=6 &s(:,3)<12 s(:,3)>12 ];
s4=[s(:,4)>=70 s(:,4)<70&s(:,4)>=65 s(:,4)<65 s(:,5)>0 &s(:,5)<=0.1 s(:,5)>0.1 ];
s5=[s(:,6)>14 s(:,6)>12 &s(:,6)<=14 s(:,6)>9 &s(:,6)<=12 s(:,6)>5 &s(:,6)<=9 s(:,6)<=5 ];
s6=[s(:,7)<1.2 s(:,7)>=1.2 &s(:,7)<2 s(:,7)>=2 &s(:,7)<3.5 (s(:,7)>=3.5 &s(:,7)<5)|(s(:,8)<84) (s(:,7)>5)|(s(:,8)<34) ];

nrcol=size(reformat4,2);   %nr of variables in data
reformat4(1,nrcol+1:nrcol+7)=0;
for i=1:size(reformat4,1)
    t=max(p(s1(i,:)))+max(p(s2(i,:)))+max(p(s3(i,:)))+max(p(s4(i,:)))+max(p(s5(i,:)))+max(p(s6(i,:)));  %SUM OF ALL 6 CRITERIA

    if t
    reformat4(i,nrcol+1:nrcol+7)=    [max(p(s1(i,:))) max(p(s2(i,:))) max(p(s3(i,:))) max(p(s4(i,:))) max(p(s5(i,:))) max(p(s6(i,:))) t];
    end
end

% SIRS - at each timepoint |  need: temp HR RR PaCO2 WBC
a=zeros(5,1); % indices of vars used in SOFA
a(1)=find(ismember(reformat4t.Properties.VariableNames,{'Temp_C'}));
a(2)=find(ismember(reformat4t.Properties.VariableNames,{'HR'}));
a(3)=find(ismember(reformat4t.Properties.VariableNames,{'RR'}));
a(4)=find(ismember(reformat4t.Properties.VariableNames,{'paCO2'}));
a(5)=find(ismember(reformat4t.Properties.VariableNames,{'WBC_count'}));
s=table2array(reformat4t(:,a));

s1=[s(:,1)>=38| s(:,1)<=36];   %count of points for all criteria of SIRS
s2=[s(:,2)>90 ];
s3=[s(:,3)>=20|s(:,4)<=32];
s4=[s(:,5)>=12| s(:,5)<4];
reformat4(:,nrcol+8)=s1+s2+s3+s4;

% adds 2 cols for SOFA and SIRS, if necessary
if sum(ismember(reformat4t.Properties.VariableNames,{'SIRS'}))== 0
reformat4t(:,end+1:end+2)=array2table(0);
reformat4t.Properties.VariableNames(end-1:end)= {'SOFA','SIRS'};
end

% more IO corrections
 ii=reformat4t.input_total<0;
 reformat4t.input_total(ii)=0;
 ii=reformat4t.input_4hourly<0;
 reformat4t.input_4hourly(ii)=0;

% records values
reformat4t(:,end-1)=array2table(reformat4(:,end-1));
reformat4t(:,end)=array2table(reformat4(:,end));

%% ########################################################################
%                     CREATE FINAL MIMIC_TABLE
% ########################################################################

MIMICtable = reformat4t;

author: matthieukomorowski


text:
function crops_data = prepare_image(im)
% ------------------------------------------------------------------------
% caffe/matlab/+caffe/imagenet/ilsvrc_2012_mean.mat contains mean_data that
% is already in W x H x C with BGR channels
d = load('ilsvrc_2012_mean.mat');
mean_data = d.mean_data;
IMAGE_DIM = 256;
CROPPED_DIM = 224; % 224 for googLeNet , 227 for VGG and AlexNet

% Convert an image returned by Matlab's imread to im_data in caffe's data
% format: W x H x C with BGR channels
im_data = im(:, :, [3, 2, 1]);  % permute channels from RGB to BGR
im_data = permute(im_data, [2, 1, 3]);  % flip width and height
im_data = single(im_data);  % convert from uint8 to single
im_data = imresize(im_data, [IMAGE_DIM IMAGE_DIM], 'bilinear');  % resize im_data
im_data = im_data - mean_data;  % subtract mean_data (already in W x H x C, BGR)

% oversample (4 corners, center, and their x-axis flips)
crops_data = zeros(CROPPED_DIM, CROPPED_DIM, 3, 10, 'single');
indices = [0 IMAGE_DIM-CROPPED_DIM] + 1;
n = 1;
for i = indices
  for j = indices
    crops_data(:, :, :, n) = im_data(i:i+CROPPED_DIM-1, j:j+CROPPED_DIM-1, :);
    crops_data(:, :, :, n+5) = crops_data(end:-1:1, :, :, n);
    n = n + 1;
  end
end
center = floor(indices(2) / 2) + 1;
crops_data(:,:,:,5) = ...
  im_data(center:center+CROPPED_DIM-1,center:center+CROPPED_DIM-1,:);
crops_data(:,:,:,10) = crops_data(end:-1:1, :, :, 5);

author: zhoubolei

repo: CAM

save_path: github-repos/MATLAB/zhoubolei-CAM