text:
function S = prolateSurfaceArea(elli, varargin)
%PROLATESURFACEAREA  Approximated surface area of a prolate ellipsoid.
%
%   S = prolateSurfaceArea(R1,R2)
%
%   Example
%   prolateSurfaceArea
%
%   See also
%   geom3d, ellipsoidSurfaceArea, oblateSurfaceArea
%

% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2015-07-03,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2015 INRA - Cepia Software Platform.

%% Parse input argument

if size(elli, 2) == 7
    R1 = elli(:, 4);
    R2 = elli(:, 5);

elseif size(elli, 2) == 1 && ~isempty(varargin)
    R1 = elli(:, 1);
    R2 = varargin{1};
end

assert(R1 > R2, 'first radius must be larger than second radius');

% surface theorique d'un ellipsoide prolate
% cf http://fr.wikipedia.org/wiki/Ellipso%C3%AFde_de_r%C3%A9volution
e = sqrt(R1.^2 - R2.^2) ./ R1;
S = 2 * pi * R2.^2 + 2 * pi * R1 .* R2 .* asin(e) ./ e;

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%% DEMO 18: Arbitrary axis of rotation
%
%
%
% Some modenr CT geometires are starting to be a bit more complex, one of
% the common things being arbitrary axis of rotation i.e. the detector and the
% source can move not in a circular path, but in a "spherical" path.
%
% In TIGRE this has been implemented by defining the rotation with 3
% angles, specifically the ZYZ configuration of Euler angles.
%
%  This demo shows how to use it.
%
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% This file is part of the TIGRE Toolbox
% % Copyright (c) 2015, University of Bath and
%                     CERN-European Organization for Nuclear Research
%                     All rights reserved.
%
% License:            Open Source under BSD.
%                     See the full license at
%                     https://github.com/CERN/TIGRE/blob/master/LICENSE
%
% Contact:            tigre.toolbox@gmail.com
% Codes:              https://github.com/CERN/TIGRE/
% Coded by:           Ander Biguri
%--------------------------------------------------------------------------
%% Initialize

clear;
close all;
%% Define Geometry
%
% VARIABLE                                   DESCRIPTION                    UNITS
%-------------------------------------------------------------------------------------
geo.DSD = 1536;                             % Distance Source Detector      (mm)
geo.DSO = 1000;                             % Distance Source Origin        (mm)
% Detector parameters
geo.nDetector=[512; 512];					% number of pixels              (px)
geo.dDetector=[0.8; 0.8]; 					% size of each pixel            (mm)
geo.sDetector=geo.nDetector.*geo.dDetector; % total size of the detector    (mm)
% Image parameters
geo.nVoxel=[128;128;128];                   % number of voxels              (vx)

% a bit smaller than usual because the demo includes a very big detector
% angle for showcase
geo.sVoxel=[256;256;256]/1.5;               % total size of the image       (mm)

geo.dVoxel=geo.sVoxel./geo.nVoxel;          % size of each voxel            (mm)
% Offsets
geo.offOrigin =[0;0;0];                     % Offset of image from origin   (mm)
geo.offDetector=[0; 0];                     % Offset of Detector            (mm)

% Auxiliary
geo.accuracy=0.5;                           % Accuracy of FWD proj          (vx/sample)

geo.mode='cone';

%% Define angles
numProjs = 100;

anglesY=linspace(0,2*pi,numProjs);
anglesZ2=anglesY;
anglesZ1=pi*sin(linspace(0,2*pi,numProjs));
angles=[anglesZ1;anglesY;anglesZ2];
%% Get Image

head=headPhantom(geo.nVoxel);

%% Project

projections=Ax(head,geo,angles);

plotProj(projections,(1:100)*pi/180); % angle information not right in the title
%% Reconstruct:

% Note, FDK will not work.

imgSIRT = SIRT(projections,geo, angles,50);
imgCGLS = CGLS(projections,geo, angles,10);

plotImg([head imgCGLS imgSIRT] ,'dim',3)

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Max manipulability index ALONG A LINE.
% Use stomp like to optimize along a surface/line
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function experiment2B_K50_N30
close all;
global parameters

%STOMP PARAMETERS
%number of particles
K = 50;

%repeat experiment number of times
parameters.n_repeat = 3;
parameters.experiment_name = 'experiment2B_K50_N30_time.mat';
parameters.experiment_number = 2;
parameters.animate = 0;

%repeat the experiment E times
random_manips=[];
Gout = [];
ela = [];
for i=1:parameters.n_repeat
    close all
    tic;
    [pk, final_manip] = path_planning_SCO_4DOF(K);
    Gout{i}=pk;
    random_manips = [random_manips; final_manip];
    elapsed = toc;
    ela = [ela elapsed];
    save(parameters.experiment_name)
end

fprintf('TIC TOC mean: %g\n', mean(ela));
fprintf('TIC TOC cov: %g\n', cov(ela));
fprintf('TIC TOC 2sigma: %g\n', 2*sqrt(cov(ela)));

'ended'
parameters.experiment_name

repo: ARTE


text:
function rtk=udstate_rtkins(rtk,obsr,obsb,nav,ind)

global glc
tt=rtk.tt;

% update the position/velocity/attitude
rtk=udpos_rtkins(rtk,tt);

% update ionosphereic parameter
if rtk.opt.ionoopt==glc.IONOOPT_EST
    roverpos=blh2xyz(rtk.x(7:9))';
    bl=norm(roverpos-rtk.basepos);
    rtk=udion_rtkins(rtk,tt,bl,ind);
end

% update tropspheric parameter
if rtk.opt.tropopt>=glc.TROPOPT_EST
    rtk=udtrop_rtkins(rtk,tt,bl);
end

% update reciever inter-frequency bias for glonass
if rtk.opt.glomodear==2&&rtk.mask(2)==1
    rtk=udrcvbias_rtkins(rtk,tt);
end

% update the ambiguity
if rtk.opt.mode>glc.PMODE_DGNSS
    rtk=udbias_rtkins(rtk,obsr,obsb,nav,ind);
end

return

repo: GINav


text:
function vw = computeMeanMap(vw,scanList,forceSave)
% Computing the mean functional image for each tSeries
%
%  vw = computeMeanMap(vw,[scanList],[forceSave])
%
% The mean functional images are combined into a parameter map and
% calls setParameterMap to set vw.map = meanMap.
%
% scanList:
%   0 - do all scans
%   number or list of numbers - do only those scans
%   default - prompt user via selectScans dialog
%
% forceSave: 1 = true (overwrite without dialog)
%            0 = false (query before overwriting)
%           -1 = do not save
%
% If you change this function make parallel changes in:
%    computeCorAnal, computeResStdMap, computeStdMap
%
% djh, 12/30/98
% djh, 2/22/2001 updated to version 3
% ras, 01/05, added forceSave flag
% ras 10/05, checks if the meanMap file is saved already

if notDefined('forceSave'), forceSave = 0; end

% nScans = numScans(vw);
nScans = viewGet(vw,'numscans');

if strcmp(vw.mapName,'meanMap')
    % If exists, initialize to existing map
    map=vw.map;
elseif exist(fullfile(dataDir(vw),'meanMap.mat'),'file')
    % load from the mean map file
    load(fullfile(dataDir(vw),'meanMap.mat'),'map')
else
    % Otherwise, initialize to empty cell array
    map = cell(1,nScans);
end

% (Re-)set scanList
if ~exist('scanList','var')
    scanList = er_selectScans(vw);
elseif scanList == 0
    scanList = 1:nScans;
end
if isempty(scanList),  error('Analysis aborted: scan list is empty.'); end

% Compute it
waitHandle = mrvWaitbar(0,'Computing mean images from the tSeries.  Please wait...');
ncScans = length(scanList);
for iScan = 1:ncScans
    scan = scanList(iScan);
    dims = viewGet(vw, 'sliceDims', scan);
    map{scan} = NaN*ones(dataSize(vw,scan));
    switch lower(viewGet(vw, 'viewType'))
        case 'inplane'
            % since data in the inplane view is stored as a nifti, we can
            % get the whole data slab in one call to load tseries rather
            % than looping over slices.
            [~, nii] = loadtSeries(vw,scan);
            tSeries = niftiGet(nii, 'data');
            timeDim = length(size(tSeries));
            map{scan} = mean(tSeries, timeDim);

        otherwise
            for slice = sliceList(vw,scan)
                tSeries = loadtSeries(vw,scan,slice);
                nValid = sum(isfinite(tSeries));
                tSeries(isnan(tSeries(:))) = 0;

                % if there is one time point, we will have problems, so duplicate
                % the data
                if size(tSeries,1) == 1,
                    tSeries = repmat(tSeries, 3, 1);
                    nValid = sum(isfinite(tSeries));
                end

                tmp = sum(tSeries) ./ nValid;
                if strcmp(vw.viewType,'Inplane')
                    map{scan}(:,:,slice) = reshape(tmp,dims);
                else
                    map{scan}(:,:,slice) = tmp;
                end
            end
    end
    mrvWaitbar(scan/ncScans)
end
close(waitHandle);

% Set Parameter MAp
vw = setParameterMap(vw, map, 'meanMap');

% Save file
if forceSave >= 0, saveParameterMap(vw, [], forceSave); end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% function [x1,x2,P1,P2] = make_synthetic_data(N,range_model,range_image,verbose)
% makes a synthetic 3D model and makes two projective images of the model
% using arbitrary projective cameras Pi=[M|t]. The second cameras are normalized
% to Pi = Pi / Pi_{34}.
% inputs:
%               N               1x1     (optional) the maximum number of correspondences to generate.
%               range_model     2x1     (optional) defines a bounding box (starting from origin) containing the model
%               range_image     2x1     (optional) the size of the output image(scaling on correspondences)
%               verbose         1x1     (optional) plots the generated model and images
% outputs:
%               x1              3x1     homogeneous coordinates of the correspondences in image 1
%               x2              3x1     homogeneous coordinates of the correspondences in image 2
%               P1              3x4     the camera that generates x1(up to an affine transformation)
%               P2              3x4     the camera that generates x2(up to an affine transformation)
%
% Author: Omid Aghazadeh, KTH(Royal Institute of Technology), 2010/05/09
function [x1,x2,P1,P2] = make_synthetic_data(N,range_model,range_image,verbose)
if nargin < 1, N = 100; end
if nargin < 2, range_model = [10;10;10]; end
if nargin < 3, range_image = [640;480]; end
if nargin < 4, verbose = 1; end
X = rand(3,N); X = (X - repmat(min(X,[],2),1,N)) .* repmat(range_model ./ (max(X,[],2) - min(X,[],2)),1,N);
X = [X; ones(1,N)];
P1 = rand(3,4); P1(end) = 1;
P2 = rand(3,4); P2(end) = 1; %avoid cameras at infinity
x1 = P1*X; x1 = x1./repmat(x1(3,:),3,1); x2 = P2*X; x2 = x2./repmat(x2(3,:),3,1);
ixinv = sum(isinf(x1) | isnan(x1) | isinf(x2) | isnan(x2) | abs(x1) > 1e3 | abs(x2) > 1e3,1 )>0 ;
x1 = x1(:,~ixinv); x2 = x2(:,~ixinv); N = size(x1,2);
x1 = (x1-repmat([min(x1(1:2,:),[],2);0],1,N)).* repmat([range_image./(max(x1(1:2,:),[],2) - min(x1(1:2,:),[],2));1],1,N);

x2 = (x2-repmat([min(x2(1:2,:),[],2);0],1,N)).* repmat([range_image./(max(x2(1:2,:),[],2) - min(x2(1:2,:),[],2));1],1,N);
if verbose
    figure(1); subplot(1,3,1); plot3(X(1,:),X(2,:),X(3,:),'b.'); title('generated 3D model');
    subplot(1,3,2); plot(x1(1,:),x1(2,:),'rx'); title('generated first image');
    subplot(1,3,3); plot(x2(1,:),x2(2,:),'gx'); title('generated second image');
end
end

author: Sable

repo: mcbench-benchmarks


text:
%% AI Clinician Identifiying MIMIC-III sepsis cohort

% (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
% IDENTIFIES THE COHORT OF PATIENTS WITH SEPSIS in MIMIC-III

% PURPOSE:
% ------------------------------
% This creates a list of icustayIDs of patients who develop sepsis at some point
% in the ICU. records charttime for onset of sepsis. Uses sepsis3 criteria

% STEPS:
% -------------------------------
% IMPORT DATA FROM CSV FILES
% FLAG PRESUMED INFECTION
% PREPROCESSING
% REFORMAT in 4h time slots
% COMPUTE SOFA at each time step
% FLAG SEPSIS

% note: the process generates the same features as the final MDP dataset, most of which are not used to compute SOFA
% External files required: Reflabs, Refvitals, sample_and_hold (all saved in reference_matrices.mat file)

% 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

%% ########################################################################
% IMPORT ALL DATA

tic
abx=table2array(readtable('D:/exportdir/abx.csv'));
culture=table2array(readtable('D:/exportdir/culture.csv'));
microbio=table2array(readtable('D:/exportdir/microbio.csv'));
demog=(readtable('D:/exportdir/demog.csv'));
ce010=table2array(readtable('D:/exportdir/ce010000.csv'));
ce1020=table2array(readtable('D:/exportdir/ce1000020000.csv'));
ce2030=table2array(readtable('D:/exportdir/ce2000030000.csv'));
ce3040=table2array(readtable('D:/exportdir/ce3000040000.csv'));
ce4050=table2array(readtable('D:/exportdir/ce4000050000.csv'));
ce5060=table2array(readtable('D:/exportdir/ce5000060000.csv'));
ce6070=table2array(readtable('D:/exportdir/ce6000070000.csv'));
ce7080=table2array(readtable('D:/exportdir/ce7000080000.csv'));
ce8090=table2array(readtable('D:/exportdir/ce8000090000.csv'));
ce90100=table2array(readtable('D:/exportdir/ce90000100000.csv'));
labU=[ table2array(readtable('D:/exportdir/labs_ce.csv')) ; table2array(readtable('D:/exportdir/labs_le.csv'))  ];
MV=table2array(readtable('D:/exportdir/mechvent.csv'));
inputpreadm=table2array(readtable('D:/exportdir/preadm_fluid.csv'));
inputMV=table2array(readtable('D:/exportdir/fluid_mv.csv'));
inputCV=table2array(readtable('D:/exportdir/fluid_cv.csv'));
vasoMV=table2array(readtable('D:/exportdir/vaso_mv.csv'));
vasoCV=table2array(readtable('D:/exportdir/vaso_cv.csv'));
UOpreadm=table2array(readtable('D:/exportdir/preadm_uo.csv'));
UO=table2array(readtable('D:/exportdir/uo.csv'));
toc

%% ########################################################################
%                       INITIAL DATA MANIPULATIONS
% #########################################################################

ii=isnan(microbio(:,3));  %if charttime is empty but chartdate isn't
microbio(ii,3)=microbio(ii,4);   %copy time
microbio( :,4)=[];    %delete chardate
% Add empty col in microbio (# 3 and #5)
microbio(:,4)=microbio(:,3);
microbio(:,[3 5])=0;
% Combine both tables for micro events
bacterio = [microbio ; culture];
% correct NaNs in DEMOG
demog.morta_90(isnan(demog.morta_90))=0;
demog.morta_hosp(isnan(demog.morta_hosp))=0;
demog.elixhauser(isnan(demog.elixhauser))=0;

% compute normalized rate of infusion
% if we give 100 ml of hypertonic fluid (600 mosm/l) at 100 ml/h (given in 1h) it is 200 ml of NS equivalent
% so the normalized rate of infusion is 200 ml/h (different volume in same duration)
inputMV(:,8)=inputMV(:,7).*inputMV(:,6)./inputMV(:,5);

% fill-in missing ICUSTAY IDs in bacterio
for i=1:size(bacterio,1)
if bacterio(i,3)==0   %if missing icustayid
    o=bacterio(i,4);  %charttime
    subjectid=bacterio(i,1);
    hadmid=bacterio(i,2);
   ii=find(demog.subject_id==subjectid);
   jj=find(demog.subject_id==subjectid & demog.hadm_id==hadmid);
    for j=1:numel(ii)
        if o>=demog.intime(ii(j))-48*3600 & o<=demog.outtime(ii(j))+48*3600
            bacterio(i,3)=demog.icustay_id(ii(j));
        elseif numel(ii)==1   %if we cant confirm from admission and discharge time but there is only 1 admission: it's the one!!
            bacterio(i,3)=demog.icustay_id(ii(j));
        end
    end
end
end
toc

for i=1:size(bacterio,1)
if bacterio(i,3)==0   %if missing icustayid
    subjectid=bacterio(i,1);
    hadmid=bacterio(i,2);
    jj=find(demog.subject_id==subjectid & demog.hadm_id==hadmid);
    if numel(jj)==1
        bacterio(i,3)=demog.icustay_id(jj);
    end
end
end

% fill-in missing ICUSTAY IDs in ABx

for i=1:size(abx,1)
if isnan(abx(i,2))
    o=abx(i,3);  %time of event
    hadmid=abx(i,1);
    ii=find(demog.hadm_id==hadmid);   %row in table demographics
    for j=1:numel(ii)
        if o>=demog.intime(ii(j))-48*3600 & o<=demog.outtime(ii(j))+48*3600
            abx(i,2)=demog.icustay_id(ii(j));
        elseif numel(ii)==1   %if we cant confirm from admission and discharge time but there is only 1 admission: it's the one!!
            abx(i,2)=demog.icustay_id(ii(j));
        end
    end
end
end

%% ########################################################################
%    find presumed onset of infection according to sepsis3 guidelines
% ########################################################################

% METHOD:
% I loop through all the ABx given, and as soon as there is a sample present
% within the required time criteria I pick this flag and break the loop.

onset=zeros(100000,3);

for icustayid=1:100000

    ab=abx(abx(:,2)==icustayid+200000,3);   %start time of abx for this icustayid
    bact=bacterio(bacterio(:,3)==icustayid+200000,4);  %time of sample
    subj_bact=bacterio(bacterio(:,3)==icustayid+200000,1);  %subjectid

    if ~isempty(ab) & ~isempty(bact)   %if we have data for both: proceed

      D = pdist2(ab, bact)/3600;  %pairwise distances btw antibio and cultures, in hours

      for i=1:size(D,1)  % looping through all rows of AB given, from early to late
        [M,I] = min(D(i,:));   %minimum distance in this row
        ab1=ab(i);       %timestamp of this value in list of antibio
        bact1=bact(I);      %timestamp in list of cultures

        if M<=24 & ab1<=bact1      %if ab was first and delay < 24h
            onset(icustayid,1)=subj_bact(1);   %subject_id
            onset(icustayid,2)=icustayid;       % icustay_id
            onset(icustayid,3)=ab1;     %onset of infection = abx time
              icustayid
            break
        elseif M<=72 & ab1>=bact1    %elseif sample was first and delay < 72h
            onset(icustayid,1)=subj_bact(1);
            onset(icustayid,2)=icustayid;
            onset(icustayid,3)=bact1;       %onset of infection = sample time
            break
        end
      end
    end
end
toc

%sum of records found
sum(onset(:,3)>0)

%% Replacing item_ids with column numbers from reference tables

% replace itemid in labs with column number
% this will accelerate process later

tic
for i=1:size(labU,1)
[~,locb]=ismember(Reflabs,labU(i,3));
labU(i,3)=find(max(locb')');
end
toc

% replace itemid in vitals with col number

for i=1:size(ce010,1)
[~,locb]=ismember(Refvitals,ce010(i,3));ce010(i,3)=find(max(locb')');
end
for i=1:size(ce1020,1)
[~,locb]=ismember(Refvitals,ce1020(i,3));ce1020(i,3)=find(max(locb')');
end
for i=1:size(ce2030,1)
[~,locb]=ismember(Refvitals,ce2030(i,3));ce2030(i,3)=find(max(locb')');
end
for i=1:size(ce3040,1)
[~,locb]=ismember(Refvitals,ce3040(i,3));ce3040(i,3)=find(max(locb')');
end
for i=1:size(ce4050,1)
[~,locb]=ismember(Refvitals,ce4050(i,3));ce4050(i,3)=find(max(locb')');
end
for i=1:size(ce5060,1)
[~,locb]=ismember(Refvitals,ce5060(i,3));ce5060(i,3)=find(max(locb')');
end
for i=1:size(ce6070,1)
[~,locb]=ismember(Refvitals,ce6070(i,3));ce6070(i,3)=find(max(locb')');
end
for i=1:size(ce7080,1)
[~,locb]=ismember(Refvitals,ce7080(i,3));ce7080(i,3)=find(max(locb')');
end
for i=1:size(ce8090,1)
[~,locb]=ismember(Refvitals,ce8090(i,3));ce8090(i,3)=find(max(locb')');
end
for i=1:size(ce90100,1)
[~,locb]=ismember(Refvitals,ce90100(i,3));ce90100(i,3)=find(max(locb')');
end

%% ########################################################################
%           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 -48 -> +24 because that's for sepsis3 cohort defintion!!
% I need different time period for the MDP (-24 -> +48)

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

tic
for icustayid=1:100000
qst=onset(icustayid,3); %flag for presumed infection
if qst>0  % if we have a flag
d1=table2array(demog(demog.icustay_id==icustayid+200000,[11 5])); %age of patient + discharge time

if d1(1)>6574  % if older than 18 years old

    waitbar(icustayid/100000,h,icustayid/1000) %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;%(locb(:,1)); %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; %flag for presumed infection / this is time of sepsis if SOFA >=2 for this patient
%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)=d1(2); %dischargetime

end
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;

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

%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
% ########################################################################

% 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) || corrected for tonicity

        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 % -52 until +28 = 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

          % #####################   DISCUSS ADDING STUFF HERE / RANGE, MIN, MAX ETC   ################

        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);

%% ########################################################################
%    CONVERT TO TABLE AND DELETE VARIABLES WITH EXCESSIVE MISSINGNESS
% ########################################################################

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;
miss=sum(isnan(reformat2))./size(reformat2,1);

% if values have less than 70% missing values (over 30% of values present): I keep them
reformat3t=reformat2t(:,[true(1,11) miss(12:74)<0.70 true(1,11)]) ;

%% ########################################################################
%             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...
% ########################################################################

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

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

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

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

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

% 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));

% 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

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

%% ########################################################################
%                            EXCLUSION OF SOME PATIENTS
% ########################################################################

numel(unique(reformat4t.icustayid))  %count before

% check for patients with extreme UO = outliers = to be deleted (>40 litres of UO per 4h!!)
a=find(reformat4t.output_4hourly>12000);
i=unique(reformat4t.icustayid(a));
i=find(ismember(reformat4t.icustayid,i));
reformat4t(i,:)=[];

% some have bili = 999999
a=find(reformat4t.Total_bili>10000);
i=unique(reformat4t.icustayid(a));
i=find(ismember(reformat4t.icustayid,i));
reformat4t(i,:)=[];

% check for patients with extreme INTAKE = outliers = to be deleted (>10 litres of intake per 4h!!)
a=find(reformat4t.input_4hourly>10000);
i=unique(reformat4t.icustayid(a));  % 28 ids
i=find(ismember(reformat4t.icustayid,i));
reformat4t(i,:)=[];

% #### exclude early deaths from possible withdrawals ####
% stats per patient
q=reformat4t.bloc==1;
% fence_posts=find(q(:,1)==1);
num_of_trials=numel(unique(reformat4t.icustayid));%size(fence_posts,1);
a=array2table([reformat4t.icustayid reformat4t.mortality_90d reformat4t.max_dose_vaso reformat4t.SOFA]);
a.Properties.VariableNames={'id','mortality_90d','vaso','sofa'};
d=grpstats(a,'id','max');

%finds patients who match our criteria
e=zeros(num_of_trials,1);
for i=1:num_of_trials
    if d.max_mortality_90d(i) ==1
    ii=reformat4t.icustayid==d.id(i) & reformat4t.bloc==d.GroupCount(i);  %last row for this patient
    e(i)=sum((reformat4t.max_dose_vaso(ii)==0 & d.max_vaso(i)>0.3 & reformat4t.SOFA(ii)>=d.max_sofa(i)/2))>0;
    end
end
r=d.id(e==1 & d.GroupCount<20); % ids to be removed
ii=ismember(reformat4t.icustayid,r);
reformat4t(ii,:)=[];

% exclude patients who died in ICU during data collection period
ii=reformat4t.bloc==1&reformat4t.died_within_48h_of_out_time==1& reformat4t.delay_end_of_record_and_discharge_or_death<24;
ii=ismember(icustayidlist,reformat4t.icustayid(ii));
reformat4t(ii,:)=[];

numel(unique(reformat4t.icustayid))   %count after

%% #######################################################################
%                       CREATE SEPSIS COHORT
% ########################################################################

% create array with 1 row per icu admission
% keep only patients with flagged sepsis (max sofa during time period of interest >= 2)
% we assume baseline SOFA of zero (like other publications)

sepsis=zeros(30000,5);
irow=1;

tic
for icustayid=1:100000
    ii=find(ismember(reformat4t.icustayid,icustayid));
    if mod(icustayid,10000)==0;disp([num2str(icustayid/1000), ' %']);end
    if ii

         sofa=reformat4t.SOFA(ii);
         sirs=reformat4t.SIRS(ii);
         sepsis(irow,1)=icustayid+200000;
         sepsis(irow,2)=reformat4t.mortality_90d(ii(1)); % 90-day mortality
         sepsis(irow,3)=max(sofa);
         sepsis(irow,4)=max(sirs);
         sepsis(irow,5)=qstime(icustayid);   %time of onset of sepsis
         irow=irow+1;
    end
end
toc
sepsis(irow:end,:)=[];

sepsis=array2table(sepsis);
sepsis.Properties.VariableNames={'icustayid','morta_90d','max_sofa','max_sirs','sepsis_time'};

% delete all non-sepsis
sepsis(sepsis.max_sofa<2,:)=[];

% final count
size(sepsis,1)

%save cohort
writetable(sepsis,'sepsis_mimiciii.csv','Delimiter',',');

author: matthieukomorowski


text:
% Semi-NMF: Semi Non-negative Matrix Factorization (Trigeorgis et al. 2014)
% process_video('NMF', 'Semi-NMF', 'dataset/demo.avi', 'output/demo_Semi-NMF.avi');
rank = 10;
[W,H] = seminmf(M, rank);
L = W * H;
S = M - L;

author: andrewssobral

repo: lrslibrary


text:
function classifier = laprlsc(options,data)
% {laprlsc} trains a Laplacian RLSC.
%
%      classifier = laprlsc(options,data)
%
%      Note: This function is only a wrapper to 'lapsvmp', where the use of
%            a Hinge function for labeled data is disabled. See 'lapsvmp'
%            for the description of all required parameters.
%
% Author: Stefano Melacci (2009)
%         mela@dii.unisi.it

options.Hinge=0;
classifier=lapsvmp(options,data);
classifier.name='laprlsc';

author: viggin

repo: domain-adaptation-toolbox


text:
classdef QuadrilaterSubTriangulator < handle

    properties (Access = public)

    end

    properties (Access = private)
        subTriangulation
    end

    properties (Access = private)
        localIncidenceMatrix
        coord
        connec
        xLocalVertex
    end

    methods (Access = public)

        function obj = QuadrilaterSubTriangulator(cParams)
            obj.init(cParams)
        end

        function [connR,xCoordsIsoBoundaryR] = compute(obj)

            for is = 1:2
                subTriang = obj.subTriangulation(:,:,is);
                connecL             = obj.computeSubtriangulation(subTriang);
                conne               = obj.divideInSubTriangle(obj.connec,connecL);

                xB = permute(obj.xLocalVertex,[3 2 1]);
                xCoordsIsoBoundaryS = obj.divideInSubTriangle(xB,connecL);

                xCoordsIsoBoundaryS = permute(xCoordsIsoBoundaryS,[4 1 2 3]);

                conn1(:,:,is) = conne(:,:,1);
                conn2(:,:,is) = conne(:,:,2);
                %connRi(:,:,is) = cat(1,conn1,conn2);

                nodesSubCases(:,:,is,:) = permute(conne,[3 2 1]);

                xCoordsIsoBoundary1(:,:,:,is) = xCoordsIsoBoundaryS(:,:,:,1);
                xCoordsIsoBoundary2(:,:,:,is) = xCoordsIsoBoundaryS(:,:,:,2);

                %  nodesSubCases(:,:,1,:) =
            end

            s.nodesSubCases = nodesSubCases;
            s.coord = obj.coord;
            b = BestSubCellCaseSelector(s);
            imax = b.compute();

            conR1 = zeros(size(conn1(:,:,1)));
            conR2 = zeros(size(conn1(:,:,1)));
            xCoordsIsoBoundaryR1 = zeros(size(xCoordsIsoBoundary1(:,:,:,1)));
            xCoordsIsoBoundaryR2 = zeros(size(xCoordsIsoBoundary1(:,:,:,1)));
            for is = 1:2
                isIndex = imax == is;
                con1 = conn1(:,:,is);
                con2 = conn2(:,:,is);

                conR1(isIndex,:) = con1(isIndex,:);
                conR2(isIndex,:) = con2(isIndex,:);

                xR1 = xCoordsIsoBoundary1(:,:,:,is);
                xR2 = xCoordsIsoBoundary2(:,:,:,is);
                xCoordsIsoBoundaryR1(:,isIndex,:) = xR1(:,isIndex,:);
                xCoordsIsoBoundaryR2(:,isIndex,:) = xR2(:,isIndex,:);

            end
            %   connR = cat(1,conR1,conR2);
            connR = cat(1,conne(:,:,1),conne(:,:,2));

            xCoordsIsoBoundaryRN = cat(2,xCoordsIsoBoundaryR1,xCoordsIsoBoundaryR2);
            xCoordsIsoBoundaryR = permute(xCoordsIsoBoundaryRN,[1 3 2]);
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.localIncidenceMatrix = cParams.localIncidenceMatrix;
            obj.coord                = cParams.coord;
            obj.connec               = cParams.connec;
            obj.xLocalVertex         = cParams.xLocalVertex;
            obj.subTriangulation(:,:,1) = [1 2 3; 3 4 1];
            obj.subTriangulation(:,:,2) = [2 3 4; 4 1 2];
        end

        function [xCoordsIsoBoundaryR] = divideInSubTriangle(obj,xB,connecL)

     %       connB = obj.localToGlobal(connecL,connec);

            nElem = size(connecL,1);
            xCoordsIsoBoundaryR = zeros(nElem,3,2,size(xB,3));
            for idim = 1:size(xB,3)
                xBi = (xB(:,:,idim));
                xC = obj.localToGlobal(connecL,xBi);
                xCoordsIsoBoundaryR(:,:,:,idim) = xC;
            end
        end

        function [connecG] = localToGlobal(obj,connecL,dataToTriang)

            for itriang = 1:2
                for i = 1:3
                    nodei1 = connecL(:,i,itriang);
                    index1 = sub2ind(size(dataToTriang),[1:size(connecL,1)]',nodei1);
                    connecG(:,i,itriang) = dataToTriang(index1);
                end
            end

        end

        function nodes = obtainOrderedNodes(obj)
            incMatrix = obj.localIncidenceMatrix;
            nElem = size(incMatrix,1);
            incMatrix = permute(incMatrix,[1 3 2]);
            edge1 = incMatrix(:,:,1);
            edge2 = incMatrix(:,:,2);
            edge3 = incMatrix(:,:,3);
            edge4 = incMatrix(:,:,4);
            areIntersected(:,1) = obj.areEdgesIntersected(edge1,edge2);
            areIntersected(:,2) = obj.areEdgesIntersected(edge1,edge3);
            areIntersected(:,3) = obj.areEdgesIntersected(edge1,edge4);

            integerCases = obj.computeIntegerCases(areIntersected);

            nodes = zeros(nElem,4);
            cases = [3 5 6];
            nodesCases = [1 2 4 3; 1 2 3 4; 1 3 2 4];
            for icase = 1:3
                isCase = integerCases == cases(icase);
                for inode = 1:4
                    nodes(isCase,inode) = nodesCases(icase,inode);
                end
            end

        end

        function connec = computeSubtriangulation(obj,subTriang)
            nodes = obj.obtainOrderedNodes();
            connec = zeros(size(nodes,1),3,2);
            for inode = 1:3
                for itriangle = 1:2
                    connec(:,inode,itriangle) = nodes(:,subTriang(itriangle,inode));
                end
            end
        end

    end

    methods (Access = private, Static)

        function [areInt] = areEdgesIntersected(edge1,edge2)
            isNode11 = edge1(:,1) == edge2(:,1);
            isNode12 = edge1(:,1) == edge2(:,2);
            isNode21 = edge1(:,2) == edge2(:,1);
            isNode22 = edge1(:,2) == edge2(:,2);
            areInt = any([isNode11 isNode12 isNode21 isNode22],2);
        end

        function integerCases = computeIntegerCases(isNeigbour)
            nodes = isNeigbour;
            nnode = size(nodes,2);
            nodePos = (1:nnode) - 1;
            pow2vector = 2.^(nodePos);
            d = pow2vector*nodes';
            integerCases = d;
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
classdef SettingsFilterP0 < AbstractSettings
    properties
        shape
        conec
        quadr
    end
    methods
    end
end


text:
function AddSpikeTimes(spikes)

%AddSpikeTimes - Add spike timestamps.
%
%  USAGE
%
%    AddSpikeTimes(spikes)
%
%  SEE
%
%    See also ThresholdSpikes.

% Copyright (C) 2012 by Michaël Zugaro
%
% 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.

global DATA;
if isempty(DATA),
	error('No session defined (did you forget to call SetCurrentSession? Type ''help Data'' for details).');
end

% Parse parameters
if ~issamples(spikes,'#2'),
	error('Incorrect spikes (type ''help AddSpikeTimes'' for details).');
end

DATA.spikes = [DATA.spikes;spikes];

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [] = mci_compare_sensitivities (model,pars)
% Compare methods for sensitivity computation
% FORMAT [] = mci_compare_sensitivities (model,pars)
%
% model     'phase', 'nmm-r2p2'
% pars      vector indicating which sensitivities to plot
%           eg. [1,2,..,Np] (default) for all parameters
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: mci_compare_sensitivities.m 6548 2015-09-11 12:39:47Z will $

[P,M,U,Y] = mci_compare_setup (model);

if nargin < 2 | isempty(pars)
    pars=[1:M.Np];
end
Np=length(pars);

disp(' ');
disp('Plot sensitivity to parameter changes as computed by:');
disp('(1) Forward equations based on Matlab''s ODE suite (red)');
disp('(2) Forward equations based on Sundials (blue)');

[G,sy_m,st] = spm_mci_sens (P,M,U);

[G,sy,st] = spm_mci_sens_sun (P,M,U);

if st==-1
    disp('Problem with integration');
    return
end

mci_plot_outputs(M,G);

hs=figure;
set(hs,'Name','Sensitivities');
k=1; lw=2;
for i=1:M.l,
    for p=1:Np,
        j=pars(p);
        subplot(M.l,Np,k);
        plot(M.t,squeeze(sy(:,i,j)),'LineWidth',lw);
        hold on
        plot(M.t,squeeze(sy_m(:,i,j)),'r','LineWidth',lw);
        grid on
        title (sprintf('dy(%d)/dp(%d)',i,j));
        xlabel('Time');
        k=k+1;
    end
end

author: spm


text:
function stk = bin_stack(w,bin,ovr)

   %BIN_STACK: Stack input waveforms with N waveforms per stack. The number N
   %    is specified by input 'bin'. The input 'ovr' specifies the amount of
   %    overlap between bins.
   %
   %USAGE: stk = bin_stack(w,bin,ovr)
   %    If input w contains 100 waveforms then:
   %    stk = bin_stack(w,20,0) % stk has 5 stacks of 20 waveforms
   %    stk = bin_stack(w,5,0)  % stk has 20 stacks of 5 waveforms
   %    stk = bin_stack(w,20,10)  % stk has 9 stacks of 20 waveforms
   %
   %INPUTS:  w   - Input waveform
   %         bin - Bin size of stacked waveforms
   %         ovr - Number of overlapping events between bins
   %
   %OUTPUTS: stk - Array of stacked waveforms

   % Author: Dane Ketner, Alaska Volcano Observatory
   % $Date$
   % $Revision$

   nw = numel(w);
   inc = bin-ovr;
   n=1;
   while 1
      if (n-1)*inc+bin > nw
         return
      else
         stk(n) = stack(w((n-1)*inc+1:(n-1)*inc+bin));
         n=n+1;
      end
   end
end

author: geoscience-community-codes

repo: GISMO


text:
function pi = predlssvm(model,Xt,alpha,conftype)

% Construction of bias corrected 100(1-\alpha)% pointwise or
% simultaneous prediction intervals
%
% >> pi = predlssvm({X,Y,type,gam,kernel_par,kernel,preprocess}, Xt, alpha,conftype)
% >> pi = predlssvm(model, Xt, alpha, conftype)
%
% This function calculates bias corrected 100(1-\alpha)% pointwise or
% simultaneous prediction intervals. The procedure support homoscedastic
% data sets as well heteroscedastic data sets. The construction of the
% prediction intervals are based on the central limit theorem for linear
% smoothers combined with bias correction and variance estimation.
%
% 1. Using the functional interface:
%
%
% >> pi = predlssvm({X,Y,type,gam,kernel_par,kernel,preprocess})
% >> pi = predlssvm({X,Y,type,gam,kernel_par,kernel,preprocess}, alpha)
% >> pi = predlssvm({X,Y,type,gam,kernel_par,kernel,preprocess}, alpha, conftype)
%
%
%      Outputs
%        pi            : N x 2 matrix containing the lower and upper prediction intervals
%
%      Inputs
%        X             : N x d matrix with the inputs of the training data
%        Y             : N x 1 vector with the outputs of the training data
%        type          : 'function estimation' ('f') or 'classifier' ('c')
%        gam           : Regularization parameter
%        sig2          : Kernel parameter(s) (bandwidth in the case of the 'RBF_kernel')
%        kernel(*)     : Kernel type (by default 'RBF_kernel')
%        preprocess(*) : 'preprocess'(*) or 'original'
%        alpha(*)      : Significance level (by default 5%)
%        conftype(*)   : Type of prediction interval 'pointwise' or 'simultaneous' (by default 'simultaneous')
%
% 2. Using the object oriented interface:
%
%
% >> pi = predlssvm(model)
% >> pi = predlssvm(model, alpha)
% >> pi = predlssvm(model, alpha, conftype)
%
%
%      Outputs
%        pi          : N x 2 matrix containing the lower and upper confidence intervals
%
%      Inputs
%        model       : Object oriented representation of the LS-SVM model
%        alpha       : Significance level (by default 5%)
%        conftype    : Type of prediction interval 'pointwise' or 'simultaneous' (by default 'simultaneous')
%
%
%  See also:
%    trainlssvm, simlssvm, cilssvm

% Copyright (c) 2011,  KULeuven-ESAT-SCD, License & help @ http://www.esat.kuleuven.be/sista/lssvmlab

if iscell(model)
    model = initlssvm(model{:});
end

if nargin == 1
    error('Please specify a test set');
elseif nargin <= 2
    alpha = 0.05;
    conftype = 'simul';

elseif nargin <= 3
    conftype = 'simul';
end

if model.preprocess(1)=='p'
    [x,y] = postlssvm(model,model.xtrain,model.ytrain);
else
    x = model.xtrain; y = model.ytrain;
end

% train model
model = trainlssvm(model);

% smoother on test
S = smootherlssvm(model,Xt);

% bias: double smoothing with fourt order kernel RBF4
modelb = initlssvm(x,y,'f',[],[],'RBF4_kernel','o');
modelb = tunelssvm(modelb,'simplex','crossvalidatelssvm',{10,'mse'});
modelb = trainlssvm(modelb);

% output on test Yphat
Yphat = simlssvm(modelb,Xt);

% output of the smoother based on fourth order kernel on training data Yhat
Yhat = simlssvm(modelb,x);

biascorr = S*Yhat-Yphat;

%1) estimate variance nonparametrically
[sigma2,modele] = varest(model);
% smoother on training data
L = smootherlssvm(model);
sigma2t = varfun(modele,L,Xt);

%2) calculate var-cov matrix
S = S*diag(sigma2)*S';

%2b) find standardized absolute maxbias
delta = max(abs(biascorr./sqrt(diag(S))));

%3) pointwise or simultaneous
if conftype(1)=='s'
    if model.preprocess(1)=='p'
        model.xtrain = prelssvm(model,x);
    end
    z = tbform(model,alpha) + delta;
elseif conftype(1)=='p'
    z = norminv(alpha/2);
    Yphat = Yphat - biascorr;
else
    error('Wrong type of confidence interval. Please choose ''pointwise'' or ''simultaneous''');
end

V = sqrt(sigma2t+diag(S));
pi = [Yphat+z*V Yphat-z*V];

function V = varfun(modele,L,xt)
% estimate variance op new data, given the smooth of squared residuals
% "modele" and smoother matrix of the original smooth L
V = max(simlssvm(modele,xt),0);
% make estimate of V unbiased in homoscedastic case if regression
% estimate is unbiased
S = smootherlssvm(modele,xt);
V = V./(ones(size(xt,1),1)+S*diag(L*L'-2*L));

function [var,modele] = varest(model)

% if preprocessed data, construct original data
if model.preprocess(1)=='p'
    [x,y] = postlssvm(model,model.xtrain,model.ytrain);
else
    x = model.xtrain; y = model.ytrain;
end

model = trainlssvm(model);

Yh = simlssvm(model,x);

% Squared normalized residuals
e2 = (y-Yh).^2;

% Make variance model
if model.nb_data <= 200
    costfun = 'leaveoneoutlssvm'; costargs = {'mae'};
else
    costfun = 'crossvalidatelssvm'; costargs = {10,'mae'};
end
modele = initlssvm(x,e2,'f',[],[],'RBF_kernel');
modele = tunelssvm(modele,'simplex',costfun,costargs);
modele = trainlssvm(modele);

% variance model
var = max(simlssvm(modele,x),0);

% make estimate of Var unbiased in homoscedastic case if regression
% estimate is unbiased
L = smootherlssvm(model);
S = smootherlssvm(modele);

var = var./(ones(size(x,1),1)+S*diag(L*L'-L-L'));

author: peterhcharlton

repo: RRest


text:
function truncated_normal_ab_variance_test ( )

%*****************************************************************************80
%
%% TRUNCATED_NORMAL_AB_VARIANCE_TEST tests TRUNCATED_NORMAL_AB_VARIANCE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 March 2015
%
%  Author:
%
%    John Burkardt
%
  sample_num = 1000;
  seed = 123456789;
  a = 50.0;
  b = 150.0;
  mu = 100.0;
  sigma = 25.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TRUNCATED_NORMAL_AB_VARIANCE_TEST\n' );
  fprintf ( 1, '  TRUNCATED_NORMAL_AB_VARIANCE computes the variance.\n' );
  fprintf ( 1, '  of the Truncated Normal Distribution.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The "parent" normal distribution has\n' );
  fprintf ( 1, '    mean =               %g\n', mu );
  fprintf ( 1, '    standard deviation = %g\n', sigma );
  fprintf ( 1, '  The parent distribution is truncated to\n' );
  fprintf ( 1, '  the interval [%g,%g]\n', a, b );

  variance = truncated_normal_ab_variance ( mu, sigma, a, b );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PDF variance = %g\n', variance );

  for i = 1 : sample_num
    [ x(i), seed ] = truncated_normal_ab_sample ( mu, sigma, a, b, seed );
  end

  variance = var ( x );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Sample size =     %d\n', sample_num );
  fprintf ( 1, '  Sample variance = %g\n', variance );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ clmParams, pdm_mouth] = Load_CLM_params_inner()
%LOAD_CLM_PARAMS_WILD Summary of this function goes here
%   Detailed explanation goes here
    clmParams.window_size = [19,19;];
    clmParams.numPatchIters = size(clmParams.window_size,1);

    % the PDM created from in the wild data
    pdmLoc = ['../models/hierarch_pdm/pdm_51_inner.mat'];

    load(pdmLoc);

    pdm_mouth = struct;
    pdm_mouth.M = double(M);
    pdm_mouth.E = double(E);
    pdm_mouth.V = double(V);

    % the default model parameters to use
    clmParams.regFactor = 2.5;
    clmParams.sigmaMeanShift = 1.75;
    clmParams.tikhonov_factor = 2.5;

    clmParams.startScale = 1;
    clmParams.num_RLMS_iter = 5;
    clmParams.fTol = 0.01;
    clmParams.useMultiScale = true;
    clmParams.use_multi_modal = 1;

end

author: TadasBaltrusaitis

repo: OpenFace


text:
function g = sigmoidGradient(z)
%SIGMOIDGRADIENT returns the gradient of the sigmoid function
%evaluated at z
%   g = SIGMOIDGRADIENT(z) computes the gradient of the sigmoid function
%   evaluated at z. This should work regardless if z is a matrix or a
%   vector. In particular, if z is a vector or matrix, you should return
%   the gradient for each element.

g = zeros(size(z));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the gradient of the sigmoid function evaluated at
%               each value of z (z can be a matrix, vector or scalar).

g = sigmoid(z) .* (1 - sigmoid(z));

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

end

author: vugsus

repo: coursera-machine-learning


text:
function out = ctranspose( in )
%CTRANSPOSE for nested cells
%
% (c) Thomas Kuestner
% ---------------------------------------------------------------------
out = in; % TODO

end

author: thomaskuestner


text:
 function mse = ComputeMultiscaleEntropy(data, m, r, maxTau, maxVecSize,cg_moment,cg_method,constant_r)

% mse = ComputeMultiscaleEntropy(data, m, r, max_tau,cg_moment)
%
% Overview
%	Calculates multiscale entropy on a vector of input data.
%
% Input
%	data        - data to analyze; vector of doubles
%   m           - pattern length; int
%   r           - radius of similarity (% of std); double
%   maxTau      - maximum number of coarse-grainings; int
%   maxVecSize  - optional, parameter to switch from SampEn to FastSempEn
%   cg_moment   - optional, moment used to coarse-grain the time series,
%                 by default use the mean 'mean'
%                 Options: 'mean', 'varaince'
%   cg_method   - [optional] - method use to generate the coarse-grain
%                 time series. By default use the original method proposed
%                 by Costa et al. ('fir' with mean).
%                 Options: 'fir' [Costa et al.], 'butter' [Valencia et al.]
%   constant_r  - [DEFAULT] 1 : use r as function of std of original time
%                 series in this implementation correspond to zscore only
%                 the original time series
%               - 0 : r as a function of scale factor (tau), i.e.,
%                 r(i)= r*std(ScaleData(i)), in this implementation it
%                 corresponds to zscore  each coarse-grain time siries
%
%
% Output
%   mse          - vector of [max_tau, 1] doubles
%
% Example
%   data = rand(1e4, 1); % generate random data
%   m = 2; % template length
%   r = 0.2; % radius of similarity
%   maxTau = 4; % calculate sample entropy over four coarse grainings
%   mse = multiscaleEntropy(data, m, r, maxTau, 'Fast','mean');
%
%	REPO:
%       https://github.com/cliffordlab/PhysioNet-Cardiovascular-Signal-Toolbox
%
% Reference(s)
%
% Copyright (C) 2017 Erik Reinertsen <er@gatech.edu>
% All rights reserved.
%
%
% 08-23-2017 Modyfied by Giulia Da Poian to be inclused in the HRV Toolbox
% for PhysioNet Cardiovascular Signal ToolboxToolbox
% Removed the possibility to use different types of entropy, only
% fastSampen method in this version
%
% 10-10-2017 Modyfied by Giulia Da Poian, use FastSampEn in series < 34000
% otherwise use traditional SampEn that is faster for long series
%
% 10-23-2017 Modyfied by Giulia Da Poian, using scales like in Costa's
% paper instead of Coarse-grain data using halving method.
%
% This software may be modified and distributed under the terms
% of the BSD license. See the LICENSE file in this repo for details.

if nargin<5 || isempty(maxVecSize)
    maxVecSize = 34000;
end
if nargin<6 || isempty(cg_moment)
    cg_moment = 'mean';
end
if nargin<7 || isempty(cg_method)
    cg_method = 'fir';
end
if nargin<8 || isempty(constant_r)
    constant_r = 0;
end

data = zscore(data);  % (introduced by GDP) normalization of the signal that
                      % replace the common practice of expressing the
                      % tolerance as r times the standard deviation

mse = NaN(maxTau, 1); % Initialize output vector

SampEnType = 'Maxim'; % Initialize default SampEn method

% Check data length, if < 34000 use Fast Implementation (introduced GDP)
if length(data) < maxVecSize
     SampEnType = 'Fast';
end

% Loop through each window

% Loop through each timescale
% Note: i_tau == 1 is the original time series
for i_tau = 1:maxTau

    scaledData = coarsegrain(data,i_tau,cg_moment,cg_method); % Changed by GDP

    if ~constant_r
        scaledData = zscore(scaledData);
    end

    switch SampEnType
        case 'Fast'
           mse(i_tau) = fastSampen(scaledData, m, r);
        otherwise
           mse(i_tau) = sampenMaxim(scaledData, m, r);
    end

end % end for loop

% REFERENCES

% Costa et al. "Multiscale entropy analysis of complex physiologic time
% series." Physical review letters 89.6 (2002): 068102

% Valencia et al. "Refined multiscale entropy: Application to 24-h holter
% recordings of heart period variability in healthy and aortic stenosis
% subjects." IEEE Transactions on Biomedical Engineering 56, no. 9 (2009).

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function [Eft, Covft, ljpyt, Eyt, Covyt] = gpia_jpred(gp_array, x, y, varargin)
%GPIA_JPRED  Prediction with Gaussian Process GP_IA solution.
%
%  Description
%    [EFT, COVFT] = GPIA_JPRED(GP_ARRAY, X, Y, XT, OPTIONS)
%    takes a cell array of GP structures together with matrix X of
%    training inputs and vector Y of training targets, and
%    evaluates the predictive distribution at test inputs XT with
%    parameters marginalized out with IA. Returns a posterior mean
%    EFT and covariance COVFT of latent variables.
%
%    [EFT, COVFT, JPYT, EYT, COVYT] = GPIA_JPRED(GP, X, Y, XT, 'yt', YT, ...)
%    returns also logarithm of the predictive joint density PYT of
%    the observations YT at test input locations XT with parameters
%    marginalized out with IA. This can be used for example in the
%    cross-validation. Here Y has to be vector. Returns also
%    posterior predictive mean EYT and covariance COVYT.
%
%    [EF, COVF, LJPY, EY, COVY] = GPIA_JPRED(GP, X, Y, OPTIONS)
%    evaluates the predictive distribution at training inputs X
%    and logarithm of the joint predictive density JPY of the training
%    observations Y.
%
%    OPTIONS is optional parameter-value pair
%      predcf - index vector telling which covariance functions are
%               used for prediction. Default is all (1:gpcfn). See
%               additional information below.
%      tstind - a vector/cell array defining, which rows of X belong
%               to which training block in *IC type sparse models.
%               Deafult is []. In case of PIC, a cell array
%               containing index vectors specifying the blocking
%               structure for test data. IN FIC and CS+FIC a vector
%               of length n that points out the test inputs that
%               are also in the training set (if none, set TSTIND=[])
%      yt     - optional observed yt in test points (see below)
%      z      - optional observed quantity in triplet (x_i,y_i,z_i)
%               Some likelihoods may use this. For example, in case
%               of Poisson likelihood we have z_i=E_i, that is,
%               expected value for ith case.
%      zt     - optional observed quantity in triplet (xt_i,yt_i,zt_i)
%               Some likelihoods may use this. For example, in case
%               of Poisson likelihood we have z_i=E_i, that is, the
%               expected value for the ith case.
%
%    NOTE! In case of FIC and PIC sparse approximation the
%    prediction for only some PREDCF covariance functions is just
%    an approximation since the covariance functions are coupled in
%    the approximation and are not strictly speaking additive
%    anymore.
%
%    For example, if you use covariance such as K = K1 + K2 your
%    predictions Eft1 = gpia_pred(gp_array, X, Y, X, 'predcf', 1) and
%    Eft2 = gpia_pred(gp_array, x, y, x, 'predcf', 2) should sum up to
%    Eft = gpia_pred(gp_array, x, y, x). That is Eft = Eft1 + Eft2. With
%    FULL model this is true but with FIC and PIC this is true only
%    approximately. That is Eft \approx Eft1 + Eft2.
%
%    With CS+FIC the predictions are exact if the PREDCF covariance
%    functions are all in the FIC part or if they are CS
%    covariances.
%
%    NOTE! When making predictions with a subset of covariance
%    functions with FIC approximation the predictive variance can
%    in some cases be ill-behaved i.e. negative or unrealistically
%    small. This may happen because of the approximative nature of
%    the prediction.
%
%  See also
%    GP_PRED, GP_SET, GP_IA
%
% Copyright (c) 2009 Ville Pietil�inen
% Copyright (c) 2009-2010 Jarno Vanhatalo
% Copyright (c) 2011-2012 Ville Tolvanen
% Copyright (c) 2012 Aki Vehtari

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

  ip=inputParser;
  ip.FunctionName = 'GPIA_JPRED';
  ip.addRequired('gp_array', @iscell);
  ip.addRequired('x', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
  ip.addRequired('y', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
  ip.addOptional('xt', [], @(x) isempty(x) || (isreal(x) && all(isfinite(x(:)))))
  ip.addParamValue('yt', [], @(x) isreal(x) && all(isfinite(x(:))))
  ip.addParamValue('z', [], @(x) isreal(x) && all(isfinite(x(:))))
  ip.addParamValue('zt', [], @(x) isreal(x) && all(isfinite(x(:))))
  ip.addParamValue('predcf', [], @(x) isempty(x) || ...
                   isvector(x) && isreal(x) && all(isfinite(x)&x>0))
  ip.addParamValue('tstind', [], @(x) isempty(x) || iscell(x) ||...
                   (isvector(x) && isreal(x) && all(isfinite(x)&x>0)))
  if numel(varargin)==0 || isnumeric(varargin{1})
    % inputParser should handle this, but it doesn't
    ip.parse(gp_array, x, y, varargin{:});
  else
    ip.parse(gp_array, x, y, [], varargin{:});
  end
  xt=ip.Results.xt;
  yt=ip.Results.yt;
  z=ip.Results.z;
  zt=ip.Results.zt;
  predcf=ip.Results.predcf;
  tstind=ip.Results.tstind;
  if isempty(xt)
    xt=x;
    if isempty(tstind)
      if iscell(gp_array)
        gptype=gp_array{1}.type;
      else
        gptype=gp_array.type;
      end
      switch gptype
        case {'FULL' 'VAR' 'DTC' 'SOR'}
          tstind = [];
        case {'FIC' 'CS+FIC'}
          tstind = 1:size(x,1);
        case 'PIC'
          if iscell(gp)
            tstind = gp_array{1}.tr_index;
          else
            tstind = gp_array.tr_index;
          end
      end
    end
    if isempty(yt)
      yt=y;
    end
    if isempty(zt)
      zt=z;
    end
  end

  % pass these forward
  options=struct();
  if ~isempty(yt);options.yt=yt;end
  if ~isempty(z);options.z=z;end
  if ~isempty(zt);options.zt=zt;end
  if ~isempty(predcf);options.predcf=predcf;end
  if ~isempty(tstind);options.tstind=tstind;end

  if nargout > 2 && isempty(yt)
    pyt = NaN;
  end

  nGP = numel(gp_array);

  % Make predictions with different models in gp_array
  for i1=1:nGP
    P_TH(:,i1) = gp_array{i1}.ia_weight;
    if nargout <= 2
      [Efts(:,i1), Covfts(:,:,i1)]=gp_jpred(gp_array{i1},x,y,xt,options);
    else
      [Efts(:,i1), Covfts(:,:,i1), ljpyts(i1), Eyts(:,i1), Covyts(:,:,i1)]=gp_jpred(gp_array{i1},x,y,xt, options);
    end
  end

  % Calculate mean and variance of the distributions
  Eft = sum(bsxfun(@times,Efts,P_TH), 2);
  % Calculate covariances of means
  Efts = bsxfun(@minus,Efts,Eft);
  for i1=1:nGP
    CovEfts(:,:,i1)=Efts(:,i1)'*Efts(:,2);
  end
  % Covariance is E(Cov)+Cov(E)
  Covft = squeeze(sum(bsxfun(@times, Covfts, permute(P_TH,[1 3 2])), 3)) ...
          + sum(bsxfun(@times, CovEfts, permute(P_TH,[1 3 2])), 3);

  % Calculate jpyt with weight given in P_TH.
  if nargout > 2
    if isempty(yt)
      ljpyt = [];
    else
      ljpyt = log(sum(exp(ljpyts)'.*P_TH));
    end
  end

  if nargout > 3
    Eyt = sum(bsxfun(@times,Eyts,P_TH),2);
    Covyt = squeeze(sum(bsxfun(@times, Covyts, permute(P_TH,[1 3 2])), 3)) ...
            + diag(sum(bsxfun(@times,bsxfun(@minus,Eyts,Eyt).^2, P_TH),2));
  end

author: gpstuff-dev

repo: gpstuff

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


text:
function varargout=sort(varargin)
%SORT (overloaded)
%
% [t,loc] = sort(x)
%
% The variable t will be the sorted version of x.
%
% SORT is implemented in the nonlinear operator framework using a big-M
% model.

switch class(varargin{1})

    case 'double'
        error('Overloaded SDPVAR/SORT CALLED WITH DOUBLE. Report error')

    case 'sdpvar' % Overloaded operator for SDPVAR objects. Pass on args and save them.

        if 0
            error('')
        end

        y = [];

        update = varargin{1};
        for j = 1:length(update)
           y = [y;yalmip('define',mfilename,varargin{:},j)];
        end
        varargout{1} = y;

    case 'char' % YALMIP send 'graph' when it wants the epigraph or hypograph
        switch varargin{1}
            case {'milp','graph'}
                % Description using epigraphs
                t = varargin{2};
                X = varargin{3};

                % Call external to allow subsrefs in classs
                [F,vars] = pwadynamics_internal(t,varargin{3:end-1});

                varargout{1} = F;
                varargout{2} = struct('convexity','milp','monotonicity','milp','definiteness','milp');
                varargout{3} = X;

                % Inofficial way to model several nonlinear variables in
                % one call
                varargout{2}.models = vars;

            otherwise
                error('SDPVAR/SORT called with CHAR argument?');
        end
    otherwise
        error('Strange type on first argument in SDPVAR/SORT');
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function O = gs_orthonorm(M)
% O = gs_orthonorm(M)
%
% Classic Gram-Schmidt orthonormalization.
%
%                         PARAMETERS
%
% M A matrix with columns to be orthonormalized.
%
%                        RETURN VALUES
%
% O A matrix with orthonormal columns spanning the same
%   space as the columns in M. Linearly dependent vectors
%   are removed from the result.
%
% Author: Luca Vezzaro (elvezzaro@gmail.com)

	ncol = size(M, 2);

	% Zero threshold
	tol = 100 * eps;

	% Index to result column
	k = 1;

	for i=1:ncol
		% Column to orthogonalize
		v = M(:,i);

		% Subtract projections over previous vectors
		for j=1:k-1
			v = v - (O(:,j)' * v) * O(:,j);
		end

		% Only keep nonzero vectors
		n = norm(v);
		if n > tol
			% Normalize and store
			O(:,k) = v / n;
			k = k + 1;
		end
	end


text:
function [connmat] = triangle2connectivity(tri, pos)

% TRIANGLE2CONNECTIVITY computes a connectivity-matrix from a triangulation.
%
% Use as
%  [connmat] = triangle2connectivity(tri)
% or
%  [connmat] = triangle2connectivity(tri, pos)
%
% The input tri is an Mx3 matrix describing a triangulated surface,
% containing indices to connecting vertices. The output connmat is a sparse
% logical NxN matrix, with ones, where vertices are connected, and zeros
% otherwise.
%
% If you specify the vertex positions in the second input argument as Nx3
% matrix, the output will be a sparse matrix with the lengths of the
% edges between the connected vertices.
%
% See also CHANNELCONNECTIVIY

% Copyright (C) 2015, Jan-Mathijs Schoffelen
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% ensure that the vertices are indexed starting from 1
if min(tri(:))==0,
  tri = tri + 1;
end

% ensure that the vertices are indexed according to 1:number of unique vertices
tri = tri_reindex(tri);

% create the unique edges from the triangulation
edges  = [tri(:,[1 2]); tri(:,[1 3]); tri(:,[2 3])];
edges  = double(unique(sort([edges; edges(:,[2 1])],2), 'rows'));

% fill the connectivity matrix
n        = size(edges,1);
if nargin<2
  % create sparse binary matrix
  connmat = sparse([edges(:,1);edges(:,2)],[edges(:,2);edges(:,1)],true(2*n,1));
else
  % create sparse matrix with edge lengths
  dpos    = sqrt(sum( (pos(edges(:,1),:) - pos(edges(:,2),:)).^2, 2));
  connmat = sparse([edges(:,1);edges(:,2)],[edges(:,2);edges(:,1)],[dpos(:);dpos(:)]);
end

function [newtri] = tri_reindex(tri)

% this subfunction reindexes tri such that they run from 1:number of unique vertices
newtri       = tri;
[srt, indx]  = sort(tri(:));
tmp          = cumsum(double(diff([0;srt])>0));
newtri(indx) = tmp;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function y = permute( x, order )

%   Disciplined convex/geometric programming information for PERMUTE:
%      PERMUTE(A,ORDER) imposes no convexity restrictions on A. ORDER
%      must be constant.

%
% Determine the permutation
%

s = x.size_;
ndxs = 1 : prod( s );
ndx2 = permute( reshape( ndxs, s ), order );

%
% Permute the data
%

b = x.basis_;
try
    b = x.basis_( :, ndx2 );
catch
    ndxs( ndx2( : ).' ) = ndxs;
    [ r, c, v ] = find( b );
    b = sparse( r, ndxs( c ), v, size( b, 1 ), size( b, 2 ) );
    clear r c v
end
y = cvx( size( ndx2 ), b );

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function [nlcon,nlrhs,nle] = detNlcon(nonlcon,x0)
%DETNLCON  Determine OPTI format Nonlinear Constraints

%   Copyright (C) 2011 Jonathan Currie (IPL)

%Sort out nonlinear constraints (unfortunately rather inefficient)
if(isempty(nonlcon))
    nlcon = []; nlrhs = []; nle = [];
    return;
end

[in,eq] = nonlcon(x0);
nineq = length(in);
neq = length(eq);

nlcon = @(x) conval(nonlcon,x);
nlrhs = zeros(nineq+neq,1);
nle = [-1*ones(nineq,1);zeros(neq,1)];

%Dummy function
function c = conval(nonlcon,x)
[in,eq] = nonlcon(x);
c = [in;eq];

author: Arrowstar

repo: ksptot


text:
function b = r83_to_r8ge ( n, a )

%*****************************************************************************80
%
%% R83_TO_R8GE copies an R83 matrix to a R8GE matrix.
%
%  Discussion:
%
%    The R83 storage format is used for a tridiagonal matrix.
%    The superdiagonal is stored in entries (1,2:N), the diagonal in
%    entries (2,1:N), and the subdiagonal in (3,1:N-1).  Thus, the
%    original matrix is "collapsed" vertically into the array.
%
%  Example:
%
%    Here is how a R83 matrix of order 5 would be stored:
%
%       *  A12 A23 A34 A45
%      A11 A22 A33 A44 A55
%      A21 A32 A43 A54  *
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    N must be at least 2.
%
%    Input, real A(3,N), the R83 matrix.
%
%    Output, real B(N,N), the R8GE matrix.
%
  for i = 1 : n
    for j = 1 : n

      if ( j == i-1 )
        b(i,j) = a(3,j);
      elseif ( j == i )
        b(i,j) = a(2,j);
      elseif ( j == i+1 )
        b(i,j) = a(1,j);
      else
        b(i,j) = 0.0;
      end

    end
  end

  return
end


text:
function t = transmission_exponential(d, beta, varargin)
%TRANSMISSION_EXPONENTIAL  Compute transmission map using given depth map based
%on the Lambert-Beer law
%   Inputs:
%       -|d|: H-by-W matrix with values of depth for processed image in meters.
%       -|beta|: positive constant that parameterizes the density of haze.
%       Larger values indicate thicker haze.
%
%   Outputs:
%       -|t|: H-by-W matrix with medium transmission values ranging in [0, 1].

t = exp(-beta * d);

end

author: sakaridis