text:
function [C, Cinv] = gp_dcov(gp, x1, x2, predcf)
%GP_COV  Evaluate covariance matrix between two input vectors.
%
%  Description
%    C = GP_DCOV(GP, X1, X2, PREDCF) takes in Gaussian process GP and
%    two matrixes X1 and X2 that contain input vectors to GP.
%    Returns covariance matrix C between the latent values in the training
%    inputs X1 & XV and the test inputs X2. XV denotes the virtual input we
%    have (GP.XV). The covariance is computed between latent vector [f
%    df/dx_1 df/dx_2, ... df/dx_d], evaluated at X1 for f and XV for df/dx
%    respectively, and [f df/dx_1 df/dx_2, ... df/dx_d] evaluated at X2
%    for both the f and its gradient.
%
% Copyright (c) 2007-2010 Jarno Vanhatalo
% Copyright (c) 2010 Tuomas Nikoskinen
% Copyright (c) 2014 Ville Tolvanen

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

% Split the training data for normal latent input and gradient inputs
x12=x1;
x11=gp.xv;

% Derivative observations
[n,m]=size(x1);
[n4,m4]=size(x2);
ncf=length(gp.cf);
if isfield(gp, 'nvd')
  % Only specific dimensions
  ii1=abs(gp.nvd);
else
  % All dimensions
  ii1=1:m;
end
for i1=1:ncf
  gpcf = gp.cf{i1};
  if m==1
    Gset1 = gpcf.fh.ginput4(gpcf, x11, x2);
    Gset2 = gpcf.fh.ginput4(gpcf, x2, x12);
    Kff = gpcf.fh.cov(gpcf, x12, x2);
    Kdd = gpcf.fh.ginput2(gpcf, x11, x2);

    Kdf=Gset1{1};
    Kfd=Gset2{1};
    C = [Kff Kfd'; Kdf Kdd{1}];

    % Input dimension is >1
  else
    [n,m]=size(x11);
    [n2,m2]=size(x2);

    Kff = gpcf.fh.cov(gpcf, x12, x2);
    Gset1 = gpcf.fh.ginput4(gpcf, x11,x2);
    Gset2 = gpcf.fh.ginput4(gpcf, x2, x12);

    %Gather matrices from Gset (d k(x1,x2) /d x1)
    Kfd=cat(2,Gset1{ii1});
    Kdf=cat(1,Gset1{ii1});
    Kfd22=cat(2,Gset2{ii1});
    Kdf22=cat(1,Gset2{ii1})';

    % both x derivatives, same dimension (to diagonal blocks)
    D = gpcf.fh.ginput2(gpcf, x11, x2);
    % both x derivatives, different dimension (non-diagonal blocks)
    Kdf2 = gpcf.fh.ginput3(gpcf, x11 ,x2);

    % Now build up Kdd m*n x m*n2 matrix, which contains all the
    % both partial derivative" -matrices

    % Add the diagonal matrices
    Kdd=blkdiag(D{1:m});
    % Add the non-diagonal matrices to Kdd
    ii3=0;
    for j=0:m-2
      for i=1+j:m-1
        ii3=ii3+1;
        Kdd(i*n+1:(i+1)*n,j*n2+1:j*n2+n2) = Kdf2{ii3};
        Kdd(j*n+1:j*n+n,i*n2+1:(i+1)*n2) = Kdf2{ii3};
      end
    end
    if isfield(gp, 'nvd')
      % Collect the correct gradient dimensions, i.e. select the blocks
      % that correspond to the input dimensions for which we want the
      % gradients to be monotonic
      Kddtmp=[];
      for ii2=1:length(ii1)
        for ii3=ii2:length(ii1)
          Kddtmp((ii2-1)*n+1:ii2*n, (ii3-1)*n2+1:ii3*n2) = ...
            Kdd((ii1(ii2)-1)*n+1:ii1(ii2)*n,(ii1(ii3)-1)*n2+1:ii1(ii3)*n2);
          if ii2~=ii3
            Kddtmp((ii3-1)*n+1:ii3*n, (ii2-1)*n2+1:ii2*n2) = ...
              Kdd((ii1(ii3)-1)*n+1:ii1(ii3)*n,(ii1(ii2)-1)*n2+1:ii1(ii2)*n2);
          end
        end
      end
      Kdd=Kddtmp;
    end

    % Gather all the matrices into one final matrix K which is the
    % training covariance matrix
    C = [Kff Kdf22; Kdf Kdd];
    %   C = [Kff; Kdf];
  end
  if i1==1
    CC=C;
  else
    CC=CC+C;
  end
end
C=CC;

author: gpstuff-dev

repo: gpstuff

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


text:
function PlotFittedModel(protocol, MODELNAME, fittedpars, constants, h, style)
% Plots the normalized measurements predicted by the model with
% the fitted parameters against the absolute dot product of the
% gradient and fibre directions.
%
% PlotFittedModel(protocol, MODELNAME, fittedpars, h)
% adds the plot to figure handle h.
%
% protocol is the imaging protocol.
%
% MODELNAME is the name identifying the model.
%
% fittedpars and the model parameter values.
%
% constants is a structure containing fixed values required for the model.
%
% h is a figure handle to add the plot to.  If not specified, a new figure
% appears.
%
% author: Daniel C Alexander (d.alexander@ucl.ac.uk)
%         Gary Hui Zhang (gary.zhang@ucl.ac.uk)
%

if(nargin<5)
    h = figure;
end

% Get estimated fibre direction and b=0 signal.
fibredir = GetFibreOrientation(MODELNAME, fittedpars);
b0est = GetB0(MODELNAME, fittedpars);

if(nargin<6)
    style = '-';
end

linedef{1} = ['r', style];
linedef{2} = ['b', style];
linedef{3} = ['g', style];
linedef{4} = ['m', style];
linedef{5} = ['c', style];
linedef{6} = ['k', style];
linedef{7} = ['y', style];

hold on;

scale = GetScalingFactors(MODELNAME);
xsc = fittedpars(1:(length(scale)-1))./scale(1:(end-1));
S_Meas = SynthMeas(MODELNAME, xsc, protocol, fibredir, constants);
Snormdw = S_Meas./b0est;
for j=1:length(protocol.uG)
    inds = find(protocol.G == protocol.uG(j) & protocol.delta == protocol.udelta(j) & protocol.smalldel == protocol.usmalldel(j));
    dps = abs(protocol.grad_dirs(inds,:)*fibredir);
    [t tinds] = sort(dps);
    plot(dps(tinds), Snormdw(inds(tinds)), linedef{j}, 'LineWidth', 2);
end
xlabel('|n.G|/|G|_{max}');
ylabel('S/S_0');

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [data,times] = fiff_read_raw_segment(raw,from,to,sel)
%
% [data,times] = fiff_read_raw_segment(raw,from,to,sel)
%
% Read a specific raw data segment
%
% raw    - structure returned by fiff_setup_read_raw
% from   - first sample to include. If omitted, defaults to the
%          first sample in data
% to     - last sample to include. If omitted, defaults to the last
%          sample in data
% sel    - optional channel selection vector
%
% data   - returns the data matrix (channels x samples)
% times  - returns the time values corresponding to the samples (optional)
%

%
%   Author : Matti Hamalainen, MGH Martinos Center
%   License : BSD 3-clause
%
%   Revision 1.9  2007/11/22 05:04:25  msh
%   Fixed help text and some error checking
%
%   Revision 1.8  2006/04/23 15:29:40  msh
%   Added MGH to the copyright
%
%   Revision 1.7  2006/04/21 22:11:34  msh
%   Times calculation needed casting to double
%
%   Revision 1.6  2006/04/21 22:03:43  msh
%   Fixed time limit output.
%
%   Revision 1.5  2006/04/21 16:17:48  msh
%   Added handling of CTF compensation
%
%   Revision 1.4  2006/04/21 15:11:56  msh
%   Made the debug output optional.
%
%   Revision 1.3  2006/04/21 14:43:59  msh
%   Report projection operators.
%   Some more checks in raw data reading.
%
%   Revision 1.2  2006/04/21 14:23:16  msh
%   Further improvements in raw data reading
%
%

me='MNE:fiff_read_raw_segment';

if nargin == 3
    sel = [];
elseif nargin == 2
    to  = raw.last_samp;
    sel = [];
elseif nargin == 1
    from = raw.first_samp;
    to   = raw.last_samp;
    sel  = [];
elseif nargin ~= 4
    error(me,'Incorrect number of arguments');
end
%
%  Initial checks
%
from = double(from);
to   = double(to);
if from < raw.first_samp
    from = raw.first_samp;
end
if to > raw.last_samp
    to = raw.last_samp;
end
%
if from > to
    error(me,'No data in this range');
end
fprintf(1,'Reading %d ... %d  =  %9.3f ... %9.3f secs...', ...
    from,to,from/raw.info.sfreq,to/raw.info.sfreq);
%
%  Initialize the data and calibration vector
%
nchan = raw.info.nchan;
dest  = 1;
cal   = diag(raw.cals);
%
if isempty(sel)
    data = zeros(nchan,to-from+1);
    if isempty(raw.proj) && isempty(raw.comp)
        mult = [];
    else
        if isempty(raw.proj)
            mult = raw.comp*cal;
        elseif isempty(raw.comp)
            mult = raw.proj*cal;
        else
            mult = raw.proj*raw.comp*cal;
        end
    end
else
    data = zeros(length(sel),to-from+1);
    if isempty(raw.proj) && isempty(raw.comp)
        mult = [];
        cal  = diag(raw.cals(sel));
    else
        if isempty(raw.proj)
            mult = raw.comp(sel,:)*cal;
        elseif isempty(raw.comp)
            mult = raw.proj(sel,:)*cal;
        else
            mult = raw.proj(sel,:)*raw.comp*cal;
        end
    end
end
do_debug=false;
if ~isempty(cal)
    cal = sparse(cal);
end
if ~isempty(mult)
    mult = sparse(mult);
end

if isempty(fopen(raw.fid))
   fid = fopen(raw.info.filename,'rb','ieee-be');
   if (fid < 0)
      error(me,'Cannot open file %s',raw.info.filename);
   end
else
   fid = raw.fid;
end

for k = 1:length(raw.rawdir)
    this = raw.rawdir(k);
    %
    %  Do we need this buffer
    %
    if this.last > from
        if isempty(this.ent)
            %
            %  Take the easy route: skip is translated to zeros
            %
            if do_debug
                fprintf(1,'S');
            end
            doing_whole = false;
            if isempty(sel)
                one = zeros(nchan,this.nsamp);
            else
                one = zeros(length(sel),this.nsamp);
            end
        else
            tag = fiff_read_tag(fid,this.ent.pos);
            %
            %   Depending on the state of the projection and selection
            %   we proceed a little bit differently
            %
            if isempty(mult)
                if isempty(sel)
                    one = cal*double(reshape(tag.data,nchan,this.nsamp));
                else
                    one = double(reshape(tag.data,nchan,this.nsamp));
                    one = cal*one(sel,:);
                end
            else
                one = mult*double(reshape(tag.data,nchan,this.nsamp));
            end
        end
        %
        %  The picking logic is a bit complicated
        %
        if to >= this.last && from <= this.first
            %
            %    We need the whole buffer
            %
            first_pick = 1;
            last_pick  = this.nsamp;
            if do_debug
                fprintf(1,'W');
            end
        elseif from > this.first
            first_pick = from - this.first + 1;
            if to < this.last
                %
                %   Something from the middle
                %
                last_pick = this.nsamp + to - this.last;
                if do_debug
                    fprintf(1,'M');
                end
            else
                %
                %   From the middle to the end
                %
                last_pick = this.nsamp;
                if do_debug
                    fprintf(1,'E');
                end
            end
        else
            %
            %    From the beginning to the middle
            %
            first_pick = 1;
            last_pick  = to - this.first + 1;
            if do_debug
                fprintf(1,'B');
            end
        end
        %
        %   Now we are ready to pick
        %
        picksamp = last_pick - first_pick + 1;
        if picksamp >= 0
            data(:,dest:dest+picksamp-1) = one(:,first_pick:last_pick);
            dest = dest + picksamp;
        end
    end
    %
    %   Done?
    %
    if this.last >= to
        fprintf(1,' [done]\n');
        break;
    end
end

fclose(fid);

if nargout == 2
    times = [ from:to ];
    times = double(times)/raw.info.sfreq;
end

return;

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%% (Internal) Create a colormap for signal visualization
%
%   colors = my_colormap( cant_colors )
%
% Arguments:
%
%      + cant_colors: amount of colors required
%
% Output:
%
%      + colors: the colormap created
%
% Example:
%
% Author: Mariano Llamedo Soria llamedom@electron.frba.utn.edu.ar
% Version: 0.1 beta
% Last update: 14/5/2014
% Birthdate  : 21/4/2015
% Copyright 2008-2015
%
function colors = my_colormap( cant_colors )
colors = hsv(64);

%avoid yellowish colours
bAux = colors(:,1) > 200/255 & colors(:,2) > 200/255 & colors(:,3) == 0;
if( any(bAux) )
    colors( bAux,: ) = repmat([200 200 0]./255, sum(bAux),1 );
end

ncol = size(colors,1);

aux_idx = 1+rem(round(round(ncol/6)+linspace(0, (round(ncol/2)+1)*(ncol-1), ncol)), ncol);
if( cant_colors > 4 )
    colors = colors( aux_idx(round(linspace(1,64, cant_colors))) , :);
else
    colors = colors( aux_idx(1:cant_colors) , :);
end

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function [val,vec,info] = syevd(jobz,uplo,a)
%SYEV   Computing eigen information.
%   [VAL,VEC,INFO] = SYEVD(JOBZ,UPLO,A) computes eigenvalues and,
%	optionally, eigenvectors of a symmetric matrix A.
%
%	JOBZ allows the user to choose if eigenvectors are to be computed.
%	JOBZ
%          = 'N':  Compute eigenvalues only;
%          = 'V':  Compute eigenvalues and eigenvectors.
%
%	UPLO allows the user to choose which part of the matrix will be referenced.
%	UPLO
%          = 'U':  Upper triangle of A is stored;
%          = 'L':  Lower triangle of A is stored.
%
%   If only eigenvalues are desired, the QR algorithm is used.
%
%	No input parameter are optional. If you want an easy to use interface,
%	see SYEV_DRIVER.

    vec = [];
	val = [];
	info = [];

	if (nargin~=3),
		disp('Wrong number of input parameters.');
		return;
	end;

	    %Do all the boring checking
	if (~isnumeric(a)),
		disp('The matrix is not composed of numeric values.');
		return;
	elseif (isinteger(a)),
		a=double(a);
	end;

    if (size(a,1)~=size(a,2)),
        disp('The matrix must be square.');
        return;
	end;
	n=int32(size(a,1));
	lda=int32(n);

	if (isComplex(a)),
		if (isDouble(a))

			w=zeros(double(n),1);
			work=complex(zeros(1));
			rwork=zeros(1,1);
			iwork=int32(zeros(1,1));
			lrwork=int32(-1);
			liwork=int32(-1);
			lwork=int32(-1);
			info=int32(-1);
			[jobz, uplo, n, vec, lda, val, work, lwork, rwork,lrwork, iwork, liwork, info]=lapack_zheevd(jobz, uplo, n, a, lda, w, work, lwork, rwork,lrwork, iwork, liwork, info);
			lwork=int32(work(1));
			liwork=int32(iwork(1));
			lrwork=int32(rwork(1));

			work=complex(zeros(double(lwork),1));
			iwork=int32(zeros(double(lwork),1));
			rwork=zeros(double(lwork),1);
			[jobz, uplo, n, vec, lda, val, work, lwork, rwork,lrwork, iwork, liwork, info]=lapack_zheevd(jobz, uplo, n, a, lda, w, work, lwork, rwork,lrwork, iwork, liwork, info);
		else
			w=single(zeros(double(n),1));
			work=single(complex(zeros(1)));
			rwork=single(zeros(1,1));
			iwork=int32(zeros(1,1));
			lrwork=int32(-1);
			liwork=int32(-1);
			lwork=int32(-1);
			info=int32(-1);
			[jobz, uplo, n, vec, lda, val, work, lwork, rwork,lrwork, iwork, liwork, info]=lapack_cheevd(jobz, uplo, n, a, lda, w, work, lwork, rwork,lrwork, iwork, liwork, info);
			lwork=int32(work(1));
			liwork=int32(iwork(1));
			lrwork=int32(rwork(1));

			work=single(complex(zeros(double(lwork),1)));
			iwork=int32(zeros(double(liwork),1));
			rwork=single(zeros(double(lrwork),1));
			[jobz, uplo, n, vec, lda, val, work, lwork, rwork,lrwork, iwork, liwork, info]=lapack_cheevd(jobz, uplo, n, a, lda, w, work, lwork, rwork,lrwork, iwork, liwork, info);

		end;
	else
		if (isDouble(a)),
			w=zeros(double(n),1);
			work=zeros(1);
			iwork=int32(zeros(1));
			lwork=int32(-1);
			liwork=int32(-1);
			info=int32(-1);
			[jobz, uplo, n, vec, lda, val, work, lwork, iwork,liwork, info]= lapack_dsyevd(jobz, uplo, n, a, lda, w, work, lwork, iwork,liwork, info);
			lwork=int32(work(1));
			liwork=int32(iwork(1));
			work=zeros(double(lwork),1);
			iwork=int32(zeros(double(liwork),1));
			[jobz, uplo, n, vec, lda, val, work, lwork, iwork,liwork, info]= lapack_dsyevd(jobz, uplo, n, a, lda, w, work, lwork, iwork,liwork, info);
		else
			w=single(zeros(double(n),1));
			work=single(zeros(1));
			iwork=int32(zeros(1));
			lwork=int32(-1);
			liwork=int32(-1);
			info=int32(-1);
			[jobz, uplo, n, vec, lda, val, work, lwork, iwork,liwork, info]= lapack_ssyevd(jobz, uplo, n, a, lda, w, work, lwork, iwork,liwork, info);
			lwork=int32(work(1));
			liwork=int32(iwork(1));
			work=single(zeros(double(lwork),1));
			iwork=int32(zeros(double(liwork),1));
			[jobz, uplo, n, vec, lda, val, work, lwork, iwork,liwork, info]= lapack_ssyevd(jobz, uplo, n, a, lda, w, work, lwork, iwork,liwork, info);
		end;
	end;

repo: Aero-matlab


text:
%MAT_INTO_HDR  The old versions of SPM (any version before SPM5) store
%	an affine matrix of the SPM Reoriented image into a matlab file
%	(.mat extension). The file name of this SPM matlab file is the
%	same as the SPM Reoriented image file (.img/.hdr extension).
%
%	This program will convert the ANALYZE 7.5 SPM Reoriented image
%	file into NIfTI format, and integrate the affine matrix in the
%	SPM matlab file into its header file (.hdr extension).
%
%	WARNING: Before you run this program, please save the header
%	file (.hdr extension) into another file name or into another
%	folder location, because all header files (.hdr extension)
%	will be overwritten after they are converted into NIfTI
%	format.
%
%  Usage: mat_into_hdr(filename);
%
%  filename:	file name(s) with .hdr or .mat file extension, like:
%		'*.hdr', or '*.mat', or a single .hdr or .mat file.
%	e.g.	mat_into_hdr('T1.hdr')
%		mat_into_hdr('*.mat')
%

%  - Jimmy Shen (jimmy@rotman-baycrest.on.ca)
%
%-------------------------------------------------------------------------
function mat_into_hdr(files)

   pn = fileparts(files);
   file_lst = dir(files);
   file_lst = {file_lst.name};
   file1 = file_lst{1};
   [p n e]= fileparts(file1);

   for i=1:length(file_lst)
      [p n e]= fileparts(file_lst{i});
      disp(['working on file ', num2str(i) ,' of ', num2str(length(file_lst)), ': ', n,e]);
      process=1;

      if isequal(e,'.hdr')
         mat=fullfile(pn, [n,'.mat']);
         hdr=fullfile(pn, file_lst{i});

         if ~exist(mat,'file')
            warning(['Cannot find file "',mat  , '". File "', n, e, '" will not be processed.']);
            process=0;
         end
      elseif isequal(e,'.mat')
         hdr=fullfile(pn, [n,'.hdr']);
         mat=fullfile(pn, file_lst{i});

         if ~exist(hdr,'file')
            warning(['Can not find file "',hdr  , '". File "', n, e, '" will not be processed.']);
            process=0;
         end
      else
         warning(['Input file must have .mat or .hdr extension. File "', n, e, '" will not be processed.']);
         process=0;
      end

      if process
         load(mat);
         R=M(1:3,1:3);
         T=M(1:3,4);
         T=R*ones(3,1)+T;
         M(1:3,4)=T;

         [h filetype fileprefix machine]=load_nii_hdr(hdr);
         h.hist.qform_code=0;
         h.hist.sform_code=1;
         h.hist.srow_x=M(1,:);
         h.hist.srow_y=M(2,:);
         h.hist.srow_z=M(3,:);
         h.hist.magic='ni1';

         fid = fopen(hdr,'w',machine);
         save_nii_hdr(h,fid);
         fclose(fid);
      end
   end

   return;				% mat_into_hdr

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function specific_vortal_plots(up)

save_name = up.paths.filenames.hr_rr_scatter;
savepath = [up.paths.plots_save_folder, up.paths.filenames.hr_rr_scatter, '.eps'];
if ~up.analysis.redo_stats
    exist_log = check_exists(savepath, save_name);
    if exist_log
        return
    end
end

%% Make plot of feature vs filter resp sigs
if ~isempty(strfind(up.paths.paper_figures_folder, 'pc13')) & ...
        isempty(strfind(up.paths.data_save_folder, 'REC')) & isempty(strfind(up.paths.data_save_folder, 'WALK')) ...
        & isempty(strfind(up.paths.data_save_folder, 'EX'))
    % Load resp sigs
    subj = 2;
    load_path = [strrep(up.paths.data_save_folder, 'rest_and_rec', 'rest'), num2str(subj), up.paths.filenames.respSigs];
    resp_sigs = load(load_path);
    if sum(strcmp(fieldnames(resp_sigs), 'ppg_flt_Wam'))
        % Load raw sig
        load_path = [up.paths.data_save_folder, num2str(subj), up.paths.filenames.int_respSigs];
        int_resp_sigs = load(load_path);
        % Load annotated breaths
        up.paths.root_ann_folder = 'C:\Documents\Data\VORTAL\Manual_breath_annotations\';
        up.paths.observers = {'TB', 'DV'};
        up.paths.ann_folder = '_annotations\';
        obs_no = 1;
        data_folder = [up.paths.root_ann_folder, up.paths.observers{obs_no}, up.paths.ann_folder];
        % find SID
        temp = num2str(subj);
        if length(temp) == 1
            SID = ['00', temp];
        elseif length(temp) == 2
            SID = ['0', temp];
        end
        % Load annotated breaths
        rel_name = ['VORTAL' SID 'rest' up.paths.observers{obs_no} '-paw_an'];
        anns = load([data_folder, rel_name]); t = anns.PKS.t;
        % Load s ref time
        up.paths.s_ref_time_folder = 'C:\Documents\Data\VORTAL\Analysis_files\Processed_Data\';
        loadpath = [up.paths.s_ref_time_folder, SID, 's_con_rest'];
        load(loadpath, 'Sreftime', 'Send', 'Sstart');
        load(up.paths.db_data, 'db_data');
        bst_time = db_data.bst_log(db_data.sid == subj)/(24);
        % find t, measured in seconds since the start of the period
        t = 60*60*24*(t/(60*60*24) - floor( t/(60*60*24) ) - Sreftime + bst_time) - Sstart; % As the times are in seconds since 1970
        % Identify relevant section
        start_t = 92;
        rel_t = [start_t, start_t+up.paramSet.winLeng];
        rel_feat = int_resp_sigs.ppg_FMeam_FPt_PDtIMS_EHF;
        rel_filt = resp_sigs.ppg_flt_Wam;
        rel_sig = int_resp_sigs.ppg_EHF;
        rel_els = rel_feat.t >= rel_t(1) & rel_feat.t <= rel_t(2);
        rel_data.feat.t = rel_feat.t(rel_els);
        rel_data.feat.v = rel_feat.v(rel_els);
        rel_els = rel_filt.t >= rel_t(1) & rel_filt.t <= rel_t(2);
        rel_data.filt.t = rel_filt.t(rel_els);
        rel_data.filt.v = rel_filt.v(rel_els);
        rel_els = rel_sig.t >= rel_t(1) & rel_sig.t <= rel_t(2);
        rel_data.sig.t = rel_sig.t(rel_els);
        rel_data.sig.v = rel_sig.v(rel_els);

        % Make Figure
        h_fig = figure('Position', [200, 200, 600, 300]);

        % Plot
        ftsize = 12; lwidth = 2;
        h1 = plot(rel_data.sig.t-rel_t(1), detrend(rel_data.sig.v), 'b', 'LineWidth',1);
        hold on
        plot(rel_data.filt.t-rel_t(1), detrend(rel_data.filt.v), 'k', 'LineWidth',lwidth+2)
        h2 = plot(rel_data.filt.t-rel_t(1), detrend(rel_data.filt.v), 'c', 'LineWidth',lwidth);
        h3 = plot(rel_data.feat.t-rel_t(1), detrend(rel_data.feat.v), '-r', 'LineWidth',lwidth);
        h4 = plot(t-start_t, -0.05+zeros(length(t),1), '.k','MarkerSize',30);
        plot(rel_data.feat.t-rel_t(1), detrend(rel_data.feat.v), '.r', 'LineWidth',lwidth,'MarkerSize',20)
        xlim(rel_t-rel_t(1))
        temp = range(rel_data.sig.v);
        ylim([min(rel_data.sig.v)-0.05*temp, max(rel_data.sig.v)+0.05*temp])
        xlabel('Time [s]', 'FontSize', ftsize)
        ylabel('PPG', 'FontSize', ftsize)
        legend([h1, h2, h3, h4], {'PPG', 'Filter', 'Feature', 'Breaths'},'Location','northoutside','Orientation','horizontal')
        set(gca, 'FontSize', ftsize)
        set(gca, 'YTick', [])
        % Save
        set(gcf,'PaperUnits','inches');
        set(gcf,'PaperSize', [6, 3]);
        set(gcf,'PaperPosition',[0 0 6 3]);
        save_name = up.paths.filenames.feat_filt_plot;
        savepath = [up.paths.plots_save_folder, save_name];
        print(h_fig,'-depsc',savepath)
        close all
    end
end

if strfind(up.paths.data_save_folder, 'REST_AND_REC_TEMP')
    %% Create plot of improved precision with no of windows
    fprintf('\n--- Plot of improvement in precision with no of windows ');

    %% Load alg names
    load_name = up.paths.filenames.alg_names;
    loadpath = [up.paths.data_save_folder, up.paths.filenames.alg_names, '.mat'];
    load(loadpath, load_name);

    %% Load BA data
    load_name = 'BA_results';
    loadpath = [up.paths.data_save_folder, up.paths.filenames.global_BA, '.mat'];
    load(loadpath, load_name);
    rel_BA_res = BA_results.young;

    rel_els = find(strcmp(alg_names.sigs, 'ECG'));
    no_wins = alg_names.no_wins(rel_els);
    rel_prec = rel_BA_res.prec.val(rel_els);

    % Make Figure
    h_fig = figure('Position', [200, 200, 600, 300]);

    % Plot
    ftsize = 12; lwidth = 2;
    f2 = fit(no_wins,rel_prec,'exp2');
    coeffs = coeffvalues(f2);
    fitexp = ( coeffs(1)*exp(coeffs(2)*no_wins) ) + ( coeffs(3)*exp(coeffs(4)*no_wins) );
    plot(no_wins,rel_prec, 'mx', 'LineWidth',lwidth,'MarkerEdgeColor','k','MarkerFaceColor','k',...
        'MarkerSize',10)
    hold on
    plot(no_wins, fitexp, 'r', 'LineWidth', lwidth)
    ylim([8 12])
    xlim([1 max(no_wins)])
    xlabel('Number of windows used to estimate RR', 'FontSize', ftsize)
    ylabel('LOA interval [bpm]', 'FontSize', ftsize)
    set(gca, 'FontSize', ftsize)
    set(gca, 'YTick', 8:12)
    set(gcf,'PaperUnits','inches');
    set(gcf,'PaperSize', [6, 3]);
    set(gcf,'PaperPosition',[0 0 6 3]);
    save_name = up.paths.filenames.temp_prec_plot;
    savepath = [up.paths.plots_save_folder, save_name];
    print(h_fig,'-depsc',savepath)
    close all
end

fprintf('\n--- Making B-A Plots ');

%% Load data from entire study
load_name = up.paths.filenames.win_data;
loadpath = [up.paths.data_save_folder, up.paths.filenames.win_data, '.mat'];
load(loadpath, load_name);

%% Load algorithms names
load_name = up.paths.filenames.alg_names;
loadpath = [up.paths.data_save_folder, up.paths.filenames.alg_names, '.mat'];
load(loadpath, load_name);

%% Load relevant algorithm els
load_name = 'vortal_results';
loadpath = [up.paths.data_save_folder, up.paths.filenames.vortal_results, '.mat'];
load(loadpath, load_name);

%% Identify the most precise algorithm(s) to make plots for
no_plots_per_sig = 1;
for sig = {'ekg', 'ppg'}
    eval([sig{1,1}, '_els = vortal_results.top_ten.' sig{1,1}, '.alg_nos(1:no_plots_per_sig);']);
    eval([sig{1,1}, '_bias = vortal_results.top_ten.' sig{1,1}, '.bias.val(1:no_plots_per_sig);']);
    eval([sig{1,1}, '_lloa = vortal_results.top_ten.' sig{1,1}, '.lloa.val(1:no_plots_per_sig);']);
    eval([sig{1,1}, '_uloa = vortal_results.top_ten.' sig{1,1}, '.uloa.val(1:no_plots_per_sig);']);
end

rel_els = [ekg_els(:); ppg_els(:)];
rel_bias = [ekg_bias(:); ppg_bias(:)];
rel_lloa = [ekg_lloa(:); ppg_lloa(:)];
rel_uloa = [ekg_uloa(:); ppg_uloa(:)];

rel_names = alg_names.names(rel_els);
rel_sigs = alg_names.sigs(rel_els);

%% Generate plots for each algorithm
for rel_el_no = 1 : length(rel_els)

    alg_no = rel_els(rel_el_no);

    % identify data for this algorithm
    curr_els = win_data.alg_no == alg_no & win_data.snr_log;
    curr_est = win_data.est(curr_els);
    curr_ref = win_data.ref(curr_els);
    curr_bias = rel_bias(rel_el_no);
    curr_lloa = rel_lloa(rel_el_no);
    curr_uloa = rel_uloa(rel_el_no);
    curr_data.ave = nanmean([curr_est(:), curr_ref(:)],2);
    curr_data.error = [curr_est(:) - curr_ref(:)];

    % setup figure
    ftsize = 13;
    min_ave_val = 0;            % xlims
    max_ave_val = 36;
    min_error_val = -10;        % ylims
    max_error_val = 10;
    x_edge_int = 0.2;           % resolution of BA plot
    y_edge_int = 0.2;
    edges1 = [min_ave_val : x_edge_int : max_ave_val];       % locations of individual pixels
    edges2 = [min_error_val : y_edge_int : max_error_val];

    % Make color plot
    h_fig = figure('Position', [200 200 700 350]);
    A2=hist3([curr_data.ave,curr_data.error],'Edges',{edges1 edges2});
    h = fspecial('gaussian', [10 10], 2);
    y = filter2(h, A2');
    y = y./(max(max(y)));
    imagesc(y),
    c_han = colorbar;
    set(c_han,'YTickMode','manual')
    set(c_han, 'YTick',[0,1], 'YTickLabel', {'Min','Max'})

    % mymap = repmat(-1*[-1:0.01:0]', 1,3); colormap(mymap)
    hold on

    % X Ticks
    x_tick_int = 10/x_edge_int;
    x_label_int = (edges1(end)-edges1(1))*x_tick_int/length(edges1);
    xlims = xlim;
    edges1_labels = round( edges1(1) : x_label_int : edges1(end) )';

    x_spec.axis_lims = xlims;
    x_spec.real_lims = [edges1(1), edges1(end)];
    edges1_ticks = convert_real_to_axis_units(edges1_labels, x_spec);

    current_top_of_x_axis = xlims(1) + (0.5*(edges1(2)-edges1(1)));
    % adjust = - current_top_of_x_axis;
    set(gca, 'XTick', edges1_ticks, 'XTickLabel', num2str(edges1_labels))

    % Y Ticks
    y_tick_int = 10/y_edge_int;
    y_label_int = (edges2(end)-edges2(1))*y_tick_int/length(edges2);
    ylims = ylim;
    edges2_labels = round( edges2(1) : y_label_int : edges2(end) )';

    y_spec.axis_lims = ylims;
    y_spec.real_lims = [edges2(1), edges2(end)];
    edges2_ticks = convert_real_to_axis_units(edges2_labels, y_spec);

    current_top_of_y_axis = ylims(1) - (0.5*(edges2(2)-edges2(1)));
    %adjust = - current_top_of_y_axis;
    set(gca, 'YTick', edges2_ticks, 'YTickLabel', num2str(edges2_labels))

    % Axis Labels

    xlabel('Mean: 0.5(RR^{est} + RR^{ref}),  bpm', 'FontSize', ftsize);
    ylabel('Difference:  (RR^{est} - RR^{ref}),  bpm', 'FontSize', ftsize);

    % Bias and LOA Labels

    % xlim([min_ave_val, max_ave_val]), ylim([min_error_val, max_error_val])
    %xlims = xlim;
    %bias_adj = BA_data.bias + adjust;
    %loa_adj = [BA_data.lloa, BA_data.uloa] + adjust;
    axis_bias = convert_real_to_axis_units(curr_bias, y_spec);
    axis_lloa = convert_real_to_axis_units(curr_lloa, y_spec);
    axis_uloa = convert_real_to_axis_units(curr_uloa, y_spec);
    plot(xlims, [axis_bias, axis_bias], 'r')
    plot(xlims, [axis_lloa, axis_lloa], '--r')
    plot(xlims, [axis_uloa, axis_uloa], '--r')
    text(mean(xlims),0.5*axis_lloa,['Lower LOA = ' num2str(curr_lloa,2)],'HorizontalAlignment','center','BackgroundColor',[1 1 1], 'FontSize', ftsize);
    text(mean(xlims),1.2*axis_uloa,['Upper LOA = ' num2str(curr_uloa,2)],'HorizontalAlignment','center','BackgroundColor',[1 1 1], 'FontSize', ftsize);
    text(xlims(2),0.8*axis_bias,['Bias = ' num2str(curr_bias,2)],'HorizontalAlignment','right','BackgroundColor',[1 1 1], 'FontSize', ftsize);
    set(gca, 'FontSize', ftsize)

    set(gcf,'PaperUnits','inches');
    set(gcf,'PaperSize', [6, 3]);
    set(gcf,'PaperPosition',[0 0 6 3]);
    save_name = ['BA_' rel_sigs{rel_el_no}, rel_names{rel_el_no}];
    savepath = [up.paths.plots_save_folder, save_name];
    print(h_fig,'-depsc',savepath)
    close all

end

% %% CDF of abs errors plot
% fprintf('\n--- Making CDF Error Plot');
%
% % Load data (Best alg)
% load_name = up.paths.filenames.win_data;
% loadpath = [up.paths.data_save_folder, up.paths.filenames.win_data, '.mat'];
% load(loadpath, load_name);
% good_els = win_data.sqi & win_data.snr_log & win_data.alg_no == 203 & ~isnan(win_data.est);
% alg.est = win_data.est(good_els);
% alg.ref = win_data.ref(good_els);
% alg.error = alg.est - alg.ref;
%
% % Load data (Moderate alg)
% load_name = up.paths.filenames.win_data;
% loadpath = [up.paths.data_save_folder, up.paths.filenames.win_data, '.mat'];
% load(loadpath, load_name);
% good_els = win_data.sqi & win_data.snr_log & win_data.alg_no == 281 & ~isnan(win_data.est);
% mod_alg.est = win_data.est(good_els);
% mod_alg.ref = win_data.ref(good_els);
% mod_alg.error = mod_alg.est - mod_alg.ref;
%
% % Load data (Imp)
% load_name = up.paths.filenames.win_data;
% loadpath = [up.paths.data_save_folder, up.paths.filenames.win_data_imp, '.mat'];
% load(loadpath, load_name);
% good_els = win_data.sqi & win_data.snr_log & ~isnan(win_data.est) & ~isnan(win_data.est);
% imp.est = win_data.est(good_els);
% imp.ref = win_data.ref(good_els);
% imp.error = imp.est - imp.ref;
%
% % Normal data
% norm.error = normrnd(0,2.5,1,1000);
%
% % Make Figure
% ftsize = 12;
% h_fig = figure('Position', [200, 200, 1200, 550]);
% plot(sort(abs(alg.error)),linspace(0,1,length(alg.error))), hold on
% plot(sort(abs(imp.error)),linspace(0,1,length(imp.error)))
% plot(sort(abs(norm.error)),linspace(0,1,length(norm.error)))
% plot(sort(abs(mod_alg.error)),linspace(0,1,length(mod_alg.error)))
% title('Empirical CDF', 'FontSize', ftsize)
% xlabel('Absolute difference (bpm)', 'FontSize', ftsize)
% ylabel('F(x)', 'FontSize', ftsize)
% xlim([0 40])
% set(gca, 'XTick', [0:10, 15:5:40], 'YTick', [0:0.1:0.9, 0.95, 1.0], 'XTickLabel', {'0','','','','','5','','','','','10','15', '20', '25', '30', '35', '40'});
% legend({'a203', 'impedance', 'normal distribution', 'a281'}, 'Location', 'Best')
% grid on
%
% set(gca, 'FontSize', ftsize)
% save_name = 'Errors CDF';
% savepath = [up.paths.plots_save_folder, save_name];
% savefig(h_fig,savepath)
% close all

fprintf('\n--- Making HR / RR Scatter Plot');

%% Load data from entire study
load_name = up.paths.filenames.win_data;
loadpath = [up.paths.data_save_folder, up.paths.filenames.win_data, '.mat'];
load(loadpath, load_name);
hr = win_data.hr(win_data.comb_log & win_data.ecg_log & win_data.young_log & win_data.alg_no == 1);
rr = win_data.ref(win_data.comb_log & win_data.ecg_log & win_data.young_log & win_data.alg_no == 1);

% Make Figure
h_fig = figure('Position', [200, 200, 900, 400]);

% Plot Histograms
[counts,centers] = hist(rr,20); hold on
counts = counts/max(counts);
scale = 40;
bar(centers,scale*counts, 'c')
[counts,centers] = hist(hr,20);
counts = counts/max(counts);
scale = 5;
barh(centers,scale*counts, 'c')
rr_new = rr; hr_new = hr;

ftsize = 14;
plot(rr_new,hr_new, 'xk')
xlabel('Respiratory Rate [bpm]', 'FontSize', ftsize)
ylabel('Heart Rate [beats per minute]', 'FontSize', ftsize)
set(gca, 'FontSize', ftsize)
set(gca, 'YTick', 0:20:150)
set(gcf,'PaperUnits','inches');
set(gcf,'PaperSize', [9, 4]);
set(gcf,'PaperPosition',[0 0 9 4]);
save_name = up.paths.filenames.hr_rr_scatter;
savepath = [up.paths.plots_save_folder, save_name];
print(h_fig,'-depsc',savepath)
close all

%% Precision and prop of algs using filter / feature

fprintf('\n--- Making Prop of algs, filt/feat plot');

for plot_type = {'est_techs', 'sigs', 'comps'}
    % Load BA data for algs
    if ~strcmp(up.paths.root_data_folder, 'C:\Documents\Data\VORTAL_REST_AND_REC\')
        loadpath = [up.paths.data_save_folder, up.paths.filenames.global_BA, '.mat'];
        load(loadpath);
        rel_res = BA_results.young.prec;
        [~, orders] = sort(rel_res.val);
        rel_res.val = rel_res.val(orders); rel_res.val = rel_res.val(:)';
        rel_res.uci = rel_res.uci(orders); rel_res.uci = rel_res.uci(:)';
        rel_res.lci = rel_res.lci(orders); rel_res.lci = rel_res.lci(:)';
        if strcmp(plot_type{1,1}, 'est_techs')
            rel_els = alg_names.meths.ef(orders)>0;
            rel_res.pos_label = 'Frequency';
            rel_res.neg_label = 'Time';
        elseif strcmp(plot_type{1,1}, 'sigs')
            rel_els = strcmp(alg_names.sigs(orders), 'ECG');
            rel_res.pos_label = 'ECG';
            rel_res.neg_label = 'PPG';
        elseif strcmp(plot_type{1,1}, 'comps')
            rel_res.label1 = 'None';
            rel_res.label2 = 'X_{A1}';
            rel_res.label3 = 'X_{A4}';
            rel_res.label4 = 'E_{F3}';
            rel_els1 = alg_names.meths.xa~= 1 & alg_names.meths.xa~= 4 & alg_names.meths.ef~= 3;
            rel_els2 = alg_names.meths.xa== 1;
            rel_els3 = alg_names.meths.xa== 4;
            rel_els4 = alg_names.meths.ef== 3;
        end

    else
        loadpath = 'C:\Users\pc13\Dropbox\VORTAL\VORTAL_theoret_lims_yhvs\2016_Jan_Submission_to_Phys_Meas\Data\Complete_results_table.csv';
        [num, txt, raw] = xlsread(loadpath);
        headers = raw(1,:);
        % Eliminate IP
        rel_col = find(strcmp(headers, 'Signal'));
        bad_row = find(strcmp(raw(:,rel_col),'IP'));
        raw = raw([1:(bad_row-1), (bad_row+1):length(raw(:,1))], :);
        % Eliminate ones which had a problem with the random effects model (failed to converge)
        rel_col = find(strcmp(headers, 'Problems with Random Effects Model'));
        good_rows = ismember(raw(:,rel_col),'.'); good_rows(1) = 1;
        raw = raw(good_rows,:);
        % Extract relevant data
        rel_col = find(strcmp(headers, 'Overall Rank'));
        ranks = cell2mat(raw(2:end, rel_col));
        [~, rel_res.order] = sort(ranks);
        rel_col = find(strcmp(headers, 'Signal'));
        rel_res.sigs = raw(2:end, rel_col);
        rel_col = find(strcmp(headers, 'Algorithm'));
        rel_res.algs = raw(2:end, rel_col);
        rel_col = find(strcmp(headers, '2SD [bpm]'));
        rel_res.val = cell2mat(raw(2:end, rel_col));
        if strcmp(plot_type{1,1}, 'est_techs')
            rel_els = strfind(rel_res.algs, 'Ef');
            rel_els = not(cellfun('isempty', rel_els));
            rel_res.pos_label = 'Frequency';
            rel_res.neg_label = 'Time';
        elseif strcmp(plot_type{1,1}, 'sigs')
            rel_els = strcmp(rel_res.sigs(rel_res.order), 'ECG');
            rel_res.pos_label = 'ECG';
            rel_res.neg_label = 'PPG';
        elseif strcmp(plot_type{1,1}, 'comps')
            rel_res.label1 = 'None';
            rel_res.label2 = 'X_{A1}';
            rel_res.label3 = 'X_{A4}';
            rel_res.label4 = 'E_{F4}';
            temp = strfind(rel_res.algs, 'Xa1'); rel_els2 = ~cellfun(@isempty,temp);
            temp = strfind(rel_res.algs, 'Xa4'); rel_els3 = ~cellfun(@isempty,temp);
            temp = strfind(rel_res.algs, 'Ef4'); rel_els4 = ~cellfun(@isempty,temp);
            rel_els1 = ~rel_els2 & ~rel_els3 & ~rel_els4;
        end

    end
    orank = 1 : length(rel_res.val); orank = orank(:)';

    % data for hist
    tot = length(orank)+1;
    if strcmp(plot_type{1,1}, 'comps')
        bin_width = (tot/5);
    else
        bin_width = (tot/5);
    end
    bin_ends = 0:bin_width:tot;
    bin_starts = bin_ends(1:(end-1))+1;
    bin_ends = bin_ends(2:end);
    bin_mids = mean([bin_starts(:)'; bin_ends(:)']);
    clear prop*
    if ~strcmp(plot_type{1,1}, 'comps')
        for bin_no = 1 : length(bin_starts)
            prop1(bin_no) = 100*sum(orank(:)>= bin_starts(bin_no) & orank(:)< bin_ends(bin_no) & rel_els(:))/sum(orank(:)>= bin_starts(bin_no) & orank(:)< bin_ends(bin_no));
        end
        prop2 = 100 - prop1;
        y_vals = [prop1; prop2];
    else
        for bin_no = 1 : length(bin_starts)
            prop1(bin_no) = 100*sum(orank(:)>= bin_starts(bin_no) & orank(:)< bin_ends(bin_no) & rel_els1(:))/sum(orank(:)>= bin_starts(bin_no) & orank(:)< bin_ends(bin_no));
            prop2(bin_no) = 100*sum(orank(:)>= bin_starts(bin_no) & orank(:)< bin_ends(bin_no) & rel_els2(:))/sum(orank(:)>= bin_starts(bin_no) & orank(:)< bin_ends(bin_no));
            prop3(bin_no) = 100*sum(orank(:)>= bin_starts(bin_no) & orank(:)< bin_ends(bin_no) & rel_els3(:))/sum(orank(:)>= bin_starts(bin_no) & orank(:)< bin_ends(bin_no));
            prop4(bin_no) = 100*sum(orank(:)>= bin_starts(bin_no) & orank(:)< bin_ends(bin_no) & rel_els4(:))/sum(orank(:)>= bin_starts(bin_no) & orank(:)< bin_ends(bin_no));
        end
        y_vals = [prop1; prop2; prop3; prop4];
    end

    % setup fig
    h_fig = figure('Position', [200 200 900 550]);
    fontsize = 16;

    % upper plot
    xlims = [min(orank), max(orank)];
    if xlims(2) == 439
        xlims(2) = 440;
    end
    ylims = [0 50];
    h_prec = subplot(2,1,1, 'Position', [0.2 0.61 0.76 0.38]); hold on
    % fill([orank, fliplr(orank)], [rel_res.lci(:)', fliplr(rel_res.uci(:)')], 0.5*[1,1,1]);
    % hold on,
    plot(orank, rel_res.val, 'k', 'LineWidth', 3),
    ylab = ylabel({'2SD', '[bpm]'}, 'FontSize', fontsize, 'Rotation', 0);
    set(ylab, 'Units', 'Normalized', 'Position', [-0.17, 0.3, 0]);
    xlim(xlims)
    ylim(ylims)
    set(gca, 'XTick', [1, 50:50:tot, tot]);
    xlabel('Overall Rank', 'FontSize', fontsize)

    % lower plot (stacked bar)
    h_bar = subplot(2,1,2, 'Position', [0.2 0.10 0.76 0.38]); hold on
    h_bars = bar(bin_mids, y_vals', 'stacked');
    h_bars(1).FaceColor = 0.1*[1,1,1];
    h_bars(2).FaceColor = 0.5*[1,1,1];
    xlim(xlims),
    if ~strcmp(plot_type{1,1}, 'comps')
        ylim([0 100])
    else
        ylim([0 120])
    end
    if strcmp(plot_type{1,1}, 'est_techs')
        ylabel({'% using Frequency- and', 'Time-domain RR estimation'}, 'FontSize', fontsize)
    elseif strcmp(plot_type{1,1}, 'sigs')
        ylabel('% using ECG and PPG', 'FontSize', fontsize)
    elseif strcmp(plot_type{1,1}, 'comps')
        ylab = ylabel({'% using X_{A1},', 'X_{A4}, E_{F4}, or', 'none of these'}, 'FontSize', fontsize, 'Rotation', 0);
        set(ylab, 'Units', 'Normalized', 'Position', [-0.16, 0.33, 0]);
    end
    xlabel('Overall Rank', 'FontSize', fontsize)
    starts = strread(num2str(ceil(bin_starts)),'%s');
    mids = strread(num2str(ceil(bin_mids)),'%s');
    ends = strread(num2str(ceil(bin_ends)),'%s');
    quintiles = strread(num2str(1:length(bin_ends)),'%s');
    for s = 1 : length(ends)
        xtickstr{s} = [starts{s} ' - ' ends{s}];
        %xtickstr{s} = [quintiles{s}, ' (' starts{s} ' - ' ends{s} ')'];
    end
    set(h_bar, 'XTick', bin_mids, 'XTickLabel', xtickstr)

    allAxesInFigure = findall(h_fig,'type','axes');
    set(allAxesInFigure, 'Fontsize', fontsize-4, 'xgrid','on')

    % annotate
    if ~strcmp(plot_type{1,1}, 'comps')
        text(bin_mids(5),15, rel_res.pos_label, 'HorizontalAlignment', 'Center', 'FontSize', fontsize-4, 'Backgroundcolor', 'w')
        text(bin_mids(5),87, rel_res.neg_label, 'HorizontalAlignment', 'Center', 'FontSize', fontsize-4, 'Backgroundcolor', 'w')
    else
        text(bin_mids(1),50, rel_res.label1, 'HorizontalAlignment', 'Center', 'FontSize', fontsize-4, 'Backgroundcolor', 'w')
        text(bin_mids(5),22, rel_res.label2, 'HorizontalAlignment', 'Center', 'FontSize', fontsize-4, 'Backgroundcolor', 'w')
        text(bin_mids(5),62, rel_res.label3, 'HorizontalAlignment', 'Center', 'FontSize', fontsize-4, 'Backgroundcolor', 'w')
        text(bin_mids(5),95, rel_res.label4, 'HorizontalAlignment', 'Center', 'FontSize', fontsize-4, 'Backgroundcolor', 'w')
    end

    % add (a) and (b)
    dim = [.01 .91 .1 .1]; str = '(a)';
    annotation('textbox',dim,'String',str,'FitBoxToText','on','FontSize',fontsize,'LineStyle','none');
    dim = [.01 .43 .1 .1]; str = '(b)';
    annotation('textbox',dim,'String',str,'FitBoxToText','on','FontSize',fontsize,'LineStyle','none');

    set(gcf,'color','w');

    set(gcf,'PaperUnits','inches');
    set(gcf,'PaperSize', [9, 5.5]);
    set(gcf,'PaperPosition',[0 0 9 5.5]);
    save_name = [up.paths.filenames.stacked_bar, '_', plot_type{1,1}];
    savepath = [up.paths.plots_save_folder, save_name];
    print(h_fig,'-depsc',savepath)

    close all

end

%% New joint stacked bar plot

plot_type = {'comps_est_techs'};
% Load BA data for algs
if strcmp(up.paths.root_data_folder, 'C:\Documents\Data\VORTAL_REST_AND_REC\')
    loadpath = 'C:\Users\pc13\Dropbox\VORTAL\VORTAL_theoret_lims_yhvs\2016_Jan_Submission_to_Phys_Meas\Data\Complete_results_table.csv';
    [num, txt, raw] = xlsread(loadpath);
    headers = raw(1,:);
    % Eliminate IP
    rel_col = find(strcmp(headers, 'Signal'));
    bad_row = find(strcmp(raw(:,rel_col),'IP'));
    raw = raw([1:(bad_row-1), (bad_row+1):length(raw(:,1))], :);
    % Eliminate ones which had a problem with the random effects model (failed to converge)
    rel_col = find(strcmp(headers, 'Problems with Random Effects Model'));
    good_rows = ismember(raw(:,rel_col),'.'); good_rows(1) = 1;
    raw = raw(good_rows,:);
    % Extract relevant data
    rel_col = find(strcmp(headers, 'Overall Rank'));
    ranks = cell2mat(raw(2:end, rel_col));
    [~, rel_res.order] = sort(ranks);
    rel_col = find(strcmp(headers, 'Signal'));
    rel_res.sigs = raw(2:end, rel_col);
    rel_col = find(strcmp(headers, 'Algorithm'));
    rel_res.algs = raw(2:end, rel_col);
    rel_col = find(strcmp(headers, '2SD [bpm]'));
    rel_res.val = cell2mat(raw(2:end, rel_col));
    rel_els = strfind(rel_res.algs, 'Ef');
    rel_els = not(cellfun('isempty', rel_els));
    rel_res.pos_label = 'Frequency';
    rel_res.neg_label = 'Time';
    rel_res.label1 = 'None';
    rel_res.label2 = 'X_{A1}';
    rel_res.label3 = 'X_{A4}';
    rel_res.label4 = 'E_{F4}';
    temp = strfind(rel_res.algs, 'Xa1'); rel_els2 = ~cellfun(@isempty,temp);
    temp = strfind(rel_res.algs, 'Xa4'); rel_els3 = ~cellfun(@isempty,temp);
    temp = strfind(rel_res.algs, 'Ef4'); rel_els4 = ~cellfun(@isempty,temp);
    rel_els1 = ~rel_els2 & ~rel_els3 & ~rel_els4;
    orank = 1 : length(rel_res.val); orank = orank(:)';

    % Create first histogram
    tot = length(orank)+1;
    if strcmp(plot_type{1,1}, 'comps')
        bin_width = (tot/5);
    else
        bin_width = (tot/5);
    end
    bin_ends = 0:bin_width:tot;
    bin_starts = bin_ends(1:(end-1))+1;
    bin_ends = bin_ends(2:end);
    bin_mids = mean([bin_starts(:)'; bin_ends(:)']);
    clear prop*
    for bin_no = 1 : length(bin_starts)
        prop1(bin_no) = 100*sum(orank(:)>= bin_starts(bin_no) & orank(:)< bin_ends(bin_no) & rel_els1(:))/sum(orank(:)>= bin_starts(bin_no) & orank(:)< bin_ends(bin_no));
        prop2(bin_no) = 100*sum(orank(:)>= bin_starts(bin_no) & orank(:)< bin_ends(bin_no) & rel_els2(:))/sum(orank(:)>= bin_starts(bin_no) & orank(:)< bin_ends(bin_no));
        prop3(bin_no) = 100*sum(orank(:)>= bin_starts(bin_no) & orank(:)< bin_ends(bin_no) & rel_els3(:))/sum(orank(:)>= bin_starts(bin_no) & orank(:)< bin_ends(bin_no));
        prop4(bin_no) = 100*sum(orank(:)>= bin_starts(bin_no) & orank(:)< bin_ends(bin_no) & rel_els4(:))/sum(orank(:)>= bin_starts(bin_no) & orank(:)< bin_ends(bin_no));
    end
    y_vals = [prop1; prop2; prop3; prop4];

    % setup fig
    h_fig = figure('Position', [200 200 900 800]);
    fontsize = 16;

    % upper plot
    xlims = [min(orank), max(orank)];
    if xlims(2) == 439
        xlims(2) = 440;
    end
    ylims = [0 50];
    h_prec = subplot(3,1,1, 'Position', [0.2 0.74 0.76 0.25]); hold on
    plot(orank, rel_res.val, 'k', 'LineWidth', 3),
    ylab = ylabel({'2SD', '[bpm]'}, 'FontSize', fontsize, 'Rotation', 0);
    set(ylab, 'Units', 'Normalized', 'Position', [-0.17, 0.3, 0]);
    xlim(xlims)
    ylim(ylims)
    set(gca, 'XTick', [1, 50:50:tot, tot]);
    xlabel('Overall Rank', 'FontSize', fontsize)

    % mid plot (stacked bar)
    h_bar = subplot(3,1,2, 'Position', [0.2 0.41 0.76 0.25]); hold on
    h_bars = bar(bin_mids, y_vals', 'stacked');
    h_bars(1).FaceColor = 0.1*[1,1,1];
    h_bars(2).FaceColor = 0.5*[1,1,1];
    xlim(xlims),
    ylim([0 120])
    ylab = ylabel({'% using X_{A1},', 'X_{A4}, E_{F4}, or', 'none of these'}, 'FontSize', fontsize, 'Rotation', 0);
    set(ylab, 'Units', 'Normalized', 'Position', [-0.16, 0.33, 0]);
    xlabel('Overall Rank', 'FontSize', fontsize)
    starts = strread(num2str(ceil(bin_starts)),'%s');
    mids = strread(num2str(ceil(bin_mids)),'%s');
    ends = strread(num2str(ceil(bin_ends)),'%s');
    quintiles = strread(num2str(1:length(bin_ends)),'%s');
    for s = 1 : length(ends)
        xtickstr{s} = [starts{s} ' - ' ends{s}];
    end
    set(h_bar, 'XTick', bin_mids, 'XTickLabel', xtickstr)
    % annotate
    text(bin_mids(1),50, rel_res.label1, 'HorizontalAlignment', 'Center', 'FontSize', fontsize-4, 'Backgroundcolor', 'w')
    text(bin_mids(5),25, rel_res.label2, 'HorizontalAlignment', 'Center', 'FontSize', fontsize-4, 'Backgroundcolor', 'w')
    text(bin_mids(5),65, rel_res.label3, 'HorizontalAlignment', 'Center', 'FontSize', fontsize-4, 'Backgroundcolor', 'w')
    text(bin_mids(5),96, rel_res.label4, 'HorizontalAlignment', 'Center', 'FontSize', fontsize-4, 'Backgroundcolor', 'w')

    % lower plot (stacked bar)
    for bin_no = 1 : length(bin_starts)
        prop1(bin_no) = 100*sum(orank(:)>= bin_starts(bin_no) & orank(:)< bin_ends(bin_no) & rel_els(:))/sum(orank(:)>= bin_starts(bin_no) & orank(:)< bin_ends(bin_no));
    end
    prop2 = 100 - prop1;
    y_vals = [prop1; prop2];
    h_bar = subplot(3,1,3, 'Position', [0.2 0.08 0.76 0.25]); hold on
    h_bars = bar(bin_mids, y_vals', 'stacked');
    h_bars(1).FaceColor = 0.1*[1,1,1];
    h_bars(2).FaceColor = 0.5*[1,1,1];
    xlim(xlims),
    ylim([0 100])
    ylab = ylabel({'% using', 'Frequency- or', 'Time-domain', 'RR estimation'}, 'FontSize', fontsize, 'Rotation', 0);
    set(ylab, 'Units', 'Normalized', 'Position', [-0.16, 0.28, 0]);
    xlabel('Overall Rank', 'FontSize', fontsize)
    starts = strread(num2str(ceil(bin_starts)),'%s');
    mids = strread(num2str(ceil(bin_mids)),'%s');
    ends = strread(num2str(ceil(bin_ends)),'%s');
    quintiles = strread(num2str(1:length(bin_ends)),'%s');
    for s = 1 : length(ends)
        xtickstr{s} = [starts{s} ' - ' ends{s}];
    end
    set(h_bar, 'XTick', bin_mids, 'XTickLabel', xtickstr)
    % annotate
    text(bin_mids(5),15, rel_res.pos_label, 'HorizontalAlignment', 'Center', 'FontSize', fontsize-4, 'Backgroundcolor', 'w')
    text(bin_mids(5),87, rel_res.neg_label, 'HorizontalAlignment', 'Center', 'FontSize', fontsize-4, 'Backgroundcolor', 'w')

    % add (a) and (b) and (c)
    dim = [.01 .91 .1 .1]; str = '(a)';
    annotation('textbox',dim,'String',str,'FitBoxToText','on','FontSize',fontsize,'LineStyle','none');
    dim = [.01 .58 .1 .1]; str = '(b)';
    annotation('textbox',dim,'String',str,'FitBoxToText','on','FontSize',fontsize,'LineStyle','none');
    dim = [.01 .25 .1 .1]; str = '(c)';
    annotation('textbox',dim,'String',str,'FitBoxToText','on','FontSize',fontsize,'LineStyle','none');

    allAxesInFigure = findall(h_fig,'type','axes');
    set(allAxesInFigure, 'Fontsize', fontsize-4, 'xgrid','on')

    set(gcf,'color','w');
    set(gcf,'PaperUnits','inches');
    set(gcf,'PaperSize', [9, 8.0]);
    set(gcf,'PaperPosition',[0 0 9 8.0]);
    save_name = [up.paths.filenames.stacked_bar, '_', plot_type{1,1}];
    savepath = [up.paths.plots_save_folder, save_name];
    print(h_fig,'-depsc',savepath)
    close all

end

end

function axis_units = convert_real_to_axis_units(real_units, spec)

axis_units = spec.axis_lims(1) + ( (real_units - spec.real_lims(1)) * (spec.axis_lims(end) - spec.axis_lims(1))/ (spec.real_lims(end)-spec.real_lims(1)));

end

author: peterhcharlton

repo: RRest


text:
function ft_plot_montage(dat, varargin)

% FT_PLOT_MONTAGE makes a montage of a 3-D array by selecting slices at regular distances
% and combining them in one large 2-D image.  Note that the montage of MRI slices is not to
% be confused with the EEG montage, which is a way of specifying the reference scheme
% between electrodes.
%
% Use as
%   ft_plot_montage(dat, ...)
% where dat is a 3-D array.
%
% Additional options should be specified in key-value pairs and can be
%   'transform'     = 4x4 homogeneous transformation matrix specifying the mapping from voxel space to the coordinate system in which the data are plotted.
%   'location'      = 1x3 vector specifying a point on the plane which will be plotted, the coordinates are expressed in the coordinate system in which the data will be plotted. location defines the origin of the plane
%   'orientation'   = 1x3 vector specifying the direction orthogonal through the plane which will be plotted (default = [0 0 1])
%   'srange'        =
%   'slicesize'     =
%   'nslice'        = scalar, number of slices
%   'maskstyle'     = string, 'opacity' or 'colormix', defines the rendering
%   'background'    = needed when maskstyle is 'colormix', 3D-matrix with
%                     the same size as the data matrix, serving as
%                     grayscale image that provides the background
%
% See also FT_PLOT_ORTHO, FT_PLOT_SLICE, FT_SOURCEPLOT

% undocumented, these are passed on to FT_PLOT_SLICE
%   'intersectmesh'  = triangulated mesh through which the intersection of the plane will be plotted (e.g. cortical sheet)
%   'intersectcolor' = color for the intersection

% Copyrights (C) 2012, 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$

transform       = ft_getopt(varargin, 'transform', eye(4));
loc             = ft_getopt(varargin, 'location');
ori             = ft_getopt(varargin, 'orientation');
srange          = ft_getopt(varargin, 'slicerange');
slicesize       = ft_getopt(varargin, 'slicesize');
nslice          = ft_getopt(varargin, 'nslice');
backgroundcolor = ft_getopt(varargin, 'backgroundcolor', [0 0 0]);

% the intersectmesh and plotmarker options are passed on to FT_PLOT_SLICE
dointersect = ~isempty(ft_getopt(varargin, 'intersectmesh'));
domarker    = ~isempty(ft_getopt(varargin, 'plotmarker'));

% set the location if empty
if isempty(loc) && (isempty(transform) || isequal(transform, eye(4)))
  % go to the middle of the volume if the data seem to be in voxel coordinates
  loc = size(dat)./2;
elseif isempty(loc)
  % otherwise take the origin of the coordinate system
  loc = [0 0 0];
end

% check compatibility of inputs
if size(loc, 1) == 1 && isempty(nslice)
  nslice = 20;
elseif size(loc, 1) == 1 && ~isempty(nslice)
  % this is not a problem, slice spacing will be determined
elseif size(loc, 1) > 1 && isempty(nslice)
  % this is not a problem, number of slices is determined by loc
  nslice = size(loc, 1);
elseif size(loc, 1) > 1 && ~isempty(nslice)
  if size(loc, 1) ~= nslice
    ft_error('you should either specify a set of locations or a single location with a number of slices');
  end
end

% set the orientation if empty
if isempty(ori)
  ori = [0 0 1];
end

% ensure the ori to have unit norm
for k = 1:size(ori,1)
  ori(k,:) = ori(k,:)./norm(ori(k,:));
end

% determine the slice range
if size(loc, 1) == 1 && nslice > 1
  if isempty(srange) || (ischar(srange) && strcmp(srange, 'auto'))
    srange = [-50 70];
  else
  end
  loc = repmat(loc, [nslice 1]) + linspace(srange(1),srange(2),nslice)'*ori;
end

% ensure that the ori has the same size as the loc
if size(ori,1)==1 && size(loc,1)>1
  ori = repmat(ori, size(loc,1), 1);
end

div     = [ceil(sqrt(nslice)) ceil(sqrt(nslice))];
optarg  = varargin;
corners = [inf -inf inf -inf inf -inf]; % get the corners for the axis specification

for k = 1:nslice
  % define 'x' and 'y' axis in projection plane, the definition of x and y is more or less arbitrary
  [x, y] = projplane(ori(k,:)); % z = ori

  % get the transformation matrix to project onto the xy-plane
  T  = [x(:) y(:) ori(k,:)' loc(k,:)'; 0 0 0 1];

  optarg = ft_setopt(optarg, 'location',    loc(k,:));
  optarg = ft_setopt(optarg, 'orientation', ori(k,:));
  ix     = mod(k-1, div(1));
  iy     = floor((k-1)/div(1));
  h(k)   = ft_plot_slice(dat, optarg{:}); % FIXME is it safe to pass all optinoal inputs?

  xtmp = get(h(k), 'xdata');
  ytmp = get(h(k), 'ydata');
  ztmp = get(h(k), 'zdata');
  siz  = size(xtmp);
  if k==1 && isempty(slicesize)
    slicesize = siz;
  end

  % project the positions onto the xy-plane
  pos = [xtmp(:) ytmp(:) ztmp(:)];
  pos = ft_warp_apply(inv(T), pos);

  xtmp = reshape(pos(:,1), siz);
  ytmp = reshape(pos(:,2), siz);
  ztmp = reshape(pos(:,3), siz);

  % add some offset in the x and y directions to create the montage
  offset(1) = iy*(slicesize(1)-1);
  offset(2) = ix*(slicesize(2)-1);

  % update the specification of the corners of the montage plot
  if ~isempty(xtmp)
    c1 = offset(1) + min(xtmp(:));
    c2 = offset(1) + max(xtmp(:));
    c3 = offset(2) + min(ytmp(:));
    c4 = offset(2) + max(ytmp(:));
    c5 = min(ztmp(:));
    c6 = max(ztmp(:));
  end
  corners = [min(corners(1),c1) max(corners(2),c2) min(corners(3),c3) max(corners(4),c4) min(corners(5),c5) max(corners(6),c6)];

  % update the positions
  set(h(k), 'ydata', offset(1) + xtmp);
  set(h(k), 'xdata', offset(2) + ytmp);
  set(h(k), 'zdata',         0 * ztmp);

  if dointersect || domarker
    if ~exist('pprevious', 'var'), pprevious = []; end
    p = setdiff(findobj(gcf, 'type', 'patch'), pprevious);
    for kk = 1:numel(p)
      xtmp = get(p(kk), 'xdata');
      ytmp = get(p(kk), 'ydata');
      ztmp = get(p(kk), 'zdata');
      siz2 = size(xtmp);

      pos = [xtmp(:) ytmp(:) ztmp(:)];
      pos = ft_warp_apply(inv(T), pos);

      xtmp = reshape(pos(:,1), siz2);
      ytmp = reshape(pos(:,2), siz2);
      ztmp = reshape(pos(:,3), siz2);

      % update the positions
      set(p(kk), 'ydata', offset(1) + xtmp);
      set(p(kk), 'xdata', offset(2) + ytmp);
      set(p(kk), 'zdata',         0.0001 * ztmp);
    end
    pprevious = [pprevious(:);p(:)];
  end
  % drawnow; %this statement slows down the process big time on some file
  %systems. I don't know what's going on there, but the statement is not
  %really necessary, so commented out.
end
set(gcf, 'color', backgroundcolor);
set(gca, 'zlim', [0 1]);
%axis equal;
axis off;
view([0 90]);
axis(corners([3 4 1 2]));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [x, y] = projplane(z)
[u, s, v] = svd([eye(3) z(:)]);
x = u(:, 2)';
y = u(:, 3)';


% 3D Example for Multilevel Mass-Preserving Image Registration using VAMPIRE
%
% (c) Fabian Gigengack and Lars Ruthotto 2011/02/04, see FAIR.2 and FAIRcopyright.m.
% http://www.uni-muenster.de/EIMI/
% http://www.mic.uni-luebeck.de/
%
%   - data                 Cardiac gated PET images of a mouse heart.
%                          (level 4:6, full resolution: 40x40x40)
%                          Images show the heart in systole and diastole.
%   - viewer               viewImage3D
%   - interpolation        splineInterMex (regularizer='moments', theta=0.01)
%   - distance             SSD
%   - pre-registration     none
%   - regularizer          mfHyperElastic
%   - optimizer            Gauss-Newton with ArmijoBacktrack linesearch
%
% Acknowledgement:
% ----------------
% Thanks to the European Institute for Molecular Imaging (EIMI) and SFB 656,
% University of Muenster, Germany for supplying this interesting data.

setup3DmouseData

% prepare the plot
FAIRplots('clear')
Dshow = @(T,R,omega,m) viewIP(abs(T-R),omega,m,'colormap',gray(256));
FAIRplots('set','Dshow',Dshow);

% set the regularizer
alpha       = 1;
alphaLength = 1;
alphaArea   = 0.1;
alphaVolume = 2;
regularizer('reset', 'regularizer', 'mfHyperElastic', 'alpha', alpha, ...
    'alphaLength', alphaLength, ...
    'alphaArea',   alphaArea, ...
    'alphaVolume', alphaVolume);

% finally: run the MultiLevel Non-Parametric Image Registration
NPIRpara            = optPara('NPIR-GN');
NPIRpara.lineSearch = @ArmijoDiffeomorphic;
NPIRpara.solver     = @VAMPIREsolveGN_PCG;

[yc,wc,his] = MLIR(ML, 'NPIRobj', @VAMPIRENPIRobjFctn, ...
    'parametric', false, 'NPIRpara', NPIRpara);

%% Plot results
% Compute resulting image: dataT(yc) * det(D(yc))
Topt = reshape(linearInter(dataT,omega,center(yc,m)) .* geometry(yc,m,'Jac','omega',omega), m);

figure(1); clf;
subplot(2,2,1)
viewSlices(dataT,omega,m);
title('dataT, template');

subplot(2,2,2)
viewSlices(dataR,omega,m);
title('dataR, reference');

subplot(2,2,3)
viewSlices(Topt,omega,m);
title('dataT(y).*detDy');

subplot(2,2,4)
viewSlices((Topt-dataR),omega,m);
title('final residual');


text:
function polys = readPolygonSet(filename)
%READPOLYGONSET Read a set of simple polygons stored in a file.
%
%   POLY = readPolygonSet(FILENAME);
%   Returns the polygon stored in the file FILENAME.
%   Polygons are assumed to be stored in text files, without headers, with
%   x and y coordinates packed in two separate lines:
%     X11 X12 X13 ... X1N
%     Y11 Y12 Y13 ... Y1N
%     X21 X22 X23 ... X2N
%     Y21 Y22 Y23 ... Y2N
%
%   Each polygon may have a different number of vertices. The result is a
%   cell array of polygon, each cell containing a N-by-2 array representing
%   the vertex coordinates.
%
%   See also
%   polygons2d

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

% the set of polygons (no pre-allocation, as we do not know how many
% polygons are stored)
polys = {};

% index of polygon
p = 0;

% open file for reading
fid = fopen(filename, 'rt');

% use an infinite loop, terminated in case of EOF
while true
    % set of X, and Y coordinates
    line1 = fgetl(fid);
    line2 = fgetl(fid);

    % break loop if end of file is reached
    if line1 == -1
        break;
    end

    % create a new polygon by concatenating vertex coordinates
    p = p + 1;
    polys{p} = [str2num(line1)' str2num(line2)']; %#ok<AGROW,ST2NM>
end

% close file
fclose(fid);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function okay = vgg_mexcompile(funcName, varargin)
%VGG_MEXCOMPILE  Mex compile helper function
%
%   okay = vgg_mexcompile(funcName)
%   okay = vgg_mexcompile(..., sourceList)
%   okay = vgg_mexcompile(..., lastCompiled)
%
% Compile mexed function, given an optional list of source files. Can
% optionally check if source files have been modified since the last
% compilation, and only compile if they have.
%
%IN:
%   funcName - string containg the name of the function to compile
%   sourceList - cell array of source files to be compiled. Default:
%                {[funcName '.cxx']}.
%   lastCompiled - datenum of the current mex file. Default: 0 (i.e. force
%                  compilation).
%
%OUT:
%   okay - 1: function compiled; 0: up-to-date, no need to compile; -1:
%          compilation failed.

% $Id: vgg_mexcompile.m,v 1.3 2009/09/13 20:34:58 ojw Exp $

% Set defaults for optional inputs
sourceList = [funcName '.cxx'];
lastCompiled = 0;
% Parse inputs
for a = 1:numel(varargin)
   if iscell(varargin{a}) && ischar(varargin{a}{1})
       sourceList = varargin{a};
   elseif isnumeric(varargin{a}) && isscalar(varargin{a})
       lastCompiled = varargin{a};
   end
end

% Go to the directory containing the file
currDir = cd;
sourceDir = fileparts(which(funcName));
cd(sourceDir);

if lastCompiled
    compile = false;
    okay = 0;
    % Compile if current mex file is older than any of the source files
    for a = 1:numel(sourceList)
        dirSource = dir(sourceList{a});
        if ~isempty(dirSource) && datenum(dirSource.date) > lastCompiled
            compile = true;
            break;
        end
    end
else
    compile = true;
end

% Compile if we need to
if compile
    % Set the compiler flags
    flags = ['-O -I"' sourceDir '"'];
    switch mexext
        case 'mexsol'
            flags = [flags ' CC=gcc CFLAGS=-fPIC'];
        case {'mexglx', 'mexa64'}
            str = '"-O3 -ffast-math -funroll-loops"';
            flags = sprintf('%s CXXOPTIMFLAGS=%s LDCXXOPTIMFLAGS=%s LDOPTIMFLAGS=%s', flags, str, str, str);
        otherwise
    end

    % Call mex to compile the code
    cmd = sprintf('mex %s%s', flags, sprintf(' "%s"', sourceList{:}));
    disp(cmd);
    try
        eval(cmd);
        okay = 1;
    catch
        okay = -1;
        fprintf('ERROR while compiling %s\n', funcName);
    end
end

% Return to the original directory
cd(currDir);
return

repo: SelfExSR


text:
function pdf = gamma_pdf ( x, a, b, c )

%*****************************************************************************80
%
%% GAMMA_PDF evaluates the Gamma PDF.
%
%  Discussion:
%
%    PDF(X)(A,B,C) = EXP ( - ( X - A ) / B ) * ( ( X - A ) / B )^(C-1)
%      / ( B * GAMMA ( C ) )
%
%    GAMMA_PDF(A,B,C), where C is an integer, is the Erlang PDF.
%    GAMMA_PDF(A,B,1) is the Exponential PDF.
%    GAMMA_PDF(0,2,C/2) is the Chi Squared PDF with C degrees of freedom.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the argument of the PDF.
%    A <= X.
%
%    Input, real A, B, C, the parameters of the PDF.
%    A controls the location of the peak;  A is often chosen to be 0.0.
%    B is the "scale" parameter; 0.0 < B, and is often 1.0.
%    C is the "shape" parameter; 0.0 < C, and is often 1.0.
%
%    Output, real PDF, the value of the PDF.
%
  if ( x <= a )

    pdf = 0.0;

  else

    y = ( x - a ) / b;

    pdf = y^( c - 1.0 ) / ( b * gamma ( c ) * exp ( y ) );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [p1,p2,P1_l,P2_l] = ahmLin2ahmPnts(l)

% AHMLIN2AHMPNTS  AHM line to two AHM points conversion.
%   [P1,P2] = AHMLIN2AHMPNTS(L) extracts the two endpoints of the AHM line
%   L in the form of two HMG points with the same anchor of that of the
%   line.
%
%   [p1,p2,P1_l,P2_l] = AHMLIN2AHMPNTS(...) returns the Jacobians wrt L.

%   Copyright 2009 Teresa Vidal.

p1 = l(1:7,:);
p2 = l([1:3 8:11],:);

if nargout > 2

    P1_l = [eye(7) zeros(7,4)];
    P2_l = [eye(3) zeros(3,8) ; zeros(4,7) eye(4)];

end

return

%% jac

syms x y z a1 b1 c1 n1 a2 b2 c2 n2 real
l = [x y z a1 b1 c1 n1 a2 b2 c2 n2]';

[p1,p2,P1_l,P2_l] = ahmLin2ahmPnts(l);

simplify(P1_l - jacobian(p1,l))
simplify(P2_l - jacobian(p2,l))

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

%   # START GPL LICENSE

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

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

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function [ nodes_per_layer, n ] = hex_grid_approximate_n ( box, n_goal )

%*****************************************************************************80
%
%% HEX_GRID_APPROXIMATE_N seeks a hex grid of about N nodes.
%
%  Discussion:
%
%    The parameter NODES_PER_LAYER controls the number of nodes, but
%    in a somewhat obscure way.  This routine experiments with various
%    values until it is convinced it has the value of NODES_PER_LAYER
%    that comes as close as possible to producing N nodes.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 March 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real BOX(2,2), the lower and upper corners
%    of the rectangular region.
%
%    Input, integer N_GOAL, the desired number of nodes.
%
%    Output, integer NODES_PER_LAYER, the number of nodes per layer
%    which produces a mesh with about N_GOAL nodes.
%
%    Output, integer N, the number of nodes in the mesh.
%
  m = 2;

  if ( n_goal <= 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'HEX_GRID_APPROXIMATE_N - Fatal error!\n' );
    fprintf ( 1, '  Illegal input value of N_GOAL = %d\n', n_goal );
    error ( 'HEX_GRID_APPROXIMATE_N - Fatal error!' );
  end

  nodes_per_layer_low = 0;
  n_low = 0;

  nodes_per_layer = round ( sqrt ( n_goal ) );

  nodes_per_layer_high = n_goal;
  n_high = n_goal * n_goal;

  while ( 1 )

    n = hex_grid_n ( nodes_per_layer, box );

    if ( n == n_goal )
      break
    end

    if ( n < n_goal )
      nodes_per_layer_low = nodes_per_layer;
      n_low = n;
    else
      nodes_per_layer_high = nodes_per_layer;
      n_high = n;
    end

    if ( nodes_per_layer_low + 1 == nodes_per_layer_high )
      if ( n - n_low <= n_high - n )
        nodes_per_layer = nodes_per_layer_high;
        n = n_high;
      else
        nodes_per_layer = nodes_per_layer_low;
        n = n_low;
      end
      break
    end

    nodes_per_layer = ...
      round ( ( nodes_per_layer_low + nodes_per_layer_high ) / 2 );

  end

  return
end


text:
function [validBeacon, mpdu, validPacket] = nonHTBeaconRxMPDUDecode(payload)
% nonHTBeaconRxMPDUDecode Checks if the payload input is from a beacon
% packet and returns its MPDU information. This is a helper function for
% the non-HT beacon receiver example.
%
%   [VALIDBEACON, MPDU] = nonHTBeaconRxMPDUDecode(PAYLOAD) checks the
%   PAYLOAD binary input and sets VALIDBEACON to be true if PAYLOAD is from
%   a beacon packet. The MPDU information is also returned, but it is
%   meaningful only when VALIDBEACON is true.

% Copyright 2015-2016 The MathWorks, Inc.

%#codegen

coder.extrinsic('warning')

persistent crcDetector defaultMPDU

% Frame check sequence (FCS)
if isempty(crcDetector)
    crcDetector = comm.CRCDetector( ...
        [32 26 23 22 16 12 11 10 8 7 5 4 2 1 0], ...
        'InitialConditions', 1, ...
        'DirectMethod',      true, ...
        'FinalXOR',          1);
    defaultMPDU.Header = parseMPDUHeader(zeros(64,1), 1);
    defaultMPDU.FrameBody = unknownMPDU();
end

% Predefine output
validBeacon = false;
mpdu = defaultMPDU;
validPacket = false;

% If the CRC for the header decodes correctly (i.e. CRC flag is 0) and the
% received MPDU payload length is greater than or equal to the expected
% length, then turn on the MPDU decoder
if (length(payload) >= 32+16)
    % Compute CRC based on IEEE Std 802.11(TM)-2007 section 7.1.3.7.
    [frame, mpduCRCd] = crcDetector(double(payload));
    mpduCRC = mpduCRCd > 0; % Convert to logical

    % Extract header information
    [mpdu.Header, mpduHeaderLen] = parseMPDUHeader(frame, mpduCRC);

    validPacket = mpduCRC == 0;

    % Evaluate frame body
    if (mpduCRC == 1) || ((mpdu.Header.FrameCtrl.Type ~= 0) ...
            || (mpdu.Header.FrameCtrl.Subtype ~= 8))
        mpdu.FrameBody = unknownMPDU();
    else
        validBeacon = true;
        mpdu.FrameBody = parseBeaconMPDU(frame(mpduHeaderLen+1:end, 1), mpduCRC);
    end
end

end

function [mpduHeader, len] = parseMPDUHeader(frame, mpduCRC)
% Based on IEEE Std 802.11-2007 chapter 7.

mpduHeader.FrameCtrl       = parseFrameControl(frame(1:16,1), mpduCRC);
mpduHeader.DurationID      = zeros(1,2,'uint8');
mpduHeader.Address1        = zeros(1,6,'uint8');
mpduHeader.Address2        = zeros(1,6,'uint8');
mpduHeader.Address3        = zeros(1,6,'uint8');
mpduHeader.SequenceControl = zeros(1,2,'uint8');
mpduHeader.Address4        = zeros(1,6,'uint8');
len = (2+6+6+6+2+6)*8 + 16;

if (mpduCRC == 0)
    p = 1;
    mpduHeader.FrameCtrl  = parseFrameControl(frame(p:p+15,1), mpduCRC);
    p = p+16;
    mpduHeader.DurationID = bit2Octet(frame(p:p+15,1),2);
    p = p+16;
    mpduHeader.Address1   = bit2Octet(frame(p:p+47,1),6);
    p = p+48;

    isManagementFrame = (mpduHeader.FrameCtrl.Type == 0);
    isControlFrame = (mpduHeader.FrameCtrl.Type == 1);
    isDataFrame = (mpduHeader.FrameCtrl.Type == 2);
    isCTSFrame = isControlFrame && (mpduHeader.FrameCtrl.Subtype == 12);
    isACKFrame = isControlFrame && (mpduHeader.FrameCtrl.Subtype == 13);

    if ~(isCTSFrame || isACKFrame)
        mpduHeader.Address2 = bit2Octet(frame(p:p+47,1),6);
        p = p+48;
    end
    if isManagementFrame
        mpduHeader.Address3 = bit2Octet(frame(p:p+47,1),6);
        p = p+48;
    end
    if ~isControlFrame
        mpduHeader.SequenceControl = bit2Octet(frame(p:p+15,1),2);
        p = p+16;
    end
    if isDataFrame
        mpduHeader.Address4 = bit2Octet(frame(p:p+47,1),6);
        p = p+48;
    end
    len = p-1;
end
end

function frameCtrl = parseFrameControl(fcBits, mpduCRC)
if (mpduCRC == 0)
    frameCtrl.ProtocolVersion = uint8(bi2de(fcBits(1:2,1)'));
    frameCtrl.Type            = uint8(bi2de(fcBits(3:4,1)'));
    frameCtrl.Subtype         = uint8(bi2de(fcBits(5:8,1)'));
    frameCtrl.ToDS            = uint8(bi2de(fcBits(9,1)'));
    frameCtrl.FromDS          = uint8(bi2de(fcBits(10,1)'));
    frameCtrl.MoreFragments   = uint8(bi2de(fcBits(11,1)'));
    frameCtrl.Retry           = uint8(bi2de(fcBits(12,1)'));
    frameCtrl.PowerManagement = uint8(bi2de(fcBits(13,1)'));
    frameCtrl.MoreData        = uint8(bi2de(fcBits(14,1)'));
    frameCtrl.ProtectedFrame  = uint8(bi2de(fcBits(15,1)'));
    frameCtrl.Order           = uint8(bi2de(fcBits(16,1)'));
else
    frameCtrl.ProtocolVersion = uint8(0);
    frameCtrl.Type            = uint8(0);
    frameCtrl.Subtype         = uint8(0);
    frameCtrl.ToDS            = uint8(0);
    frameCtrl.FromDS          = uint8(0);
    frameCtrl.MoreFragments   = uint8(0);
    frameCtrl.Retry           = uint8(0);
    frameCtrl.PowerManagement = uint8(0);
    frameCtrl.MoreData        = uint8(0);
    frameCtrl.ProtectedFrame  = uint8(0);
    frameCtrl.Order           = uint8(0);
end
end

function frameBody = parseBeaconMPDU(frame, mpduCRC)
if (mpduCRC == 0)
    bitCount = uint32(1);
    frameBody.TimeStamp      = bit2Octet(frame(bitCount:bitCount+63,1),8);
    bitCount = bitCount + 64;
    frameBody.BeaconInterval = bit2Octet(frame(bitCount:bitCount+15,1),2);
    bitCount = bitCount + 16;
    frameBody.Capability     = parseCapabilities(frame(81:96), mpduCRC);
    bitCount = bitCount + 16;
else
    frameBody.TimeStamp      = zeros(1,8,'uint8');
    frameBody.BeaconInterval = zeros(1,2,'uint8');
    frameBody.Capability     = parseCapabilities(frame(81:96), mpduCRC);
    bitCount = uint32(97);
end

frameLength = length(frame);

frameBody.NumInfoElements = uint8(0);
macElementStruct = struct('ID', uint8(255), ...
    'Length', uint8(0), ...
    'Value', zeros(1,257,'uint8'));
frameBody.InfoElements = repmat(macElementStruct, 1, 60);

if (mpduCRC == 0)
    while bitCount < frameLength
        frameBody.NumInfoElements = frameBody.NumInfoElements + 1;
        frameBody.InfoElements(frameBody.NumInfoElements).ID = ...
            bit2Octet(frame(bitCount:bitCount+7,1), 1);
        bitCount = bitCount + 8;
        elemLength = bit2Octet(frame(bitCount:bitCount+7,1), 1);
        frameBody.InfoElements(frameBody.NumInfoElements).Length = elemLength;
        bitCount = bitCount + 8;
        for p=1:elemLength
            frameBody.InfoElements(frameBody.NumInfoElements).Value(p) = ...
                bit2Octet(frame(bitCount:bitCount+7),1);
            bitCount = bitCount+8;
        end
    end
end
end

function capability = parseCapabilities(frame, mpduCRC)
if (mpduCRC == 0)
    capability.ESS                = frame(1,1);
    capability.IBSS               = frame(2,1);
    capability.CFPollable         = frame(3,1);
    capability.CFPollRequest      = frame(4,1);
    capability.Privacy            = frame(5,1);
    capability.ShortPreamble      = frame(6,1);
    capability.PBCC               = frame(7,1);
    capability.ChannelAgility     = frame(8,1);
    capability.SpectrumManagement = frame(9,1);
    capability.QoS                = frame(10,1);
    capability.ShortSlotTime      = frame(11,1);
    capability.APSD               = frame(12,1);
    capability.Reserved           = frame(13,1);
    capability.DSSOFDM            = frame(14,1);
    capability.DelayedBlockAck    = frame(15,1);
    capability.ImmediateBlockAck  = frame(16,1);
else
    capability.ESS                = double(false);
    capability.IBSS               = double(false);
    capability.CFPollable         = double(false);
    capability.CFPollRequest      = double(false);
    capability.Privacy            = double(false);
    capability.ShortPreamble      = double(false);
    capability.PBCC               = double(false);
    capability.ChannelAgility     = double(false);
    capability.SpectrumManagement = double(false);
    capability.QoS                = double(false);
    capability.ShortSlotTime      = double(false);
    capability.APSD               = double(false);
    capability.Reserved           = double(false);
    capability.DSSOFDM            = double(false);
    capability.DelayedBlockAck    = double(false);
    capability.ImmediateBlockAck  = double(false);
end
end

function frameBody = unknownMPDU()
frameBody.TimeStamp      = zeros(1,8,'uint8');
frameBody.BeaconInterval = zeros(1,2,'uint8');
frameBody.Capability     = parseCapabilities(zeros(16,1), 1);

frameBody.NumInfoElements = uint8(0);
macElementStruct = struct('ID', uint8(255), ...
    'Length', uint8(0), ...
    'Value', zeros(1,257,'uint8'));
frameBody.InfoElements = repmat(macElementStruct, 1, 60);
end

function octets = bit2Octet(bits,numOctets)
octets = zeros(1, numOctets, 'uint8');
for p=1:numOctets
    octets(1,p) = (bits((p-1)*8+1:p*8)')*2.^(0:7)';
end
end

author: analogdevicesinc


text:
function [ point, seed ] = p00_sample ( test, m, n, seed )

%*****************************************************************************80
%
%% P00_SAMPLE samples points from the region in a problem.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer TEST, the index of the test problem
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real POINT(M,N), the coordinates of the points.
%
%    Output, integer SEED, a seed for the random number generator.
%
  if ( test == 1 )
    [ point, seed ] = p01_sample ( m, n, seed );
  elseif ( test == 2 )
    [ point, seed ] = p02_sample ( m, n, seed );
  elseif ( test == 3 )
    [ point, seed ] = p03_sample ( m, n, seed );
  elseif ( test == 4 )
    [ point, seed ] = p04_sample ( m, n, seed );
  elseif ( test == 5 )
    [ point, seed ] = p05_sample ( m, n, seed );
  elseif ( test == 6 )
    [ point, seed ] = p06_sample ( m, n, seed );
  elseif ( test == 7 )
    [ point, seed ] = p07_sample ( m, n, seed );
  elseif ( test == 8 )
    [ point, seed ] = p08_sample ( m, n, seed );
  elseif ( test == 9 )
    [ point, seed ] = p09_sample ( m, n, seed );
  elseif ( test == 10 )
    [ point, seed ] = p10_sample ( m, n, seed );
  elseif ( test == 11 )
    [ point, seed ] = p11_sample ( m, n, seed );
  elseif ( test == 12 )
    [ point, seed ] = p12_sample ( m, n, seed );
  elseif ( test == 13 )
    [ point, seed ] = p13_sample ( m, n, seed );
  elseif ( test == 14 )
    [ point, seed ] = p14_sample ( m, n, seed );
  elseif ( test == 15 )
    [ point, seed ] = p15_sample ( m, n, seed );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P00_SAMPLE - Fatal error!\n' );
    fprintf ( 1, '  Input value of TEST is out of range.\n' );
    error ( 'P00_SAMPLE - Fatal error!' );
  end

  return
end


text:
function A = warmUpExercise()
%WARMUPEXERCISE Example function in octave
%   A = WARMUPEXERCISE() is an example function that returns the 5x5 identity matrix

A = [];
% ============= YOUR CODE HERE ==============
% Instructions: Return the 5x5 identity matrix
%               In octave, we return values by defining which variables
%               represent the return values (at the top of the file)
%               and then set them accordingly.

% ===========================================
A = eye(5);

end

author: zhouxc

repo: Stanford-Machine-Learning-Course


text:
function s = Simpson2(a, b)
%% Simpson Rule for a prestablished f(x)
% a: is the lower limit of integration
% b: is the upper limit of integration
% for this excersice f(x)=x^2

h=(b-a)/4;
x(1)=a;
x(2)=a+(b-a)*1/4;
x(3)=a+(b-a)*1/2;
x(4)=a+(b-a)*3/4;
x(5)=b;

f(1)=myfunction(x(1));
f(2)=myfunction(x(2));
f(3)=myfunction(x(3));
f(4)=myfunction(x(4));
f(5)=myfunction(x(5));

s=h/3*sum([1 4 2 4 1].*f);


text:
function [qp]=da2(input)
% [qp] = da2(input)
% Dahlin's controller for processes of 2nd order.
% This function computes parameters of the controller (q0, q1, q2, p1, p2).
% Transfer function of the controller is as follows:
%
%            q0 + q1*z^-1 + q2*z^-2     q0 + q1*z^-1 + q2*z^-2
% G(z^-1) = ------------------------ = ------------------------
%                  1 - z^-1              1 + p1*z^-1 + p2*z^-2
%
% where p1=-1 and p2=0.
%
% Transfer function of the controlled system is:
%
%                  b1*z^-1
% Gs(z^-1) = -----------------------
%             1 + a1*z^-1 + a2*z^-2
%
% Input: input ... input parameters
%                  input(1) ... a1
%                  input(2) ... b1
%                  input(3) ... a2
%                  input(4) ... sample time T0
%                  input(5) ... adjustment factor B
% Output: qp ... controller parameters
%                qp(1) ... q0
%                qp(2) ... q1
%                qp(3) ... q2
%                qp(4) ... -1   (p1 of the controller)
%                qp(5) ... 0    (p2 of the controller)

a1 = input(1);
b1 = input(2);
a2 = input(3);
T0 = input(4);
B = input(5);

% check inputs
if (T0 <= 0)
   disp('da2.m: Input parameter T0 (sample time) must be greater than 0');
   disp('       parameter value has been changed to 1');
   T0 = 1;
end;
if (B <= 0)
   disp('da2.m: Input parameter B (adjustment factor) must be greater than 0');
   disp('       parameter value has been changed to 1e-6');
   B = 1e-6;
end;

Q = 1-exp(-T0/B);
Kp = -(a1+2*a2)*Q/b1;
Td = T0*a2*Q/(Kp*b1);
Ti = -T0/((1/(a1+2*a2))+1+Td/T0);

q0 = Kp*(1+T0/Ti+Td/T0);
q1 = -Kp*(1+2*Td/T0);
q2 = Kp*Td/T0;
p1 = -1;
p2 = 0;

qp=[q0; q1; q2; p1; p2];

author: Sable

repo: mcbench-benchmarks


text:
function [phi, theta, psi] = rotMatXYZ2Eul(rotMat)
    % Function to convert rotation matrix to Euler angles (extrinsic x-y-z convention).
    % Used in yaw to convert between global and local reference frames.
    %
    % Parameters
    % ------------
    %     rotMat : float array [3,3]
    %         Rotational matrix to transform vectors to the new frame
    %
    % Returns
    % ------------
    %     phi : float [rad]
    %         Angle of rotation about the x-axis
    %     theta : float [rad]
    %         Angle of rotation about the y-axis
    %     psi : float [rad]
    %         Angle of rotation about the z-axis

    phi = atan2(-rotMat(2,3), rotMat(3,3));
    theta = asin(rotMat(1,3));
    psi = atan2(-rotMat(1,2), rotMat(1,1));
end

author: WEC-Sim

repo: WEC-Sim

save_path: github-repos/MATLAB/WEC-Sim-WEC-Sim


text:
% SP_EVALUATE_ELEMENT_LIST: compute the basis functions in a given list of elements.
%
%     sp = sp_evaluate_element_list (space, msh_elems, 'option1', value1, ...)
%
% INPUTS:
%
%    space:     object defining the space of discrete functions (see sp_vector)
%    msh_elems: msh structure containing the information of quadrature or
%               visualization points, for a given list of elements
%               (see msh_cartesian/msh_evaluate_element_list)
%   'option', value: additional optional parameters, currently available options are:
%
%              Name     |   Default value |  Meaning
%           ------------+-----------------+----------------------------------
%            value      |      true       |  compute shape_functions
%            gradient   |      false      |  compute shape_function_gradients
%            divergence |      false      |  compute shape_function_divs
%            curl       |      false      |  compute shape_function_curls
%
% OUTPUT:
%
%    sp: struct representing the discrete function space, with the following fields:
%              (see the article for a detailed description)
%
%    FIELD_NAME      (SIZE)                                     DESCRIPTION
%    ncomp           (scalar)                                   number of components of the functions of the space
%    ndof            (scalar)                                   total number of degrees of freedom
%    ndof_dir        (ncomp x ndim matrix)                      for each component, number of degrees of freedom along each direction
%    nsh_max         (scalar)                                   maximum number of shape functions per element
%    nsh             (1 x msh_col.nel vector)                   actual number of shape functions per each element
%    connectivity    (nsh_max x msh_col.nel vector)             indices of basis functions that do not vanish in each element
%    shape_functions (ncomp x msh_col.nqn x nsh_max x msh_col.nel)  basis functions evaluated at each quadrature node in each element
%    shape_function_gradients
%       (ncomp x rdim x msh_col.nqn x nsh_max x msh_col.nel)    basis function gradients evaluated at each quadrature node in each element
%    shape_function_divs (msh_col.nqn x nsh_max x msh_col.nel)  basis function divergence evaluated at each quadrature node in each element
%    shape_function_curls
%         2D:  (msh_col.nqn x nsh_max x msh_col.nel)            basis function curl evaluated at each quadrature node in each element
%         3D:  (3 x msh_col.nqn x nsh_max x msh_col.nel)
%
% Copyright (C) 2009, 2010, 2011 Carlo de Falco
% Copyright (C) 2011, 2015, 2019 Rafael Vazquez
%
%    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/>.

function sp = sp_evaluate_element_list (space, msh, varargin)

value = true;
gradient = false;
divergence = false;
curl = false;
hessian = false;
if (~isempty (varargin))
  if (~rem (length (varargin), 2) == 0)
    error ('sp_evaluate_element_list: options must be passed in the [option, value] format');
  end
  for ii=1:2:length(varargin)-1
    if (strcmpi (varargin {ii}, 'value'))
      value = varargin {ii+1};
    elseif (strcmpi (varargin {ii}, 'gradient'))
      gradient = varargin {ii+1};
    elseif (strcmpi (varargin {ii}, 'curl'))
      curl = varargin {ii+1};
    elseif (strcmpi (varargin {ii}, 'divergence'))
      divergence = varargin {ii+1};
    elseif (strcmpi (varargin {ii}, 'hessian'))
      hessian = varargin {ii+1};
    else
      warning ('Ignoring unknown option %s', varargin {ii});
    end
  end
end

grad_param = gradient || divergence || curl || hessian;
value_param = value || grad_param;
div_param = false; curl_param = false;
switch (lower (space.transform))
  case {'curl-preserving'}
    curl_param = curl;
  case {'div-preserving'}
    div_param = divergence;
end

sp = sp_evaluate_element_list_param (space, msh, 'value', value_param, 'gradient', grad_param, 'divergence', div_param, 'curl', curl_param, 'hessian', hessian);

if (isempty (msh.elem_list))
  return
end

switch (lower (space.transform))
  case {'grad-preserving'}
    sp = sp_vector_grad_preserving_transform (sp, msh, value, gradient, curl, divergence, hessian);
  case {'curl-preserving'}
    sp = sp_vector_curl_preserving_transform (sp, msh, value, curl);
    if (gradient || divergence || hessian)
      warning ('Gradient, divergence and hessian not implemented for curl-preserving transformation')
    end
  case {'div-preserving'}
    sp = sp_vector_div_preserving_transform (sp, msh, value, gradient, curl, divergence);
    if (hessian)
      warning ('Hessian not implemented for div-preserving transformation')
    end
end

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
%% Analyzing Investment Strategies with CVaR Portfolio Optimization in MATLAB
%
% Robert Taylor
% The MathWorks, Inc.

% Copyright (C) 2012 The MathWorks, Inc.

%% Introduction

% This script is a "superscript" that organizes the scripts for the webinar "Analyzing Investment
% Strategies with CVaR Portfolio Optimization in MATLAB." It describes what each script does and
% shows the order in which the scripts should be examined.

%% Instructions

% The file structure for these scripts has a top-level folder that contains these scripts and should
% have two folders 'data' and 'source' with data and analytics. To add these folders to the path,
% start in the folder with the scripts and execute the commands

setlocalpaths

% The script cvarwebinar_scenarios.m, which simulates scenarios, must be run before subsequent
% scripts can be run because it generates a file BuyWriteScenarios.mat that is needed for these
% scripts. The script to generate scenarios can take about one hour on a typical computer and
% requires that the computer be a 64-bit machine. It creates a 12MB file in the ./data folder. Make
% sure that the script to generate scenarios is run in the folder that contains this script so that
% the scenarios file ends up in the correct data folder.

%% Theory

% This script illustrates basic features of covered-call strategies and provides a theoretical
% analysis of issues regarding slippage due to assignment and re-investment.

cvarwebinar_theory

%% From Theory to Reality

% This script illustrates an event-driven simulation of total returns for uncovered and covered
% positions based on a single realization of an underlying stock. It moves from the simplicity of
% theory to the messiness of reality as it models and simulates various contributions to slippage.

cvarwebinar_reality

%% Calibration

% This script is the first of the sequence of scripts to illustrate a complete workflow to analyze a
% covered-call strategy for a universe of 26 stocks. Given total return price data, this script
% illustrates maximum likelihood calibration of the assumed geometric Brownian motion process for
% the universe of stocks.

cvarwebinar_calibration

%% Scenario Generation

% This script is the second of the sequence of scripts that generates scenarios for uncovered and
% covered positions to be used for subsequent portfolio optimization and analysis. This is the
% slowest script since it generates scenarios by simulation of investment actions during the course
% of an investment period.

cvarwebinar_scenarios

%% Normality Tests

% This script is the third of the sequence of scripts that examines the statistical properties of
% the scenarios. As the probability of early exercise increases during an investment period, the
% distribution of the log of covered-call returns becomes increasingly non-normal.

cvarwebinar_normality

%% Optimization

% This script is the final of the sequence of scripts that performs several portfolio optimization
% steps with both CVaR and mean-variance portfolio optimization. The difference in results between
% the two types of optimization is examined to provide greater insights into the normative
% implications of covered-call strategies.

cvarwebinar_optimization

%% References
%
% # P. Bernstein (1998), _Against the Gods: The Remarkable Story of Risk_, Wiley.
% # F. Black (1975), "Fact and Fantasy in the Use of Options," _Financial Analysts Journal_, Vol. 31,
% No. 4, pp. 36-41 and 61-72.
% # P. Glassermann (1991), _Monte Carlo Methods in Financial Engineering_, Springer.
% # I. Karatzas and S. Shreve (1991), _Brownian Motion and Stochastic Calculus_, 2nd ed., Springer.
% # H. Markowitz (1952), "Portfolio Selection," _Journal of Finance_, Vol. 7, No. 1, pp. 77-91.
% # R. Merton, M. Scholes, and M. Gladstein (1978), "The Returns and Risk of Alternative Call Option
% Portfolio Investment Strategies," _Journal of Business_, Vol. 51, No. 2, pp. 183-242.
% # R. T. Rockafellar and S. Uryasev (2002). "Conditional Value-at-Risk for General Loss
% Distributions," _Journal of Banking and Finance_, Vol. 26, pp. 1443-1471.


text:
% Riemannian Dictionary Learning and Sparse Coding
%  if you plan to use this code, please cite
% Cherian, Anoop, and Suvrit Sra. "Riemannian Dictionary Learning and Sparse Coding for Positive Definite Matrices." arXiv preprint arXiv:1507.02772 (2015).
% This code is released on the BSD3 license.
% This code should not be used for non-commercial purposes.
% The authors are not liable to any loss or damage caused by running this code.
% For any issues/bugs, please contact anoop.cherian@gmail.com.

% Copyright (c) 2016, Anoop Cherian
% All rights reserved.
% Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
% 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
% 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
% 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
% PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
% INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
% OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
% OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

addpath ./spalgos/;
addpath(genpath('./manopt/'));
addpath ./tools;
key = 1234;
randn('state', key); rand('twister', key);

n=100; num_atoms=20; active_size=5; d=5;
[data, B, A] = generate_data(n, d, num_atoms, active_size); data=data';
data = data(randperm(n));

% sparsity for varying lambda:
lambda = 0.1;
[BB, alpha, obj] = Fast_Riem_DL(data, num_atoms, active_size, 'sim', lambda, 1234);

author: hiroyuki-kasai

repo: ClassifierToolbox


text:
function H = elec_meshplot(X,Y,Z,interp,grid_res)

% elec_meshplot - Plot a mesh and points for 3D electrode positions.
%
% Usage: H = elec_meshplot(X,Y,Z [,interp] [,grid_resolution])
%
% Uses an interpolated mesh at 0.25 cm resolution, unless
% user specifies alternative grid resolution (optional).
%
% Interpolation can be 'linear' or 'cubic'.  Unless specified,
% cubic is the default.
%
% Returns a handle to the figure created.
%

% $Revision: 1.1 $ $Date: 2009-04-28 22:13:55 $

% Licence: Gnu GPL
% Author: Darren.Weber_at_radiology.ucsf.edu
% Created:  18/05/00
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  H = figure('NumberTitle','off','Name','Electrode Mesh Plot','Position',[200 200 650 500],'PaperUnits','centimeters','PaperType','A4','Units','centimeters');

  colormap(gray);

  if ~exist( 'grid_res', 'var' )  grid_res = 0.25;  end

  xi =   (min(X) - 2):grid_res:(max(X) + 2);
  yi =  [(min(Y) - 2):grid_res:(max(Y) + 2)]';

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Use 'linear' for a linear grid or 'cubic' for a polynomial grid

  if ~exist( 'interp', 'var' )  interp = 'cubic';  end

  [Xi,Yi,Zi] = griddata(X,Y,Z,xi,yi,interp);

  surf(Xi,Yi,Zi);

  brighten(0.75);

  mesh(Xi,Yi,Zi);

  hold on;

  plot3(X,Y,Z,'ro');

  rotate3d;

  clear xi yi Xi Yi Zi X Y Z;

%  view(2);  % top view
%  view(   0,0);  % back view
%  view(-180,0);  % front view
%  view(-90, 0);  % left side view
%  view( 90, 0);  % right side view

author: PatternRecognition

repo: OpenBMI


text:
%compute dxheadsm_mm

function [data,units]=compute_dxheadsm_mm(trx,n)

larvae=trx.exp2flies{n};
numlarvae=numel(larvae);
dxheadsm_mm=cell(1,numlarvae);
for i=1:numlarvae
    larva=larvae(i);
    dxheadsm_mm{1,i}=(trx(larva).xheadsm_mm(2:end)-trx(larva).xheadsm_mm(1:end-1))./trx(larva).dt;
end

units=parseunits('mm/s');
data=dxheadsm_mm;

author: kristinbranson

repo: JAABA


text:
function [ y, m, d, f ] = frac_borrow_common ( y, m, d, f )

%*****************************************************************************80
%
%% FRAC_BORROW_COMMON borrows fractions from days in a Common YMDF date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 December 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, integer M, integer D, real F,
%    a YMDF date.
%
%    Output, integer Y, integer M, integer D, real F,
%    the YMDF date has been adjusted.
%
  while ( f < 0.0 )

    f = f + 1.0;

    d = d - 1;

  end

  [ y, m, d ] = day_borrow_common ( y, m, d );

  return
end