text:
function [M, Mx, My] = diffusionTermRadial2D(D)
% This function uses the central difference scheme to discretize a 2D
% diffusion term in the form \grad . (D \grad \phi) where u is a face vactor
% It also returns the x and y parts of the matrix of coefficient.
%
% SYNOPSIS:
%
%
% PARAMETERS:
%
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

% extract data from the mesh structure
Nr = D.domain.dims(1);
Ntheta = D.domain.dims(2);
G=reshape(1:(Nr+2)*(Ntheta+2), Nr+2, Ntheta+2);
DR = repmat(D.domain.cellsize.x, 1, Ntheta);
DTHETA = repmat(D.domain.cellsize.y', Nr, 1);
dr = 0.5*(DR(1:end-1,:)+DR(2:end,:));
dtheta = 0.5*(DTHETA(:,1:end-1)+DTHETA(:,2:end));
rp = repmat(D.domain.cellcenters.x, 1, Ntheta);
rf = repmat(D.domain.facecenters.x, 1, Ntheta);

% define the vectors to store the sparse matrix data
iix = zeros(3*(Nr+2)*(Ntheta+2),1);	iiy = zeros(3*(Nr+2)*(Ntheta+2),1);
jjx = zeros(3*(Nr+2)*(Ntheta+2),1);	jjy = zeros(3*(Nr+2)*(Ntheta+2),1);
sx = zeros(3*(Nr+2)*(Ntheta+2),1);	sy = zeros(3*(Nr+2)*(Ntheta+2),1);
mnx = Nr*Ntheta;	mny = Nr*Ntheta;

% reassign the east, west, north, and south diffusivity vectors for the
% code readability
De = rf(2:Nr+1,:).*D.xvalue(2:Nr+1,:)./(rp.*dr(2:Nr+1,:).*DR(2:Nr+1,:));
Dw = rf(1:Nr,:).*D.xvalue(1:Nr,:)./(rp.*dr(1:Nr,:).*DR(2:Nr+1,:));
Dn = D.yvalue(:,2:Ntheta+1)./(rp.*rp.*dtheta(:,2:Ntheta+1).*DTHETA(:,2:Ntheta+1));
Ds = D.yvalue(:,1:Ntheta)./(rp.*rp.*dtheta(:,1:Ntheta).*DTHETA(:,2:Ntheta+1));

% calculate the coefficients for the internal cells
AE = reshape(De,mnx,1);
AW = reshape(Dw,mnx,1);
AN = reshape(Dn,mny,1);
AS = reshape(Ds,mny,1);
APx = reshape(-(De+Dw),mnx,1);
APy = reshape(-(Dn+Ds),mny,1);

% build the sparse matrix based on the numbering system
rowx_index = reshape(G(2:Nr+1,2:Ntheta+1),mnx,1); % main diagonal x
iix(1:3*mnx) = repmat(rowx_index,3,1);
rowy_index = reshape(G(2:Nr+1,2:Ntheta+1),mny,1); % main diagonal y
iiy(1:3*mny) = repmat(rowy_index,3,1);
jjx(1:3*mnx) = [reshape(G(1:Nr,2:Ntheta+1),mnx,1); reshape(G(2:Nr+1,2:Ntheta+1),mnx,1); reshape(G(3:Nr+2,2:Ntheta+1),mnx,1)];
jjy(1:3*mny) = [reshape(G(2:Nr+1,1:Ntheta),mny,1); reshape(G(2:Nr+1,2:Ntheta+1),mny,1); reshape(G(2:Nr+1,3:Ntheta+2),mny,1)];
sx(1:3*mnx) = [AW; APx; AE];
sy(1:3*mny) = [AS; APy; AN];

% build the sparse matrix
kx = 3*mnx;
ky = 3*mny;
Mx = sparse(iix(1:kx), jjx(1:kx), sx(1:kx), (Nr+2)*(Ntheta+2), (Nr+2)*(Ntheta+2));
My = sparse(iiy(1:ky), jjy(1:ky), sy(1:ky), (Nr+2)*(Ntheta+2), (Nr+2)*(Ntheta+2));
M = Mx + My;

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


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

%STOMP PARAMETERS
%number of particles
K = 50;

%repeat experiment number of times
parameters.n_repeat = 30;
parameters.experiment_name = 'experiment2B_K50_N30.mat';
parameters.animate = 0;

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

'ended'
parameters.experiment_name

repo: ARTE


text:
% Calculate the dot product of coefficient vectors of two polynomials
%
%  Syntax:  (pdot is the shortened alias for PolynomialDotProduct)
%    >> d = PolynomialDotProduct(f,g)
%    >> d = pdot(f,g)
%
%   Input:   f --- (string)            polynomial
%            g --- (string)            polynomial
%
%  Output:   d --- (numeric)  dot product of coeffcient vectors of f and g
%
%  Example:  >> d = pdot('3+x*y-y^3*x^2+x^5','5-4*x^2*y^3+9*x*y')
%            d =
%                  28

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function site = db_get_site_info(stalist,datetime,dbname)

% SITE = DB_GET_SITE_INFO(STATION_LIST,DATE,DATABASE_NAME) This
% This function takes a list of stations and returns a structure containing
% relavent fields from the site table. The database must include a
% site table. The DATE field is required to resolve conflicts when a
% station configuration has changed through time. Typically the date will
% match waveform, origin or arrival data that is being extracted from the
% database. The date should be given is a string format recognized by
% antelope, or as a matlab-format numeric time.
%
%
% Example:
%    site = db_get_site_info({'AUL' 'AUH' 'AUW'},'7/1/2007','dbtmp')
%
%

% Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
% $Date$
% $Revision$

% FORMAT INPUT
if ischar(stalist)
   stalist = cellstr(stalist);
end
site.sta = reshape( stalist , numel(stalist), 1 );

% PREP DAY OF YEAR
if isreal(datetime)
   datetime = datestr(datetime,'mm/dd/yy HH:MM:SS.FFF');
end
jday = num2str(yearday(str2epoch(datetime)));

% GET SITE INFO
db = dbopen(dbname,'r');
db = dblookup_table(db,'site');
nlist = [];
for n = 1:length(site.sta);
	db1 = dbsubset(db,['sta==''' site.sta{n} '''']);
    db1 = dbsubset(db1,['ondate<=''' jday ''' ']);
    db1 = dbsubset(db1,['offdate>=''' jday ''' || offdate==NULL']);
    recnum = dbquery(db1,'dbRECORD_COUNT');
	[ site.ondate(n) , site.offdate(n) , site.lat(n) , site.lon(n) , site.elev(n) ] = dbgetv(db1,'ondate','offdate','lat','lon','elev');
end
dbclose(db);
site.ondate = reshape( site.ondate , numel(site.ondate), 1 );
site.offdate = reshape( site.offdate , numel(site.offdate), 1 );
site.lat = reshape( site.lat , numel(site.lat), 1 );
site.lon = reshape( site.lon , numel(site.lon), 1 );
site.elev = reshape( site.elev , numel(site.elev), 1 );

% plot(correlation(w));
%figure
%plot(w(1),'b')
%hold on;
%plot(w1(1),'g')
%plot(w2(1),'r')

author: geoscience-community-codes

repo: GISMO


text:
function [ h ] = my_drawArrow(x1,y1,x2,y2,headwidth,props)

h = annotation('arrow');
set(h,'parent', gca, ...
    'position', [x1,y1,x2-x1,y2-y1], ...
    'HeadLength', 10, 'HeadWidth', headwidth, 'HeadStyle', 'cback1', ...
    props{:} );
end

author: angeloyeo

repo: gongdols

save_path: github-repos/MATLAB/angeloyeo-gongdols


text:
function test_bug1412

% MEM 3gb
% WALLTIME 00:10:00
% DEPENDENCY ft_read_header loadcnt

hdr = ft_read_header(dccnpath('/home/common/matlab/fieldtrip/data/test/bug1412/Sub1_1.cnt'));

assert(hdr.nSamples == 1505200);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
classdef CourseraClient < handle

% Copyright (C) 2013, Georgia Tech Research Corporation
% see the LICENSE file included with this software

    properties
        root_url

        part_layout
        layout
        parent

        ui_colors
        ui_size

        login_field
        password_field
        submit_button

        part_list
        course_url
        course_week
    end

    methods

        function obj = CourseraClient(root_url)
            obj.root_url = root_url;
            obj.ui_colors = struct('gray',  [220 220 220]/255, ...
                                   'green', [ 57 200  67]/255, ...
                                   'red',   [221  23  31]/255, ...
                                   'dgray', [183 183 183]/255, ...
                                   'lgray', [242 242 242]/255, ...
                                   'black', [ 48  48  48]/255, ...
                                   'white', [255 255 255]/255);
            obj.part_list = simiam.containers.ArrayList(10);
        end

        function load_ui(obj)
            obj.populate_assignment();
            obj.create_ui();
        end

        function populate_assignment(obj)

            % Read in XML file
            blueprint = xmlread(fullfile(obj.root_url, 'coursera.xml'));

            % Parse XML file for robot configurations
            e = blueprint.getElementsByTagName('course').item(0);
            obj.course_url = ['http://class.coursera.org/' char(e.getAttribute('identifier'))];
            week_t = str2double(e.getAttribute('week'));

            assignment = blueprint.getElementsByTagName('assignment').item(week_t-1);
            obj.course_week = str2double(assignment.getAttribute('week'));

            part_list_t = assignment.getElementsByTagName('part');
            for k = 0:(part_list_t.getLength-1)
                part_i = part_list_t.item(k);

                identifier = char(part_i.getAttribute('identifier'));
                title = char(part_i.getAttribute('title'));
                test = str2func(['simiam.test.' char(part_i.getAttribute('test')) '.run_test']);

                part_i = struct('identifier', identifier, ...
                                'title', title, ...
                                'test', test, ...
                                'checkbox', 0);

                obj.part_list.appendElement(part_i);
            end
        end

        function create_ui(obj)
            obj.parent = figure('MenuBar', 'none', ...
                                'NumberTitle', 'off', ...
                                'Name', 'Sim.I.am: Programming Assignment Submission', ...
                                'Color', obj.ui_colors.gray);

            n_parts = obj.part_list.length();
            obj.ui_size = [0 0 3*128 256+n_parts*(48+10)];

            ui_new_size = get(obj.parent, 'Position');
            ui_new_size(3:4) = obj.ui_size(3:4);
            screen_size = get(0, 'ScreenSize');
            margins = (screen_size(3:4)-obj.ui_size(3:4))/2;
            ui_new_size(1:2) = margins;
            set(obj.parent, 'Position', ui_new_size);

            obj.layout = GridLayout(obj.parent, 'RowHeight', {96,32,32,'*',32}, ....
                                                'ColWidth', {128, '*'}, ...,
                                                'LMargin', 10, 'RMargin', 6, ...
                                                'Gap', 5, ...
                                                'CellColor', obj.ui_colors.gray);

            set(obj.layout.Container, 'BackgroundColor', obj.ui_colors.gray);
            MergeCells(obj.layout, 4, [1 2]);
            MergeCells(obj.layout, 5, [1 2]);
            Update(obj.layout);

            button_string = 'Login:';
            ui_args = {'Style','pushbutton', 'String', button_string, 'ForegroundColor', 'k', 'FontWeight', 'bold', 'BackgroundColor', obj.ui_colors.gray, 'Callback', ''};
            ui_parent = obj.layout.Cell(2,1);
            aButtonLabel = uicontrol(ui_parent, ui_args{:});
            set(findjobj(aButtonLabel), 'BorderPainted', 0);
            set(aButtonLabel, 'Enable', 'inactive');

            ui_args = {'Style','edit', 'String', '', 'ForegroundColor', 'k', 'FontWeight', 'normal', 'BackgroundColor', obj.ui_colors.white, 'Callback', ''};
            ui_parent = obj.layout.Cell(2,2);
            obj.login_field = uicontrol(ui_parent, ui_args{:});

            button_string = 'Password:';
            ui_args = {'Style','pushbutton', 'String', button_string, 'ForegroundColor', 'k', 'FontWeight', 'bold', 'BackgroundColor', obj.ui_colors.gray, 'Callback', ''};
            ui_parent = obj.layout.Cell(3,1);
            aButtonLabel = uicontrol(ui_parent, ui_args{:});
            set(aButtonLabel, 'Enable', 'inactive');
            set(findjobj(aButtonLabel), 'BorderPainted', 0);

            ui_args = {'Style','edit', 'String', '', 'ForegroundColor', 'k', 'FontWeight', 'normal', 'BackgroundColor', obj.ui_colors.white, 'Callback', ''};
            ui_parent = obj.layout.Cell(3,2);
            obj.password_field = uicontrol(ui_parent, ui_args{:});

            icon_file = fullfile(obj.root_url, 'resources/simiam-round-small.png');
            if(isunix)
                icon_url = ['file://' icon_file];
            else
                icon_url = strrep(['file:/' icon_file],'\','/');
            end
            button_string = ['

'];
            ui_args = {'Style','pushbutton', 'String', button_string, 'ForegroundColor', 'k', 'FontWeight', 'bold', 'BackgroundColor', obj.ui_colors.gray, 'Callback', ''};
            ui_parent = obj.layout.Cell(1,1);
            aButtonLabel = uicontrol(ui_parent, ui_args{:});
            set(findjobj(aButtonLabel), 'BorderPainted', 0);
            set(aButtonLabel, 'Enable', 'inactive');

            button_string = 'Remember to use your submission login and password!';
            ui_args = {'Style','pushbutton', 'String', button_string, 'ForegroundColor', 'k', 'FontWeight', 'bold', 'BackgroundColor', obj.ui_colors.gray, 'Callback', ''};
            ui_parent = obj.layout.Cell(1,2);
            aButtonLabel = uicontrol(ui_parent, ui_args{:});
            set(findjobj(aButtonLabel), 'BorderPainted', 0);
            set(aButtonLabel, 'Enable', 'inactive');

            n_parts = obj.part_list.length();
            part_panel = uipanel( ...
                                'Parent', obj.layout.Cell(4,1), ...
                                'Title', sprintf('Assignment parts:'), ...
                                'Units', 'pixels', ...
                                'BackgroundColor', obj.ui_colors.gray);
            obj.part_layout = GridLayout(part_panel, 'NumRows', n_parts, 'RowHeight', 48,...
                                                     'ColWidth', {32, '*', 32}, ...,
                                                     'Gap', 5, ...
                                                     'CellColor', obj.ui_colors.gray);

            for i = 1:n_parts
                part = obj.part_list.elementAt(i);

                ui_parent = obj.part_layout.Cell(i,1);
                ui_args = {'Style','checkbox', 'String', '', 'BackgroundColor', obj.ui_colors.gray, 'Callback', ''};
                part_i = uicontrol(ui_parent, ui_args{:});
                set(part_i, 'Value', 1);
                part.checkbox = part_i;

%                 button_string = ['
' list_items{i} '
'];
                button_string = part.title;
                ui_args = {'Style','edit', 'String', button_string, 'ForegroundColor', 'k', 'FontWeight', 'bold', 'BackgroundColor', obj.ui_colors.gray, 'HorizontalAlignment', 'left'};
                ui_parent = obj.part_layout.Cell(i,2);
                aButtonLabel = uicontrol(ui_parent, ui_args{:});
%                 set(findjobj(aButtonLabel), 'BorderPainted', 0);
                set(aButtonLabel, 'Enable', 'inactive');
                set(aButtonLabel, 'Max', 10);
%                 set(aButtonLabel, 'HorizontalAlignment', 'left');
                part.label = aButtonLabel;

                ui_args = {'Style','pushbutton', 'String', button_string, 'ForegroundColor', 'k', 'FontWeight', 'bold', 'BackgroundColor', obj.ui_colors.gray, 'HorizontalAlignment', 'left'};
                ui_parent = obj.part_layout.Cell(i,3);
                aButtonLabel = uicontrol(ui_parent, ui_args{:});
                set(findjobj(aButtonLabel), 'BorderPainted', 0);
                set(aButtonLabel, 'Enable', 'inactive');
                obj.ui_set_button_icon(aButtonLabel, 'ui_status_unknown.png');
                part.status = aButtonLabel;

                obj.part_list.replaceElementAt(part, i);
            end

            button_string = 'Submit to Coursera for Grading';
            ui_args = {'Style','pushbutton', 'String', button_string, 'ForegroundColor', 'k', 'FontWeight', 'bold', 'BackgroundColor', obj.ui_colors.dgray, 'Callback', @obj.ui_pressed_submit};
            ui_parent = obj.layout.Cell(5,1);
            obj.submit_button = uicontrol(ui_parent, ui_args{:});

            Update(obj.layout);
            Update(obj.part_layout);

            set(obj.parent, 'CloseRequestFcn', @obj.ui_close);

            if exist('coursera_login_data.mat','file')
                load('coursera_login_data.mat');
                set(obj.login_field, 'String', login)
                set(obj.password_field, 'String', password);
            end
        end

        function ui_set_button_icon(obj, ui_button, icon)
            icon_file = fullfile(obj.root_url, 'resources/icons', icon);
            if(isunix)
                icon_url = ['file://' icon_file];
            else
                icon_url = strrep(['file:/' icon_file],'\','/');
            end
            button_string = [''];
            set(ui_button, 'String', button_string);
        end

        function ui_pressed_submit(obj, src, event)

            % save login and password information
            login = get(obj.login_field, 'String');
            password = get(obj.password_field, 'String');
            save('coursera_login_data.mat', 'login', 'password');

            for i = 1:length(obj.part_list)
                part_i = obj.part_list.elementAt(i);
                obj.ui_set_button_icon(part_i.status, 'ui_status_unknown.png');
                set(part_i.label, 'String', part_i.title);
            end
            drawnow;

            for i = 1:length(obj.part_list)
                part_i = obj.part_list.elementAt(i);
                if get(part_i.checkbox, 'Value')

                    [login, ch, signature, auxstring] = obj.get_challenge_from_coursera(login, part_i.identifier);
                    if isempty(login) || isempty(ch) || isempty(signature)
                        % Some error occured, error string in first return element.
                        fprintf('\n!! Error: %s\n\n', login);
                        return
                    end

                    % Attempt Submission with Challenge
                    ch_resp = obj.generate_challenge_response(login, password, ch);

                    [result, str] = obj.submit_solution_to_coursera(login, ch_resp, part_i.identifier, ...
                           obj.evaluate_test_fcn(part_i.test, auxstring), obj.read_fcn_source(part_i.test), signature);

                    if (result == 0)
                        obj.ui_set_button_icon(part_i.status, 'ui_status_ok.png');
                    else
                        obj.ui_set_button_icon(part_i.status, 'ui_status_error.png');
%                         fprintf('Feedback: %s\n', str);
                        label = get(part_i.label, 'String');
                        set(part_i.label, 'String', sprintf('%s',[label ': <Feedback> ' str]));
                    end
                end
            end

        end

        function ui_close(obj, src, event)
            delete(src);
        end

        %% Coursera-specific functionality (derived from assignment-api-examples.zip)
%
        function [email,ch,signature,auxstring] = get_challenge_from_coursera(obj, login, part_id)
            challenge_url = [obj.course_url '/assignment/challenge'];
            str = urlread(challenge_url(), 'post', {'email_address', login, 'assignment_part_sid', part_id, 'response_encoding', 'delim'});

            str = strtrim(str);
            r = struct;
            while(numel(str) > 0)
                [f, str] = strtok (str, '|');
                [v, str] = strtok (str, '|');
                r = setfield(r, f, v);
            end

            email = getfield(r, 'email_address');
            ch = getfield(r, 'challenge_key');
            signature = getfield(r, 'state');
            auxstring = getfield(r, 'challenge_aux_data');
        end

        function [str] = generate_challenge_response(obj,email, passwd, challenge)
            str = simiam.util.SecureHashAlgorithm.generate_sha1_hash([challenge passwd]);
        end

        function out = evaluate_test_fcn(obj, fcn, aux_string)
          out = sprintf('%s', fcn(aux_string, obj.root_url));
        end

        function src = read_fcn_source(obj, filename)
            src = '';
        end

        function [result, str] = submit_solution_to_coursera(obj, email, ch_resp, sid, output, ...
                source, signature)

            params = {'assignment_part_sid', sid, ...
                'email_address', email, ...
                'submission', simiam.util.EncodingScheme.encode_in_base64(output, ''), ...
                'submission_aux', simiam.util.EncodingScheme.encode_in_base64(source, ''), ...
                'challenge_response', ch_resp, ...
                'state', signature};

            submit_url = [obj.course_url '/assignment/submit'];
            [str, status] = urlread(submit_url, 'post', params);

            % Parse str to read for success / failure
            result = -1;

            if strcmp(str, 'Fantastic!')
                result = 0;
            else

            end

        end

    end

end

author: jdelacroix

repo: simiam

save_path: github-repos/MATLAB/jdelacroix-simiam


text:
%%
%   MATLAB code for DeepFool
%
%   adversarial_DeepFool_caffe(x,net):
%   computes the adversarial perturbations for a Caffe's model
%
%   INPUTS
%   x: image in W*H*C format
%   net: Caffe's network (without loss layer - do not forget to enable 'force_backward')
%   opts: A struct contains parameters (see README)
%   OUTPUTS
%   r_hat: minimum perturbation
%   l_hat: adversarial label
%   l: classified label
%   itr: number of iterations
%
%   please cite: S. Moosavi-Dezfooli, A. Fawzi, P. Frossard: DeepFool: a simple and accurate method to fool deep neural networks.
%                In Computer Vision and Pattern Recognition (CVPR 2016), IEEE, 2016.
%%
function [r_hat,l_hat,l,itr] = adversarial_DeepFool_caffe(x,net,opts)
size_x = size(x);
x = reshape(x,numel(x),1);
l = f(x,1);

if(nargin==3)
    adv = adversarial_perturbation(x,l,@Df,@f,opts);
else
    adv = adversarial_perturbation(x,l,@Df,@f);
end

l_hat = adv.new_label;
r_hat = reshape(adv.r,size_x);
itr = adv.itr;

    function out = f(y,flag)
        y = reshape(y,size_x);

        out = net.forward({y}); %do forward pass
        out = out{1}'; %convert 'out' from a cell array to a matrix

        %flag==0:compute the outputs
        %flag==1:compute the label
        if flag==1
            [~,out] = max(out);
        end
    end

    function dzdx = Df(y,label,idx)
        y = reshape(y,size_x);
        net.forward({y}); %do forward pass

        for i=1:numel(idx)
            dzdy = zeros(net.blobs(net.blob_names{end}).shape,'single');

            dzdy(idx(i)) = 1;

            res = net.backward({dzdy}); %do backward pass
            dzdx(:,i) = reshape(res{1},numel(y),1);
        end
        dzdx = dzdx-repmat(dzdx(:,idx==label),1,numel(idx));
    end
end

repo: DeepFool


text:
function im_scale = prep_im_for_blob_size(im_size, target_size, max_size)

    im_size_min = min(im_size(1:2));
    im_size_max = max(im_size(1:2));
    im_scale = double(target_size) / im_size_min;

    % Prevent the biggest axis from being more than MAX_SIZE
    if round(im_scale * im_size_max) > max_size
        im_scale = double(max_size) / double(im_size_max);
    end
end

author: ShaoqingRen


text:
function displayMobject(mObject,fieldname,note)
%   M-object information display
%   displayMobject(mObject,fieldname,note);

%   Designed and coded by Hideki Kawahara
%   27/Feb./2005
%   Copyright(c) 2005, Hideki Kawahara
%   05/Oct./2005 minor bug fix

fs = mObject.samplingFrequency;
tFrame = mObject.frameUpdateInterval;
switch fieldname
    case 'spectrogram'
        figure
        [nrow,ncolumn]=size(mObject.spectrogram);
        timeSpan = [0 (ncolumn-1)*tFrame];
        dBsgram = 20*log10(mObject.spectrogram);
        maxSgramdB = max(max(dBsgram));
        imagesc(timeSpan, [0 fs/2],max(dBsgram,maxSgramdB-70));
        axis('xy');
        set(gca,'fontsize',14);
        xlabel('time (ms)');
        ylabel('frequency (Hz)');
        title([note ' time span 0 ' num2str(timeSpan(2),10) ' (ms) ' datestr(now)]);
    case 'waveform'
        figure
        x = mObject.waveform;
        timeSpan = (0:length(x)-1)/fs*1000;
        plot(timeSpan,x);grid on;
        axis([timeSpan(1) timeSpan(end) 1.1*[min(x) max(x)]]);
        set(gca,'fontsize',14);
        xlabel('time (ms)');
        title([note ' time span 0 ' num2str(round(timeSpan(end)),8) ' (ms) ' datestr(now)]);
    case {'anchorFrequency', 'anchorTimeLocation'}
        figure
        [nrow,ncolumn]=size(mObject.spectrogram);
        timeSpan = [0 (ncolumn-1)*tFrame];
        dBsgram = 20*log10(mObject.spectrogram);
        maxSgramdB = max(max(dBsgram));
        imagesc(timeSpan, [0 fs/2],max(dBsgram,maxSgramdB-70));
        axis('xy');
        set(gca,'fontsize',14);
        xlabel('time (ms)');
        ylabel('frequency (Hz)');
        title([note ' time span 0 ' num2str(timeSpan(2),10) ' (ms) ' datestr(now)]);
        if length(mObject.anchorTimeLocation)>0
            hold on;
            for ii=1:length(mObject.anchorTimeLocation)
                hh = plot(mObject.anchorTimeLocation(ii)*[1 1],[0 fs/2],'w:');
                set(hh,'linewidth',2);
                if sum(mObject.anchorFrequency(ii,:)>0)>0
                    nFrequency = sum(mObject.anchorFrequency(ii,:)>0);
                    anchorFrequencyVector = mObject.anchorFrequency(ii,mObject.anchorFrequency(ii,:)>0); % 05/Oct./2005 HK
                    for jj=1:nFrequency
                        hh=plot(mObject.anchorTimeLocation(ii),anchorFrequencyVector(jj),'ok');
                        set(hh,'markersize',9);
                        set(hh,'linewidth',2);
                        hh=plot(mObject.anchorTimeLocation(ii),anchorFrequencyVector(jj),'.w');
                        set(hh,'markersize',7);
                        set(hh,'linewidth',4);
                    end;
                end;
            end;
            hold off;
        end;
end;

author: HidekiKawahara


text:
function str = sim_dampedOscillator(f0,tau,fs,varnum)
% return a string corresponding to the equation of motion of a damped
% oscillator. This is meant to be compatible with format used in
% sim_genVARModelFromEq()
%
% f0:   fundamental frequency (Hz)
% tau:  damping time (e.g. 100). Larger values --> more sinusoidal
% fs:   samping rate
% varnum: the index of the channel/variable being model
% (e.g. x2 --> varnum = 2)
%
% (C) Tim Mullen, May, 2011. SCCN/INC UCSD

str = sprintf('{2*exp(-1/%f)*cos(2*pi*%f/%f)}*x%d(t-1) + -exp(-2/%f)*x%d(t-2)',tau,f0,fs,varnum,tau,varnum);

repo: REST


text:
function y=binocdf(x,n,p);

% BINOCDF binomial cumulative distribution function
%
% Y=BINOCDF(X,N,P) returns the binomial cumulative distribution
% function with parameters N and P at the values in X.
%
% See also BINOPDF and STATS (Matlab statistics toolbox)

% compute the cumulative probability for all values up to the maximum
c = cumsum(binopdf(0:max(x(:)),n,p));
y = c(x+1);

% fix rounding errors
y(y<0) = 0;
y(y>1) = 1;


text:
% create a sub network that takes in waveforms and produces fourier
% coefficients

function layer = genNetworkSTFT2LogSpec(stftLayer, useLog, logConst)
layer = stftLayer;

layer{end+1} = PowerNode(stftLayer{end}.dim(1));
if useLog
    if nargin<3
        logConst = 0.00;
    end
    layer{end+1} = LogarithmNode(layer{end}.dim(1), logConst);
end

layer = ConnectLinearGraph(layer);
layer = FinishLayer_obj(layer);

end

author: singaxiong

repo: SignalGraph


text:
classdef World < handle

% Copyright (C) 2013, Georgia Tech Research Corporation
% see the LICENSE file included with this software

    properties
        robots
        obstacles
        parent
        apps
        root_path
    end

    methods
        function obj = World(parent)
            obj.parent = parent;
            obj.robots = simiam.containers.ArrayList(10);
            obj.obstacles = mcodekit.list.dl_list();
            obj.apps = simiam.containers.ArrayList(10);
            obj.root_path = '';
        end

        function build_from_file(obj, root, file, origin)

            % Read in XML file
            blueprint = xmlread(file);
            obj.root_path = root;

            % Parse XML file for robot configurations
            app_list = blueprint.getElementsByTagName('app').item(0);
            app = char(app_list.getAttribute('type'));

            r = str2func(strcat('simiam.app.', app));
            obj.apps.appendElement(r(root));

            if(strcmp(origin, 'launcher') || strcmp(origin, 'testing') || strcmp(origin, 'hardware'))

                robot_list = blueprint.getElementsByTagName('robot');

                for k = 0:(robot_list.getLength-1)
                   robot = robot_list.item(k);

                   type = char(robot.getAttribute('type'));

                   s = robot.getElementsByTagName('supervisor').item(0);
                   spv = char(s.getAttribute('type'));

                   pose = robot.getElementsByTagName('pose').item(0);
                   x = str2double(pose.getAttribute('x'));
                   y = str2double(pose.getAttribute('y'));
                   theta = str2double(pose.getAttribute('theta'));

                   r = obj.add_robot(type, spv, x, y, theta);

                   driver = robot.getElementsByTagName('driver').item(0);
                   if(~isempty(driver) && strcmp(origin, 'hardware'))
                       hostname = char(driver.getAttribute('ip'));
                       port = str2double(driver.getAttribute('port'));
                       r.add_hardware_link(hostname,port);
                       r.open_hardware_link();
                   end
                end
            end

            % Parse XML file for obstacle configurations
            obstacle_list = blueprint.getElementsByTagName('obstacle');

            for i = 0:(obstacle_list.getLength-1)
               obstacle = obstacle_list.item(i);

               pose = obstacle.getElementsByTagName('pose').item(0);
               x = str2double(pose.getAttribute('x'));
               y = str2double(pose.getAttribute('y'));
               theta = str2double(pose.getAttribute('theta'));

               geo = obstacle.getElementsByTagName('geometry').item(0);
               point_list = geo.getElementsByTagName('point');

               obstacle_geometry = zeros(point_list.getLength, 2);
               for j=0:(point_list.getLength-1)
                  point = point_list.item(j);
                  obstacle_geometry(j+1,1) = str2double(point.getAttribute('x'));
                  obstacle_geometry(j+1,2) = str2double(point.getAttribute('y'));
               end

               obj.add_obstacle(x, y, theta, obstacle_geometry);
            end
        end

        function aRobot = add_robot(obj, type, spv, x, y, theta)
            pose = simiam.ui.Pose2D(x, y, theta);

            r = str2func(strcat('simiam.robot.', type));
            robot = r(obj.parent, pose);

            r = str2func(strcat('simiam.controller.', spv));
            supervisor = r();

            supervisor.attach_robot(robot, pose);

            %            parameter_file = fullfile(obj.root_path, 'parameters.xml');
            %            [file, path, extension] = uigetfile(parameter_file, 'Pick a XML file with the parameters for the PID regulator.');
            %
            %            parameter_file = fullfile(path, file);
            %            supervisor.configure_from_file(parameter_file);

            s = struct('robot', robot, 'supervisor', supervisor, 'pose', pose);
            %            obj.robots.append_key(s);

            obj.robots.appendElement(s);
            anApp = obj.apps.elementAt(1);
            anApp.supervisors.appendElement(supervisor);

            aRobot = s.robot;
        end

        function add_obstacle(obj, x, y, theta, geometry)
           pose = simiam.ui.Pose2D(x, y, theta);
           obj.obstacles.append_key(struct('obstacle', simiam.simulator.Obstacle(obj.parent, pose, geometry), 'pose', pose));
        end

    end

end


text:
% This example shows how to calculate the TM_01 mode of a step-index
% fiber.  The magnetic field for the TM_01 mode is azimuthally
% symmetric, which means that both Hx and Hy are of the same
% magnitude, making this an excellent test case for the
% full-vector modesolver.  It also compares the numerically
% calculated result with the exact solution.  Note that the
% modes look virtually identical, but there is a discrepancy
% betweeen the numerically-computed eigenvalue and the exact
% eigenvalue.  This is thought to be caused by the the
% inaccuracy of representing a curved core-clad interface with
% a staircase function.

% Refractive indices:
nco = 2.5;          % core index
ncl = 1.5;          % cladding index
r = 0.30;           % core radius (um)

side = 0.2;         % space on side (um)

dx = 0.002;         % grid size (horizontal)
dy = 0.002;         % grid size (vertical)

lambda = 1;         % wavelength
nmodes = 1;         % number of modes to compute

% Boundary conditions for antisymmetric mode
boundary = '0A0S';

% Set up finite difference mesh:
[x,y,xc,yc,nx,ny,eps] = fiber([nco,ncl],[r],side,dx,dy);

% Now we stretch out the mesh at the boundaries:
[x,y,xc,yc,dx,dy] = stretchmesh(x,y,[96,0,96,0],[4,1,4,1]);

% Solve for mode (using transverse H modesolver)
[Hx,Hy,neff] = wgmodes (lambda, nco, nmodes, dx, dy, eps, ...
                         boundary);
fprintf(1,'neff (finite difference) = %8.6f\n',neff);

% Now solve for the exact eigenmodes
V = 2*pi*r/lambda*sqrt(nco^2-ncl^2);
U = fzero(@(U) ...
           nco^2*besselj(1,U)/(U*besselj(0,U)) + ...
           ncl^2*besselk(1,sqrt(V^2-U^2))/ ...
           (sqrt(V^2-U^2)*besselk(0,sqrt(V^2-U^2))), 3.2);
W = sqrt(V^2-U^2);
neff0 = sqrt(nco^2 - (U/(2*pi*r/lambda))^2);
fprintf(1,'neff (exact solution)    = %8.6f\n',neff0);

rho = sqrt(x.^2*ones(size(y)) + ones(size(x))*y.^2);
sinphi = ones(size(x))*y./rho;
cosphi = x*ones(size(y))./rho;
Hx0 = zeros(size(rho));
Hy0 = zeros(size(rho));
kv = find(rho < r);
Hx0(kv) = -sinphi(kv).*besselj(1,U*rho(kv)./r)/besselj(1,U);
Hy0(kv) = cosphi(kv).*besselj(1,U*rho(kv)./r)/besselj(1,U);
kv = find(rho >= r);
Hx0(kv) = -sinphi(kv).*besselk(1,W*rho(kv)./r)/besselk(1,W);
Hy0(kv) = cosphi(kv).*besselk(1,W*rho(kv)./r)/besselk(1,W);
[hmax,kv] = max(abs(Hx0(:)));
hmax = Hx0(kv);
Hx0 = Hx0/hmax;
Hy0 = Hy0/hmax;

figure(1);

subplot(221);
contourmode(x,y,Hx);
title('Hx (finite difference)');

subplot(222);
contourmode(x,y,Hy);
title('Hy (finite difference)');

subplot(223);
contourmode(x,y,Hx0);
title('Hx (exact)');

subplot(224);
contourmode(x,y,Hy0);
title('Hy (exact)');

author: Sable

repo: mcbench-benchmarks


text:
function dUpsilonS = lfmaaGradientSigmaUpsilonMatrix(gamma, sigma2, ...
    t1, t2, mode)

% LFMAAGRADIENTSIGMAUPSILONMATRIX Gradient of upsilon matrix aa wrt sigma
% FORMAT
% DESC computes the gradient wrt sigma of a portion of the LFMAA kernel.
% ARG gamma : Gamma value for system.
% ARG sigma2 : length scale of latent process.
% ARG t1 : first time input (number of time points x 1).
% ARG t2 : second time input (number of time points x 1).
% ARG mode : operation mode, according to the derivative (mode 0,
% derivative wrt t1, mode 1 derivative wrt t2)
% RETURN upsilon : result of this subcomponent of the kernel for the given values.
%
% COPYRIGHT : Mauricio Alvarez, 2010
%
% SEEALSO : lfmaaComputeUpsilonMatrix.m

% KERN

sigma = sqrt(sigma2);
gridt1 = repmat(t1, 1, length(t2));
gridt2 = repmat(t2', length(t1), 1);
timeGrid = gridt1 - gridt2;

dUpsilon = lfmapGradientSigmaUpsilonMatrix(gamma, sigma2, t1, t2, 1-mode);

if mode == 0
    dUpsilonS = (gamma^2)*dUpsilon - (2/(sqrt(pi)*sigma))*exp(-(timeGrid.^2)./sigma2).* ...
        ((2/sigma2- (2*timeGrid/sigma2).^2).*((gamma + 2*timeGrid/sigma2).* ...
        (1/sigma - 2*(timeGrid.^2)/sigma^3) + 4*timeGrid/sigma^3) ...
        - (gamma + 2*timeGrid/sigma2).*(-4/sigma^3 + 16*timeGrid.^2/sigma^5)) ...
        - (16/(sqrt(pi)*sigma^6))*timeGrid.*exp(-(timeGrid.^2)./sigma2).* ...
        (5 - 2*timeGrid.^2/sigma2);
else
    dUpsilonS = (gamma^2)*dUpsilon - (2/(sqrt(pi)*sigma))*exp(-(timeGrid.^2)./sigma2).* ...
        ((2/sigma2- (2*timeGrid/sigma2).^2).*((gamma + 2*timeGrid/sigma2).* ...
        (1/sigma - 2*(timeGrid.^2)/sigma^3) + 4*timeGrid/sigma^3) ...
        - (gamma + 2*timeGrid/sigma2).*(-4/sigma^3 + 16*timeGrid.^2/sigma^5)) ...
        - (16/(sqrt(pi)*sigma^6))*timeGrid.*exp(-(timeGrid.^2)./sigma2).* ...
        (5 - 2*timeGrid.^2/sigma2) - (2*gamma^2/(sqrt(pi)*sigma))*exp(-gamma*t1)* ...
        (((gamma-2*t2/sigma2).*(1/sigma - 2*t2.^2/sigma^3) - 4*t2/sigma^3).*exp(-t2.^2/sigma2)).';
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function obj=stat_summary(obj,varargin)
% stat_summary Display summarized data for each value of X
%
% Example syntax (default arguments): gramm_object.stat_summary('type','ci','geom','lines','setylim',false)
% For each unique value of x, this can display various estimates of
% the location and variability of the corresponding y distribution.
% The optional 'name',value pairs can ne the following:
% - 'type':
%       - 'ci' : display the mean and the 95% confidence
%       interval of the mean (based on the assumption of a normal
%       distribution)
%       - 'bootci' : display the mean and the 95% confidence
%       interval of the mean computed by bootstrap
%       - 'sem' : display the mean and the standard error of
%       the mean
%       - 'std' : display the mean and the standard deviation
%       and the mean
%       - 'quartile': display the 25% 50% (median) and 75%
%       percentiles
%       - '95percentile': display the median and 2.5 and 97.5
%       percentiles
%       - 'fitnormalci'
%       - 'fitpoissonci'
%       - 'fitbinomialci'
%       - 'fit95percentile'
%       - @function : provide a the handle to a custom function that takes
%       y values (as an n_repetitions x n_data_points matrix) and returns
%       both the central value and the CI with a matrix [y_central ; yc_CI_low ; y_CI_high]
%       (with y_central, and y_CIs being 1 x n_data_points arrays).
%       Example that uses the trimmed mean instead of regular mean:
%
%       custom_statfun = @(y)([trimmean(y,2.5);bootci(500,{@(ty)trimmean(ty,2.5),y},'alpha',0.05)]);
%       gramm_object.stat_summary('type', custom_statfun)
%
% - 'geom' (possibility to combine them using a cellstr, e.g. 'geom',{'bar','black_errorbar'} ):
%       - 'line': displays a line that connects the central locations
%       (mean,median)
%       - 'lines': displays a line that connects the central locations
%       and lighter lines for the variabilities
%       - 'area': displays a line that connects the central locations
%       and a transparent area for the variabilities. WARNING:
%       this changes the renderer to opengl and disables proper
%       vector output on older matlab versions
%       - 'area_only': displays the variabilities only, using a transparent area
%       - 'solid_area': displays a line that connects the locations
%       and a solid area for the variabilities. Use this for
%       export to vector output.
%       - 'errorbar': displays error bars for variabilities.
%       - 'black_errorbar': displays black error bars for variabilities.
%       - 'bar': displays the locations as bars
%       - 'edge_bar': displays the locations as bars with black edge
%       - 'point': displays the locations as points
% - 'setylim': set to true if you want the y axis limits to be
% set by the summarized data instead of the underlying data
% points.
% - 'interp': Use to interpolate the output, takes the same parameters as
% interp1 in order to specify the interpolation type. When the
% polar mode is specified as closed, or when 'interp','polar' is used
% the interpolation uses interpft, which supposes regular sampling around the circle.
% - 'interp_in': Use to (linearly) interpolate the input. This is intended
% for input given as cells when the x value is different for
% each cell and not aligned. The argument corresponds to the
% number of points used to generate the interpolation. Ideally
% for this the number of number of points should be higher than
% the x resolution of the data, otherwise some data will be
% unused.
% - 'bin_in': Use to bin the input. This is intended for input
% given as a 1-D array, creates bins over x and computes the summary
% over the binned data. Argument corresponds to the number of
% bins
% - 'dodge': use to dodge the plotted elements depending on
% color (recommended for 'bar', 'errorbar', 'black_errorbar').
% A value of 0 deactivates dodging. Other values set the space
% between the dodged elements as ratio of the x intervals.
% - 'width': use to set the width of bars and error bars (error
% bars are 1/4 th the width of bars).
%
% A setting of 'dodge',1,'width',1 will create bars that are
% fully dodged (ie don't overlap, but are not separated) and where
% all bars occupy the full interval between the x values, e.g.:
% * __    *
% *|  |__ *
% *|  |  |*
% *|__|__|*
%
% A setting of 'dodge',0.8,'width',0.8 will have fully dodged bars that are
% not separated, but only occupy 50% of the space between x
% values, e.g.:
% *  _    *
% * | |_  *
% * | | | *
% * |_|_| *
%
%A setting of 'dodge',0.8,'width',0.6 will add some
% spacing between the bars, e.g.:
% *  _      *
% * | |  _  *
% * | | | | *
% * |_| |_| *
%
% A setting of 'dodge',0.8,'width',1 will create dodged but overlapping
% bars, e.g.:
% *  ___    *
% * |  _|_  *
% * | |   | *
% * |_|___| *

p=inputParser;
my_addParameter(p,'type','ci'); %'95percentile'
my_addParameter(p,'geom','area');
my_addParameter(p,'dodge',[]);
my_addParameter(p,'width',[]);
my_addParameter(p,'setylim',false);
my_addParameter(p,'interp','none');
my_addParameter(p,'interp_in',-1);
my_addParameter(p,'bin_in',-1);
parse(p,varargin{:});

obj.geom=vertcat(obj.geom,{@(dobj,dd)my_summary(dobj,dd,p.Results)});
obj.results.stat_summary={};
end

function hndl=my_summary(obj,draw_data,params)

%Advanced defaults
if isempty(params.dodge)
    if sum(strcmp(params.geom,'bar'))>0 && draw_data.n_colors>1 %If we have a bar as geom, we dodge
        params.dodge=0.6;
    else
        params.dodge=0;
    end
end

if isempty(params.width) %If no width given
    if params.dodge>0 %Equal to dodge if dodge given
        params.width=params.dodge*0.8;
    else
        params.width=0.5;
    end
end

if iscell(draw_data.x) || iscell(draw_data.y) %If input was provided as cell/matrix

    if params.interp_in>0
        %If requested we interpolate the input
        uni_x=linspace(obj.var_lim.minx,obj.var_lim.maxx,params.interp_in);
        [x,y]=cellfun(@(x,y)deal(uni_x,interp1(x,y,uni_x,'linear')),draw_data.x,draw_data.y,'UniformOutput',false,'ErrorHandler',@(st,a,b)deal([],[]));
        y=padded_cell2mat(y);
        if isempty(y) %likely because we had only single points
            y = nan(size(uni_x));
            disp('Error in summary input interpolation... nothing plotted')
        end
    else
        %If not we just make a padded matrix for fast
        %computations (we'll assume that X are roughly at the
        %same location for the same indices)
        x=padded_cell2mat(draw_data.x);
        y=padded_cell2mat(draw_data.y);
        uni_x=nanmean(x);
        %Add a check for X alignment
        x_diff=max(x)-min(x);
        if  size(x,1)>1 && any(x_diff(1:end-1)>diff(uni_x)/10) %More than a tenth of delta x variation
            warning(['some repeated X values are misaligned (max ' num2str(max(x_diff)) '), use ''interp_in'' in stat_summary() or live with the consequences'])
        end
    end

    if params.bin_in>0
        warning('bin_in in stat_summary() not supported for Matrix/Cell X/Y inputs');
    end

    if ischar(params.type) && ~isempty(strfind(params.type,'fit'))
        %If we have a params.type using distributions fits we
        %can't vectorize the call to computeci so we do it in a for
        %loop
        ymean=zeros(length(uni_x),1);
        yci=zeros(length(uni_x),2);
        for ind_x=1:length(uni_x)
            [ymean(ind_x),yci(ind_x,:)]=computeci(y(:,ind_x),params.type,obj.stat_options.alpha,obj.stat_options.nboot);
        end
    else
        if size(y,1)==1
            ymean=y;
            yci=nan(length(uni_x),2);
        else
            [ymean,yci]=computeci(y,params.type,obj.stat_options.alpha,obj.stat_options.nboot);
        end
    end

else %If input was provided as 1D array

    x=comb(draw_data.x);
    y=comb(draw_data.y);

    if params.bin_in>0
        %If X binning was requested we do it
        binranges=linspace(obj.var_lim.minx,obj.var_lim.maxx,params.bin_in+1);
        bincenters=(binranges(1:(end-1))+binranges(2:end))/2;
        [~,binind]=my_histcounts(x,binranges,'count');
        uni_x=bincenters;
        sel=binind~=0; %histcounts can return zero as bin index if NaN data we remove them here
        x=bincenters(binind(sel));
        y=y(sel);
    else
        %         if sum(strcmp(params.geom,'area'))>0 || sum(strcmp(params.geom,'line'))>0 ||...
        %                 sum(strcmp(params.geom,'lines'))>0  || sum(strcmp(params.geom,'solid_area'))>0
        %             %To avoid interruptions in line and area plots we
        %             %compute uniques over current data only
        %             uni_x=unique(x); %Sorted is the default
        %         else
        %             if obj.x_factor
        %                 %If x is a factor we space everything as one
        %                 uni_x=obj.var_lim.minx:1:obj.var_lim.maxx;
        %             else
        %                 %compute unique Xs at the facet level to avoid
        %                 %weird bar sizing issues when dodging and when
        %                 %colors are missing
        %                 facet_x=comb(draw_data.facet_x);
        %                 uni_x=unique(facet_x);
        %             end
        %         end

        uni_x=unique(x);

        %Here we need to implement a loose 'unique' because of
        %potential numerical errors
        uni_x(diff(uni_x)<1e-10)=[];
    end

    if params.interp_in>0
        warning('interp_in in stat_summary() not supported for non Matrix/Cell X/Y inputs');
    end

    ymean=nan(length(uni_x),1);
    yci=nan(length(uni_x),2);

    %Loop over unique X values
    for ind_x=1:length(uni_x)
        %And here we have a loose selection also because of
        %potential numerical errors
        ysel=y(abs(x-uni_x(ind_x))<1e-10);

        if ~isempty(ysel)
            [ymean(ind_x),yci(ind_x,:)]=computeci(ysel,params.type,obj.stat_options.alpha,obj.stat_options.nboot);
        end
    end
end

%Do we set the y limits according to the smoothed curves or to
%the original data ?
if params.setylim
    if sum(sum(isnan(yci)))~=numel(yci) %We only do this if yci is not weird
        if obj.firstrun(obj.current_row,obj.current_column) %Initialize for the first run in the subplot
            obj.plot_lim.maxy(obj.current_row,obj.current_column)=max(max(yci));
            obj.plot_lim.miny(obj.current_row,obj.current_column)=min(min(yci));
        else %Update for subsequent runs in the subplot
            if max(max(yci))>obj.plot_lim.maxy(obj.current_row,obj.current_column)
                obj.plot_lim.maxy(obj.current_row,obj.current_column)=max(max(yci));
            end
            if min(min(yci))<obj.plot_lim.miny(obj.current_row,obj.current_column)
                obj.plot_lim.miny(obj.current_row,obj.current_column)=min(min(yci));
            end
        end
    end
end

%When we do bar plots we want to have zero in the y axis anyway
if sum(strcmp(params.geom,'bar'))>0
    if obj.plot_lim.miny(obj.current_row,obj.current_column)>0 %Values above zero -> change miny
        obj.plot_lim.miny(obj.current_row,obj.current_column)=0;
    end
    if obj.plot_lim.maxy(obj.current_row,obj.current_column)<0 %Values below zero -> change maxy
        obj.plot_lim.maxy(obj.current_row,obj.current_column)=0;
    end
end

%Do we interpolate the summary results for display ?
if ~strcmp(params.interp,'none')
    if size(yci,1)>2
        yci=yci';
    end

    if obj.polar.is_polar && obj.polar.is_polar_closed && ~strcmp(params.interp,'polar')
        disp([params.interp ' interpolation overriden, ''polar'' used']);
        params.interp='polar'; %%If the plot is polar we override the interpolation type do an optimal fft interpolation
    end

    if strcmp(params.interp,'polar')
        %Perform checks on uni_x
        dx=unique_no_nan(diff([uni_x ; uni_x(1)+2*pi])); %compute step
        if any(abs(diff(dx))>1e-12) %handle numerical precision problems (exact version would be length(dx)>1 )
            disp('ERROR: ''polar'' interpolation requires periodic sampling, displayed results are incorrect');
        end
        uni_x=uni_x(1):pi/50:uni_x(1)+99*pi/50;
        ymean=interpft(ymean,100);
        tmp_yci1=interpft(yci(1,:),100);
        tmp_yci2=interpft(yci(2,:),100);
        yci=[tmp_yci1 ; tmp_yci2];
    else
        %For non polar plots we do a regular interpolation
        new_x=linspace(min(uni_x),max(uni_x),100);
        ymean=interp1(uni_x,ymean,new_x,params.interp);
        tmp_yci1=interp1(uni_x,yci(1,:),new_x,params.interp);
        tmp_yci2=interp1(uni_x,yci(2,:),new_x,params.interp);
        yci=[tmp_yci1 ; tmp_yci2];
        uni_x=new_x;
    end
end

%If X were modified
if params.bin_in>0 || params.interp_in>0 || ~strcmp(params.interp,'none')
    %We reinitialize dodging parameters for plotci to work correctly
    draw_data.dodge_avl_w=uni_x(2)-uni_x(1);
    draw_data.dodge_fallback=true;
    draw_data.dodge_x=1;
    draw_data.dodge_n=draw_data.n_colors;
    draw_data.dodge_ind=draw_data.color_index;

    %draw_data.dodge_x=shiftdim(uni_x);
    %draw_data.dodge_n=repmat(draw_data.n_colors,length(uni_x),1);
    %draw_data.dodge_ind=repmat(draw_data.color_index,length(uni_x),1);
end

%Store results
obj.results.stat_summary{obj.result_ind,1}.x=uni_x;
obj.results.stat_summary{obj.result_ind,1}.y=ymean;
obj.results.stat_summary{obj.result_ind,1}.yci=yci;

%Do the actual plotting
hndl=plotci(obj,uni_x,ymean,yci,draw_data,params.geom,params.dodge,params.width);

%Copy handles
if isstruct(hndl)
    hnames=fieldnames(hndl);
    for k=1:length(hnames)
        obj.results.stat_summary{obj.result_ind,1}.(hnames{k})=hndl.(hnames{k});
    end
else
    disp('Nothing plotted... Error in summary computation?')
end
end

function [ymean,yci]=computeci(y,type,alpha,nboot)

ymean=nanmean(y);

%Check number of samples
nsamp=size(y,1);
if nsamp<3 && strcmp(type,'bootci')
    disp('Less than 3 samples for bootstrap CI computation...Skipping')
    yci=repmat([NaN NaN],size(y,2));
    return;
end
if nsamp<2
    disp('Less than 2 samples for CI computation...Skipping')
    yci=repmat([NaN NaN],size(y,2));
    return;
end

if isa(type,'function_handle')
    try
        temp=type(y);
        ymean=temp(1,:);
        yci=temp(2:3,:);
    catch ME
        disp(['Error in custom summary computation: ' ME.message ' ...skipping']);
        yci=repmat([NaN NaN],size(y,2));
        ymean=nan(size(y,2),1);
    end
    return;
end

try
    switch type
        case 'bootci'
            yci=bootci(nboot,{@(ty)nanmean(ty),y},'alpha',alpha);
        case 'ci'
            %ci=1.96*nanstd(y)./sqrt(sum(~isnan(y)));
            %Correction for small samples (equivalent to fitnormalci)
            ci=tinv(1-alpha/2,sum(~isnan(y))-1).*nanstd(y)./sqrt(sum(~isnan(y)));
            yci=bsxfun(@plus,ymean,[-ci;ci]);
        case 'std'
            ci=nanstd(y);
            yci=bsxfun(@plus,ymean,[-ci;ci]);
        case 'sem'
            ci=nanstd(y)./sqrt(sum(~isnan(y)));
            yci=bsxfun(@plus,ymean,[-ci;ci]);
        case 'quartile'
            ymean=nanmedian(y);
            yci=prctile(y,[25 75]);
        case '95percentile'
            ymean=nanmedian(y);
            yci=prctile(y,[2.5 97.5]);
        case 'fitnormalci'
            pd=fitdist(y,'Normal');
            ymean=pd.mean();
            ci=paramci(pd,alpha);
            yci=ci(:,1)';
        case 'fitpoissonci'
            pd=fitdist(y,'Poisson');
            ymean=pd.mean();
            ci=paramci(pd,alpha);
            yci=ci(:,1)';
        case 'fitnegbinomialci'
            pd=fitdist(y,'NegativeBinomial');
            ymean=pd.mean();
            ci=paramci(pd,alpha);
            yci=ci(:,1)';
        case 'fitbinomialci'
            pd=fitdist(y,'Binomial');
            ymean=pd.mean;
            ci=paramci(pd,alpha);
            yci=ci(:,2)';
        case 'fit95percentile'
            pd=fitdist(y,'Normal');
            ymean=pd.icdf(0.5);
            yci=pd.icdf([0.025 0.975]);
        otherwise
            warning(['Unknown CI type ' type]);
    end
catch ME
    disp(['Error in CI computation: ' ME.message ' ...skipping']);
    yci=repmat([NaN NaN],size(y,2));
end
end

author: piermorel

repo: gramm

save_path: github-repos/MATLAB/piermorel-gramm


text:
% To run unit tests:
%
% testCase = TestBspline;
% res = run(testCase);
%
classdef TestBspline < matlab.unittest.TestCase
    %
    % Test the spline scripts
    %

    properties
       cpts
    end

    methods (TestClassSetup)

       function setup(this)
            cpts = [0 2;
                        0 4;
                        0 6;
                        0 8;
                        0 10];

            cpts = cpts + 10*rand*randn(size(cpts));
            this.cpts = cpts;
       end

    end

    methods (Test)

        function recover(this)
            nland = size(this.cpts,1);
            sval = bspline_gen_s(nland);
            X = bspline_eval(sval,this.cpts);

            % Test fit of control points
            P = bspline_fit(sval,X,nland);
            if aeq(P,this.cpts)
               fprintf(1,'Landmarks recovered exactly! \n');
            else
               fprintf(1,'FAILED TO RECOVER landmarks \n');
            end

            % plot the spline
            figure(1);
            clf
            hold on
            x = X(:,1);
            y = X(:,2);
            plot(x,y,'b-');
            plot(this.cpts(:,1),this.cpts(:,2),'r.','MarkerSize',16);

        end
    end

end

author: brendenlake

repo: BPL

save_path: github-repos/MATLAB/brendenlake-BPL


text:
function T= opt_substruct(S, fld_list)
%OPT_SUBSTRUCT - Create substruct of specified fields
%
%Synopsis:
%  T= opt_substruct(S, <FLD_LIST>)
%
%Arguments:
%  S:  STRUCT - Source, from which fields are taken
%  FLD_LIST: CELL of CHAR - Names of the fields that should be copied from
%      S into the newly generated substruct T
%
%Returns:
%  T: STRUCT - Substruct of S with specified fields

% Let us be gracious:
if ischar(fld_list),
  fld_list= {fld_list};
elseif iscell(fld_list) && ismember(size(fld_list,2), [2 3]),
  fld_list= fld_list(:,1);
end

fldS= fieldnames(S);
[dmy, dmy, fld_idx]= intersect(lower(fld_list), lower(fldS), 'legacy');
fld_list= fldS(fld_idx);

T= [];
for ii= 1:length(fld_list),
  fld= fld_list{ii};
  [T.(fld)]= S.(fld);
end

author: bbci


text:
% navego_example_real_mpu6000: post-processing integration of MPU-6000
% IMU and Ekinox GNSS data.
%
% The main goal is to integrate MPU-6000 IMU and Ekinox-D GNSS measurements.
%
% Sensors dataset was generated driving a car through the streets of
% Turin city (Italy).
%
%   Copyright (C) 2014, Rodrigo Gonzalez, all rights reserved.
%
%   This file is part of NaveGo, an open-source MATLAB toolbox for
%   simulation of integrated navigation systems.
%
%   NaveGo is free software: you can redistribute it and/or modify
%   it under the terms of the GNU Lesser General Public License (LGPL)
%   version 3 as published by the Free Software Foundation.
%
%   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 Lesser General Public License for more details.
%
%   You should have received a copy of the GNU Lesser General Public
%   License along with this program. If not, see
%   <http://www.gnu.org/licenses/>.
%
% References:
%
%   SBG Systems. SBG Ekinox-D High Accuracy Inertial System Brochure,
% Tactical grade MEMS Inertial Systems, v1.0. February 2014.
%
%   InvenSense Inc. MPU-6000/MPU-6050 Product Specification. Document
% Number: PS-MPU-6000A-00. Revision: 3.4. Release Date: 08/19/2013.
%
%   R. Gonzalez and P. Dabove. Performance Assessment of an Ultra Low-Cost
% Inertial Measurement Unit for Ground Vehicle Navigation. Sensors 2019,
% 19(18). https://www.mdpi.com/530156.
%
% Version: 006
% Date:    2021/12/08
% Author:  Rodrigo Gonzalez <rodralez@frm.utn.edu.ar>
% URL:     https://github.com/rodralez/navego

% NOTE: NaveGo assumes that IMU is aligned with respect to body-frame as
% X-forward, Y-right and Z-down.
%
% NOTE: NaveGo assumes that yaw angle (heading) is positive clockwise.

clc
close all
clear
matlabrc

addpath ../../ins/
addpath ../../ins-gnss/
addpath ../../conversions/
addpath ../../performance-analysis/
addpath ../../misc/
addpath ../../plot/

navego_print_version;

fprintf('\nNaveGo: starting real INS/GNSS integration... \n')

%% PARAMETERS

% Comment any of the following parameters in order to NOT execute a
% particular portion of code

INS_GNSS = 'ON';
PLOT     = 'ON';

if (~exist('INS_GNSS','var')), INS_GNSS = 'OFF'; end
if (~exist('PLOT','var')),     PLOT     = 'OFF'; end

%% CONVERSION CONSTANTS

G =  9.80665;       % Gravity constant, m/s^2
G2MSS = G;          % g to m/s^2
MSS2G = (1/G);      % m/s^2 to g

D2R = (pi/180);     % degrees to radians
R2D = (180/pi);     % radians to degrees

KT2MS = 0.514444;   % knot to m/s
MS2KMH = 3.6;       % m/s to km/h

%% REF DATA

% Reference dataset was obtained by processing Ekinox IMU and Ekinox GNSS
% with tighly-coupled integration by Inertial Explorer software package.

fprintf('NaveGo: loading reference data... \n')

load ref

%% MPU-6000 IMU

fprintf('NaveGo: loading MPU-6000 IMU data... \n')

load mpu6000_imu

%% EKINOX GNSS

fprintf('NaveGo: loading Ekinox GNSS data... \n')

load ekinox_gnss

% ekinox_gnss contains the lever arm with respect to Ekinox IMU.
% ekinox_gnss.larm has to be changed for MPU-6000 IMU.
ekinox_gnss.larm = [-0.369, 0.0, -0.219]';

ekinox_gnss.eps = mean(diff(mpu6000_imu.t)) / 2; %  A rule of thumb for choosing eps.

%% NAVIGATION TIME

to = (ref.t(end) - ref.t(1));

fprintf('NaveGo: navigation time is %.2f minutes or %.2f seconds. \n', (to/60), to)

%% INS/GNSS INTEGRATION

if strcmp(INS_GNSS, 'ON')

    fprintf('NaveGo: processing INS/GNSS integration... \n')

    % Execute INS/GNSS integration
    % ---------------------------------------------------------------------
    nav_mpu6000 = ins_gnss(mpu6000_imu, ekinox_gnss, 'dcm');
    % ---------------------------------------------------------------------

    save nav_mpu6000.mat nav_mpu6000
else

    load nav_mpu6000
end

%% TRAVELED DISTANCE

distance = gnss_distance (nav_mpu6000.lat, nav_mpu6000.lon);

fprintf('NaveGo: distance traveled by the vehicle is %.2f meters or %.2f km. \n', distance, distance/1000)

%% ANALYSIS OF PERFORMANCE FOR A CERTAIN PART OF THE INS/GNSS DATASET

tmin_rmse = ref.t(1);
tmax_rmse = ref.t(end);

% Sincronize REF data to tmin and tmax
idx  = find(ref.t > tmin_rmse, 1, 'first' );
fdx  = find(ref.t < tmax_rmse, 1, 'last' );
if(isempty(idx) || isempty(fdx))
    error('ref: empty index')
end

ref.t       = ref.t    (idx:fdx);
ref.roll    = ref.roll (idx:fdx);
ref.pitch   = ref.pitch(idx:fdx);
ref.yaw     = ref.yaw  (idx:fdx);
ref.lat     = ref.lat  (idx:fdx);
ref.lon     = ref.lon  (idx:fdx);
ref.h       = ref.h    (idx:fdx);
ref.vel     = ref.vel  (idx:fdx, :);

%% INTERPOLATION OF INS/GNSS DATASET

% INS/GNSS estimates and GNSS data are interpolated according to the
% reference dataset.

[nav_i,  ref_n] = navego_interpolation (nav_mpu6000, ref);
[gnss_i, ref_g] = navego_interpolation (ekinox_gnss, ref);

%% NAVIGATION RMSE

rmse_v = print_rmse (nav_i, gnss_i, ref_n, ref_g, 'MPU-6000 INS/GNSS');

%% RMSE TO CVS FILE

csvwrite('mpu6000.csv', rmse_v);

%% NAVIGATION DATA TO CSV FILE

fprintf('\n');
navego_nav2csv(nav_mpu6000);

%% PLOTS

if (strcmp(PLOT,'ON'))

   navego_plot_main (ref, ekinox_gnss, nav_mpu6000, gnss_i, nav_i, ref_g, ref_n)
end

%% PERFORMANCE ANAYSYS OF THE KALMAN FILTER

fprintf('\nNaveGo: Kalman filter performance analysis...\n')

kf_analysis (nav_mpu6000)

author: rodralez

repo: NaveGo

save_path: github-repos/MATLAB/rodralez-NaveGo


text:
function model = addCOBRAConstraints(model, idList, d, varargin)
% Add a constraint to the model in the form c1 * v(id1) + c2 * v(id2) * ... cj * v(idj) dsense d
% where c is a vector with coefficients for each element in idList
% (default 1 for each reaction), dsense is one of lower than ('L', default), greater than
% ('G'), or equal ('E'), and d is a value.
% USAGE:
%    model = addCOBRAConstraints(model, idList, d, varargin)
%
% INPUTS:
%    model:         model structure
%    idList:        cell array of ids either from either the rxns, or the
%                   evars vectors. Can also be a a double vector of indices (in which
%                   case indices in evars have to be set off by
%                   numel(rxns).
%    d:             The right hand side of the C*v <= d constraint (or a
%                   vector, for multiple simultaneous addition;
%
%    varargin:
%                   * c:                the coefficients to use with one entry per
%                                       reaction of a constraint for multiple constraints, a matrix
%                                       (default: 1 for each element indicated in idList)
%                   * dsense:           the constraint sense ('L': <= ,
%                                       'G': >=, 'E': =), or a vector for multiple constraints
%                                       (default: ('L'))
%                   * ConstraintID:     the Name of the constraint. by
%                                       (default: 'ConstraintXYZ' with XYZ being the initial
%                                       position in the ctrs vector)
%                                       or a cell array of Strings for
%                                       multiple Constraints
%                   * checkDuplicates:  check whether the constraint already
%                                       exists, and if it does, don't add it (default: false)
%
% OUTPUT:
%    modelConstrained:      The constrained model containing the added
%                           constraints in the respective fields:
%                            * `.C` - The constraint matrix containing coefficients for reactions
%                            * `.ctrs` - The constraint IDs
%                            * `.dsense` - The constraint senses
%                            * `.d` the constraint right hand side values
%                            *  Optional: `D`, the matrix conatining coefficients for additional variables.
%
%
% EXAMPLE:
%    Add a constraint that leads to EX_glc and EX_fru not carrying a
%    combined flux higher than 5
%    model = addCOBRAConstraints(model, {'EX_glc','EX_fru'}, 5)
%    Assume Reaction 4 to be 2 A -> D and reaction 5 being A -> F. Create a
%    constraint that requires that the two reactions remove at least 4 units of A:
%    model = addCOBRAConstraints(model, model.rxns(4:5), 4, 'c', [2 1], 'dsense', 'G')
%
% NOTE:
%    This function will, if not present create the `C`, `ctrs`, `dsense`, `d` and,
%    if additional variables (i.e. the `E` matrix) is present, the `D`
%    field as defined in the model field definitions. The `C` field containts
%    the constraint coefficients referring to reactions, while the `D`
%    field contains the Constraint coefficients referring to additional
%    variables. `d` represents the right hand side of constraints and
%    `dsense` the sense (lower -'L' , equal - 'E' or greater than - 'G') of the constraint.
%    `ctrs` stores unique IDs for the constraints.
%
% Author: Thomas Pfau, Nov 2017

if ischar(idList) %This is an individual element.
    idList = {idList};
end

if iscell(idList)
    [pos,pres] = getIDPositions(model,idList,'rxns');
    if any(~pres)
        missingreactions = idList(~pres);
        error('The following ids were not found in the model:\n%s\nNo Constraint was added',strjoin(missingreactions,', '));
    end
    idList = pos;
end

multiAdd = false;
if ~all(size(idList) > 1) %if this is true, its multiple rows...
    if (length(d) > 1) %This should be a multiAdd
        dim = length(d) == size(idList);
        multiAdd = true;
        if ~all(dim)
            %make sure this is the right orientation
             cdim = find(dim);
             if cdim ~= 1
                 idList = idList';
             end
        else
            error('d has to be either a single value or a vector of doubles');
        end
    else
        %there is only one element in d and idList has a dimension of size
        %1. So we need to make sure, that idList is a row vector.
        if size(idList,1) > size(idList,2)
            idList = idList';
        end
        if ~(length(idList) == length(unique(idList)))
            error('There were duplicate reaction IDs or positions provided. No Constraint will be added.');
        end
    end
else
    multiAdd = true;
end

defaultcoefficients = ones(numel(d),numel(idList));

defaultcsense = repmat('L',numel(d),1);
if ~multiAdd
    defaultConstraintName = getConstraintName(model,1);
else
    defaultConstraintName = getConstraintName(model,length(d));
end

parser = inputParser();
parser.addRequired('model',@isstruct);
parser.addRequired('idList',@(x) isnumeric(x));
parser.addRequired('d',@isnumeric);
parser.addParamValue('c',defaultcoefficients,@(x) isnumeric(x) && (multiAdd || length(x) == length(idList)));
parser.addParamValue('dsense',defaultcsense, @ischar );
parser.addParamValue('ConstraintID',defaultConstraintName,@(x) ischar(x) || iscell(x) );
parser.addParamValue('checkDuplicates',false,@(x) islogical(x) || isnumeric(x) );

parser.parse(model,idList,d,varargin{:});

c = parser.Results.c;
d = columnVector(parser.Results.d);
dsense = columnVector(parser.Results.dsense);
ctrID = parser.Results.ConstraintID;
if ischar(ctrID)
    ctrID = {ctrID};
else
    ctrID = columnVector(ctrID);
end
checkDuplicate = parser.Results.checkDuplicates;

if ~isfield(model,'evars')
    varspresent = false;
    ConstraintFields = {'C','d','dsense','ctrs'};
else
    varspresent = true;
    ConstraintFields = {'C','D','d','dsense','ctrs'};
end
fieldsPresent = ismember(ConstraintFields,fieldnames(model));
if ~all(fieldsPresent)
    if any(fieldsPresent)
        for i = 1:numel(ConstraintFields)
            if ~isfield(model,ConstraintFields{i})
                model = createEmptyField(model,ConstraintFields{i});
            else
                if ~isempty(model.(ConstraintFields{i}))
                    error('Inconsistent Field sizes detected. Expected an empty field but %s was non empty',ConstraintFields{i})
                else
                    model = createEmptyField(model,ConstraintFields{i}); %Replace the existing to make sure the dimensions are correct.
                end
            end
        end
    else
        model = createEmptyFields(model,ConstraintFields);
    end
end

%Now, we have all fields.
%Check for duplicates:
if any(ismember(model.ctrs,ctrID))
    error('Constraints with the following IDs already exist.\n%s\n', strjoin(model.ctrs(ismember(model.ctrs,ctrID)),', '));
end
%And if multiAdd also check for duplicates in the input
if multiAdd && numel(ctrID) > numel(unique(ctrID))
    [unq,~,orig] = unique(ctrID);
    dupreacs = hist(orig,numel(unq)) > 1;
    error('Input contained the following duplicated IDs:\n %s ', strjoin(unq(dupreacs),', '));
end

if checkDuplicate && multiAdd
    %If we check for duplicates, and don't want to add them, we will first
    %filter them from the input, but only if we add multiple things.
    [sorted,order] = sort(idList,2);
    sortedC = c;
    sorted = repmat(sorted,size(c,1),1);
    for i = 1:size(idList,1)
        sortedC(i,:) = c(i,order(i,:));
    end
    %Now, concatenate all inputs (except names)
    toCompare = [sorted,sortedC,dsense]; %This will convert the dsense into doubles which is fine to get uniques.
    [~,pos] = unique(toCompare,'rows');
    %Now lets remove anything thats duplicated.
    d = d(pos,:);
    c = c(pos,:);
    dsense = dsense(pos,:);
    ctrID = ctrID(pos,:);
end

%In case nothing is left.
if isempty(ctrID)
    return
end

%Also check for duplicates in the C Matrix.
if varspresent
    mixrows = sparse(size(c,1),size(model.E,2));
end
constRow = sparse(size(c,1),size(model.C,2));
duppedRows = false(size(idList,1),1);
nRxns = numel(model.rxns);
rxnpos = idList <= nRxns;
varpos = idList > nRxns;

for i = 1:size(c,1)
    constRow(i,idList(rxnpos)) = c(i,rxnpos);
    if varspresent
         mixrows(i,idList(varpos)-nRxns) = c(i,varpos);
    end
    cRow = constRow(i,:);
    if checkDuplicate
        dupRows = all(model.C == cRow(ones(size(model.C,1),1),:),2);
    end
    if varspresent
        cmRow = mixrows(i,:);
        if checkDuplicate
            vdupRows = all(model.E == cmRow(ones(size(model.E,1),1),:),2);
            dupRows = dupRows & vdupRows;
        end
    end
    if checkDuplicate
        duppedRows(i) = any(dupRows) && (model.dsense(dupRows) == dsense(i)) && (model.d(dupRows) == d(i));
    end
end

if any(duppedRows) && checkDuplicate
    warning('Constraint(s) not added, because it already exists with ID: %s',strjoin(model.ctrs(dupRows)));
    if ~multiAdd
        return
    else
        ctrID(duppedRows) = [];
        d(duppedRows) = [];
        constRow(duppedRows,:) = [];
        dsense(duppedRows) = [];
    end
end

model.C = [model.C;constRow];
model.d = [model.d;d];
model.dsense = [model.dsense;dsense];
originalSize = numel(model.ctrs);
model.ctrs = [model.ctrs;ctrID];
if varspresent
    model.D = [model.D;mixrows];
end

model = extendModelFieldsForType(model,'ctrs','originalSize',originalSize);

end

function name = getConstraintName(model, count)
%Get a unique, not yet used constraint name.
if ~isfield(model,'ctrs')
    name = strcat('Constraint',cellfun(@num2str,num2cell(1:count)','UniformOutput',false));
    return
else
    name = cell(count,1);
    name(:) = {''}; %Need to initialize, otherwise we get a mismatch during ismember.
    constraintsCreated = 0;
    i = size(model.C,1) + 1;
    cname = ['Constraint' num2str(i)];
    while constraintsCreated < count
        while any(ismember(model.ctrs,cname)) || any(ismember(name,cname))
            i = i + 1;
            cname = ['Constraint' num2str(i)];
        end
        constraintsCreated = constraintsCreated + 1;
        name{constraintsCreated} = cname;
    end
end
end

author: opencobra

repo: cobratoolbox


text:
%
% Perform an iteration of MCMC
% for re-sampling the relations, which
% blends type-level and token-levle variables
%
function mcmc_iter_relations(MH,M,lib)

    % relations
    for sid=2:M.ns
        MH.mh_relation(sid,M,lib);
    end

    % evaluation token spots
    for sid=1:M.ns
       if strcmp(M.S{sid}.R.type,'mid')
          MH.mh_eval_spot_token(sid,M,lib);
       end
    end

end


text:
function startup()
% startup()
% --------------------------------------------------------
% R-FCN implementation
% Modified from MATLAB Faster R-CNN (https://github.com/shaoqingren/faster_rcnn)
% Copyright (c) 2016, Jifeng Dai
% Licensed under The MIT License [see LICENSE for details]
% --------------------------------------------------------

    curdir = fileparts(mfilename('fullpath'));
    addpath(genpath(fullfile(curdir, 'utils')));
    addpath(genpath(fullfile(curdir, 'functions')));
    addpath(genpath(fullfile(curdir, 'bin')));
    addpath(genpath(fullfile(curdir, 'experiments')));
    addpath(genpath(fullfile(curdir, 'imdb')));

    mkdir_if_missing(fullfile(curdir, 'datasets'));

    mkdir_if_missing(fullfile(curdir, 'external'));

    caffe_path = fullfile(curdir, 'external', 'caffe', 'matlab');
    if exist(caffe_path, 'dir') == 0
        error('matcaffe is missing from external/caffe/matlab; See README.md');
    end
    addpath(genpath(caffe_path));

    mkdir_if_missing(fullfile(curdir, 'imdb', 'cache'));

    mkdir_if_missing(fullfile(curdir, 'output'));

    mkdir_if_missing(fullfile(curdir, 'models'));

    fprintf('rfcn startup done\n');
end

repo: R-FCN


text:
function test_failed = test_libltfat_fold(varargin)

test_failed = 0;

fprintf(' ===============  %s ================ \n',upper(mfilename));

definput.flags.complexity={'double','single'};
[flags]=ltfatarghelper({},definput,varargin);
dataPtr = [flags.complexity, 'Ptr'];

Larr =     [10 11  10   10  101 101];
Lfoldarr = [3  16   3    3    1   2];
shiftarr = [0   0  -5 -101   -3  11];

for do_complex = 0:1
    complexstring = '';
    if do_complex, complexstring = 'complex'; end
    funname = makelibraryname('fold_array',flags.complexity,do_complex);

    for lId = 1:numel(Larr)
        L     = Larr(lId);
        Lfold = Lfoldarr(lId);
        shift = shiftarr(lId);

        if do_complex
            z = (1:max(L,Lfold)) + 1i*(max(L,Lfold):-1:1);
            z(L+1:end) = 0;
            zi = complex2interleaved(z);
            zout = complex2interleaved(randn(1,Lfold) + 1i*randn(1,Lfold));
        else
            z = (1:max(L,Lfold));
            z(L+1:end) = 0;
            zi = z;
            zout = randn(1,Lfold);
        end

        ziPtr = libpointer(dataPtr,zi);
        zoutPtr = libpointer(dataPtr,zout);

        periods = ceil(L/Lfold);
        fext = postpad(z,periods*Lfold);
        ffoldtrue = circshift(sum(reshape(fext, Lfold, periods),2).',[0,shift]);
        ffold2ndtrue = sum(reshape(circshift(fext,[0,shift]), Lfold, periods),2).';

        errTrues = norm(ffoldtrue - ffold2ndtrue);

        status = calllib('libltfat',funname,ziPtr,L,shift,Lfold,zoutPtr);

        if do_complex
            res = norm(ffoldtrue - interleaved2complex(zoutPtr.Value));
        else
            res = norm(ffoldtrue - zoutPtr.Value);
        end

        [test_failed,fail]=ltfatdiditfail(res+status,test_failed,0);
        fprintf(['FOLD OP L:%3i, Lfold:%3i, shift:%3i, %s %s %s %s\n'],L,Lfold,shift,flags.complexity,complexstring,ltfatstatusstring(status),fail);

        status = calllib('libltfat',funname,ziPtr,L,shift,Lfold,ziPtr);

        if do_complex
            res = norm(ffoldtrue - postpad(interleaved2complex(ziPtr.Value),Lfold));
        else
            res = norm(ffoldtrue - ziPtr.Value(1:Lfold));
        end

        [test_failed,fail]=ltfatdiditfail(res+status,test_failed,0);
        fprintf(['FOLD IP L:%3i, Lfold:%3i, shift:%3i, %s %s %s %s\n'],L,Lfold,shift,flags.complexity,complexstring,ltfatstatusstring(status),fail);
    end
end

%interleaved2complex(zoutPtr.Value)

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function maxTMap=getMaxTVal(numberOfTMaps,tMapRoot,tMapStartIndex)
% function maxTMap=getMapTVal(numberOfTMaps,tMapRoot)
% Finds the MAX of a set of T maps.
% Don't know what to do to make a retinotopy map...the best way, maybe
  % is to find the MAX of T maps of contrasts for each different phase

    if (~exist('tMapRoot','var'))
      tMapRoot='spmT_000'; % Can't have more than 9 different phases for
                           % now
    end

    if (~exist('tMapStartIndex','var'))
      tMapStartIndex=1; %

    end
    if (numberOfTMaps>8)
      error('Cannot have more than 8 TMaps at the moment');
    end

    % Load in the first one just to get the size
    firstMapVol=spm_vol([tMapRoot,int2str(tMapStartIndex)]);
    firstMap=spm_read_vols(firstMapVol);
    [x y z]=size(firstMap);
    counter=1;
    fprintf('\nSize is %d x %d x%d',x,y,z);

    % First TMap should be spm_T0001
    for thisTmap=1:(numberOfTMaps)
      fileName=[tMapRoot,int2str(thisTmap+tMapStartIndex-1)];

        mapVol(counter)=spm_vol(fileName);
        mapData(:,:,:,counter)=spm_read_vols(mapVol(counter));
        disp(thisTmap);
        counter=counter+1;
    end

    % Now for each plane we want to find the max
    maxTMap=max(mapData,[],4);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%FILE2DSET Load and convert dataset from disk file
%
%  A = FILE2DSET(NAME,M,N)
%
% The function was previously named PRDATASET. It has to be renamed for the
% renaming of DATASET into PRDATASET in PRTOOLS version 5. It may also be
% called by A = PRDATASET(NAME,M,N), both, in version 4 as well as in
% version 5.
%
% The dataset given in NAME is loaded from a .mat file and converted
% to the current PRTools definition. Objects and features requested
% by the index vectors M and N are returned.
%
% See PRDATA for loading arbitrary data into a PRTools dataset.
% See PRDATASETS for an overview of datasets.

% Copyright: R.P.W. Duin, duin@ph.tn.tudelft.nl
% Faculty of Applied Sciences, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands

% $Id: prdataset.m,v 1.3 2009/12/21 10:35:12 duin Exp $

% function a = prdataset(name,M,N,no_message)

function a = file2dset(name,varargin)

	persistent FIRST;
	if isempty(FIRST), FIRST = 1; end
  [M,N,no_message] = setdefaults(varargin,[],[],0);
	if exist([name '.mat']) ~= 2
		error([newline '---- Dataset ''' name ''' not available ----'])
	end

	s = warning;
	warning off
	b = load(name);
	warning(s);
	% Try to find out what data we actually loaded:
	names = fieldnames(b);
	eval(['a = b.' names{1} ';']);

	if ~isdataset(a) & ndims(a) > 2
		% We loaded an image?
		a = im2feat(a);
		prwarning(2,'Assumed that a feature image is loaded')
	else
		a = prdataset(a);
	end
	[m,k] = size(a);
	% Maybe a subset of the features is requested:
	if ~isempty(N)
		if min(size(N)) ~= 1
			error([newline 'Feature indices should be stored in vector'])
		end
		if max(N) > k
			error([newline 'Dataset has not the requested features'])
		end
		a = a(:,N);
	end

	% Maybe a subset of the objects are requested:
	if ~isempty(M)
		if min(size(M)) ~= 1
			error([newline 'Object indices should be stored in vector'])
		end
		if max(M) > m
			error([newline 'Dataset has not the requested objects'])
		end
		a = a(M,:);
	end

	if FIRST & ~no_message
		disp(' ')
		disp('*** You are using one of the datasets distibuted by PRTools. Most of them')
		disp('*** are publicly available on the Internet and converted to the PRTools')
		disp('*** format. See the Contents of the datasets directory (HELP PRDATASETS)')
		disp('*** or the load routine of the specific dataset to retrieve its source.')
		disp(' ')
		FIRST = 0;
	end
	return

author: marianux

repo: ecg-kit

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


text:
%
% These are the functions to import and export an SBML model into a
% MATLAB_SBML structure and back to an SBML file.
%
% FUNCTIONS include:
%
% TranslateSBML(varargin)
%
% - translates a sbml file into a matlab structure
%
% OutputSBML(varargin)
%
% - translates a appropriate matlab structure back into sbml and writes
%   out the file
%
% addLevelVersion.m
%
% - a script used by TranslateSBML to add a level and version field to all
%   parts of a structure
%
% applyUserValidation.m
%
% - a script that can be customized by the user to add additional
%   validation to the isSBML_Model function
%
% CheckAndConvert.m
%
% - a script used by TranslateSBML to change some mathematical function names
%   to those used by MATLAB
%
% ConvertFormulaToMathML.m
%
% - a script used by OutputSBML to change some mathematical function names
%   to those used by MathML
%
% getDefaultValues.m
%
% - a script used to retreive the default values for a given SBMLStructure
%
% getStructureFieldnames.m
%
% - a script used to retreive the fieldnames for a given SBMLStructure
%
% getValueType.m
%
% - a script used to retreive the value types for a given SBMLStructure
%
% isSBML_Model.m
%
% - a script used by OutputSBML to check that a structure is an appropriate
%   MATLAB_SBML structure for conversion to SBML
%
% isoctave.m
%
% - a script to determine if octave or matlab is being used
%

% Description : This is the binding to translate
%				 sbml models into a MATLAB structure
% Author(s)   : SBML Team <sbml-team@googlegroups.com>
%
% This file is part of libSBML.  Please visit http://sbml.org for more
% information about SBML, and the latest version of libSBML.
%
% Copyright (C) 2013-2018 jointly by the following organizations:
%     1. California Institute of Technology, Pasadena, CA, USA
%     2. EMBL European Bioinformatics Institute (EMBL-EBI), Hinxton, UK
%     3. University of Heidelberg, Heidelberg, Germany
%
% Copyright (C) 2009-2013 jointly by the following organizations:
%     1. California Institute of Technology, Pasadena, CA, USA
%     2. EMBL European Bioinformatics Institute (EMBL-EBI), Hinxton, UK
%
% Copyright (C) 2006-2008 by the California Institute of Technology,
%     Pasadena, CA, USA
%
% Copyright (C) 2002-2005 jointly by the following organizations:
%     1. California Institute of Technology, Pasadena, CA, USA
%     2. Japan Science and Technology Agency, Japan
%
% This library is free software; you can redistribute it and/or modify it
% under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation.  A copy of the license agreement is provided
% in the file named "LICENSE.txt" included with this software distribution
% and also available online as http://sbml.org/software/libsbml/license.html
%
% The original code contained here was initially developed by:
%
%      Sarah Keating
%      Science and Technology Research Centre
%      University of Hertfordshire
%      Hatfield, AL10 9AB
%      United Kingdom
%
%      http://www.sbml.org
%      mailto:sbml-team@googlegroups.com
%
% Contributor(s):
%
%
%


text:
function a = prolate ( alpha, n )

%*****************************************************************************80
%
%% PROLATE returns the PROLATE matrix.
%
%  Formula:
%
%    If ( I == J )
%      A(I,J) = 2 * ALPHA
%    else
%      K = abs ( I - J ) + 1
%      A(I,J) = sin ( 2 * pi * ALPHA * K ) / ( pi * K )
%
%  Example:
%
%    N = 5, ALPHA = 0.25
%
%     0.5        0.0       -0.106103   0.0        0.0636620
%     0.0        0.5        0.0       -0.106103   0.0
%    -0.106103   0.0        0.5        0.0       -0.106103
%     0.0       -0.106103   0.0        0.5        0.0
%     0.0636620  0.0       -0.106103   0.0        0.5
%
%  Properties:
%
%    A is symmetric: A' = A.
%
%    Because A is symmetric, it is normal.
%
%    Because A is normal, it is diagonalizable.
%
%    A is persymmetric: A(I,J) = A(N+1-J,N+1-I).
%
%    A is centrosymmetric: A(I,J) = A(N+1-I,N+1-J).
%
%    A is Toeplitz: constant along diagonals.
%
%    If 0 < ALPHA < 0.5, then
%      A is positive definite,
%      the eigenvalues of A are distinct,
%      the eigenvalues lie in (0,1) and cluster around 0 and 1.
%
%    The family of matrices is nested as a function of N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    James Varah,
%    The Prolate Matrix,
%    Linear Algebra and Applications,
%    Volume 187, pages 269-278, 1993.
%
%  Parameters:
%
%    Input, real ALPHA, the parameter.
%
%    Input, integer N, the order of A.
%
%    Output, real A(N,N), the matrix.
%
  a = zeros ( n, n );

  for i = 1 : n

    for j = 1 : n

      if ( i == j )
        a(i,j) = 2.0 * alpha;
      else
        k = abs ( i - j ) + 1;
        angle = 2.0 * pi * alpha * k;
        a(i,j) = sin ( angle ) / ( pi * k );
      end

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function PLOT_Estimator(Time, Trace_Cov, State_Error)

figure(6);
plot(Time,Trace_Cov)
title('Trace of the covariance matrix')
xlabel('Time (s)')
ylabel('Trace of Covariance')

figure(7)

subplot(4,1,1)
plot(Time, State_Error(1,:))
title('EKF Error: X position', 'FontSize',16)
xlabel('Time (s)')
ylabel('Distance (m)')

subplot(4,1,2)
plot(Time, State_Error(2,:))
title('EKF Error: Y position', 'FontSize',16)
xlabel('Time (s)')
ylabel('Distance (m)')

subplot(4,1,3)
plot(Time, State_Error(3,:))
title('EKF Error: Orientation', 'FontSize',16)
xlabel('Time (s)')
ylabel('Angle (rad)')

subplot(4,1,4)
plot(Time, State_Error(4,:))
title('EKF Error: Cab-Trailer Angle', 'FontSize',16)
xlabel('Time (s)')
ylabel('Angle (rad)')

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function a = r8sd_dif2 ( m, n, ndiag, offset )

%*****************************************************************************80
%
%% R8SD_DIF2 returns the DIF2 matrix in R8SD format.
%
%  Example:
%
%    N = 5
%
%    2 -1  .  .  .
%   -1  2 -1  .  .
%    . -1  2 -1  .
%    .  . -1  2 -1
%    .  .  . -1  2
%
%  Properties:
%
%    A is banded, with bandwidth 3.
%
%    A is tridiagonal.
%
%    Because A is tridiagonal, it has property A (bipartite).
%
%    A is a special case of the TRIS or tridiagonal scalar matrix.
%
%    A is integral, therefore det ( A ) is integral, and
%    det ( A ) * inverse ( A ) is integral.
%
%    A is Toeplitz: constant along diagonals.
%
%    A is symmetric: A' = A.
%
%    Because A is symmetric, it is normal.
%
%    Because A is normal, it is diagonalizable.
%
%    A is persymmetric: A(I,J) = A(N+1-J,N+1-I).
%
%    A is positive definite.
%
%    A is an M matrix.
%
%    A is weakly diagonally dominant, but not strictly diagonally dominant.
%
%    A has an LU factorization A = L * U, without pivoting.
%
%      The matrix L is lower bidiagonal with subdiagonal elements:
%
%        L(I+1,I) = -I/(I+1)
%
%      The matrix U is upper bidiagonal, with diagonal elements
%
%        U(I,I) = (I+1)/I
%
%      and superdiagonal elements which are all -1.
%
%    A has a Cholesky factorization A = L * L', with L lower bidiagonal.
%
%      L(I,I) =    sqrt ( (I+1) / I )
%      L(I,I-1) = -sqrt ( (I-1) / I )
%
%    The eigenvalues are
%
%      LAMBDA(I) = 2 + 2 * COS(I*PI/(N+1))
%                = 4 SIN^2(I*PI/(2*N+2))
%
%    The corresponding eigenvector X(I) has entries
%
%       X(I)(J) = sqrt(2/(N+1)) * sin ( I*J*PI/(N+1) ).
%
%    Simple linear systems:
%
%      x = (1,1,1,...,1,1),   A*x=(1,0,0,...,0,1)
%
%      x = (1,2,3,...,n-1,n), A*x=(0,0,0,...,0,n+1)
%
%    det ( A ) = N + 1.
%
%    The value of the determinant can be seen by induction,
%    and expanding the determinant across the first row:
%
%      det ( A(N) ) = 2 * det ( A(N-1) ) - (-1) * (-1) * det ( A(N-2) )
%                = 2 * N - (N-1)
%                = N + 1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 July 2000
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Robert Gregory, David Karney,
%    A Collection of Matrices for Testing Computational Algorithms,
%    Wiley, 1969,
%    ISBN: 0882756494,
%    LC: QA263.68
%
%    Morris Newman, John Todd,
%    Example A8,
%    The evaluation of matrix inversion programs,
%    Journal of the Society for Industrial and Applied Mathematics,
%    Volume 6, Number 4, pages 466-476, 1958.
%
%    John Todd,
%    Basic Numerical Mathematics,
%    Volume 2: Numerical Algebra,
%    Birkhauser, 1980,
%    ISBN: 0817608117,
%    LC: QA297.T58.
%
%    Joan Westlake,
%    A Handbook of Numerical Matrix Inversion and Solution of
%    Linear Equations,
%    John Wiley, 1968,
%    ISBN13: 978-0471936756,
%    LC: QA263.W47.
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns.
%
%    Input, integer NDIAG, the number of diagonals that are stored.
%    NDIAG must be at least 2.
%
%    Input, integer OFFSET(NDIAG), the offsets for the diagonal
%    storage.  It is simply assumed that OFFSET(1) = 0 and OFFSET(2) = 1.
%
%    Output, real A(N,NDIAG), the matrix.
%
  a = zeros(n,ndiag);

  a(1:n,  1) =  2.0;
  a(1:n-1,2) = -1.0;

  return
end


text:
function [elec, grad] = out_fieldtrip_channel(ChannelFile, isIncludeRef)
% OUT_FIELDTRIP_CHANNEL: Converts a channel file into elec/grad structures
%
% USAGE:  [elec, grad] = out_fieldtrip_channel(ChannelFile, isIncludeRef=1)
%         [elec, grad] = out_fieldtrip_channel(ChannelMat)

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

% ===== PARSE INPUT =====
if (nargin < 2) || isempty(isIncludeRef)
    isIncludeRef = 1;
end
if isstruct(ChannelFile)
    ChannelMat  = ChannelFile;
    ChannelFile = [];
else
    ChannelMat = [];
end

% ===== LOAD CHANNEL FILE =====
% Load channel file
if ~isempty(ChannelFile) && isempty(ChannelMat)
    ChannelMat = in_bst_channel(ChannelFile);
end
% Make sure that the channel file is defined
if isempty(ChannelMat)
    error('No channel file available.');
end
% Find MEG and EEG sensors
iEeg = channel_find(ChannelMat.Channel, 'EEG,SEEG,ECOG');
iMeg = channel_find(ChannelMat.Channel, 'MEG');
iRef = channel_find(ChannelMat.Channel, 'MEG REF');
if isIncludeRef
    iMegAll = [iMeg, iRef];
else
    iMegAll = iMeg;
end

% ===== COMPUTE PROJECTOR =====
% Compute selected SSP projectors
if ~isempty(ChannelMat.Projector)
    % Rebuild projector in the expanded form (I-UUt)
    Proj = process_ssp2('BuildProjector', ChannelMat.Projector, [1 2]);
else
    Proj = [];
end

% ===== EEG =====
if ~isempty(iEeg)
    % Create electrode structure
    elec = struct();
    elec.label = {ChannelMat.Channel(iEeg).Name};
    elec.unit  = 'm';
    % Electrode position
    elec.chanpos = zeros(length(iEeg),3);
    for i = 1:length(iEeg)
        if all(size(ChannelMat.Channel(iEeg(i)).Loc) >= [3,1])
            elec.chanpos(i,:) = ChannelMat.Channel(iEeg(i)).Loc(:,1);
        else
            elec.chanpos(i,:) = [0;0;0];
        end
    end
    elec.elecpos = elec.chanpos;
    % Default montage
    elec.tra = eye(length(iEeg));
    % Apply projectors (SSP or ICA)
    if ~isempty(Proj)
        elec.tra = Proj(iEeg,iEeg) * elec.tra;
    end
else
    elec = [];
end

% ===== MEG =====
if ~isempty(iMeg)
    % Average all the intergration points of the various coils
    chantype = cell(1,length(iMegAll));
    for i = 1:length(iMegAll)
        switch (ChannelMat.Channel(iMegAll(i)).Type)
            case 'MEG'
                chantype{i} = 'megaxial';
            case 'MEG MAG'
                chantype{i} = 'megmag';
            case 'MEG GRAD'
                chantype{i} = 'megplanar';
            case 'MEG REF'
                chantype{i} = 'megref';
        end
    end

    % Create sensor structure
    grad = struct();
    grad.label    = {ChannelMat.Channel(iMegAll).Name};
    grad.unit     = 'm';
    grad.chantype = chantype;
    % Channels positions
    grad.chanpos = figure_3d('GetChannelPositions', ChannelMat, iMegAll);
    % Coils positions
    grad.coilpos = [ChannelMat.Channel(iMegAll).Loc]';
    grad.coilori = [ChannelMat.Channel(iMegAll).Orient]';
    % Correspondance channel-coil
    grad.tra = sparse(length(iMegAll), length(grad.coilpos));
    k = 1;
    for i = 1:length(iMegAll)
        % Dealing with the multiple coils and integration points
        grad.chanori(i,:) = ChannelMat.Channel(iMegAll(i)).Orient(:,1)';
        nCoils = size(ChannelMat.Channel(iMegAll(i)).Weight,2);
        grad.tra(i,k+(0:nCoils-1)) = ChannelMat.Channel(iMegAll(i)).Weight;
        k = k + nCoils;
    end
    % Add MegRefCoef (CTF/4D 3rd order gradient compensation)
    if isIncludeRef && ~isempty(iRef)
        % Error: Not all the sensors are selected
        if (size(ChannelMat.MegRefCoef,1) ~= length(iMeg)) || (size(ChannelMat.MegRefCoef,2) ~= length(iRef))
            error('CTF compensation can be used only when using all the MEG sensors.');
        end
        % Apply compensation matrix to ".tra" matrix
        GradComp = eye(length(ChannelMat.Channel));
        GradComp(iMeg,iRef) = -ChannelMat.MegRefCoef;
        grad.tra = GradComp(iMegAll,iMegAll) * grad.tra;
        % grad.tra is later applied to the leadfield in ft_compute_leadfield:
        % lf = sens.tra * lf;
    end
    % Apply projectors (SSP or ICA)
    if ~isempty(Proj)
        grad.tra = Proj(iMegAll,iMegAll) * grad.tra;
    end
else
    grad = [];
end

author: brainstorm-tools


text:
function SetColorbar
cbar = colorbar;
% Dimensions of the colorbar
% cpos = get(cbar,'position');
% cpos(3) = cpos(3)/4 ;   % Reduce the width of colorbar by half
% cpos(2) = cpos(2)+.1 ;
% cpos(4) = cpos(4)-.2 ;
%set(cbar,'Position',cpos) ;
brighten(0.5);

% Title of the colorbar
%set(get(cbar,'title'),'string','VAL');
%locate = get(cbar,'title');
%tpos = get(locate,'position');
%tpos(3) = tpos(3)+5. ;
%set(locate,'pos',tpos);

% Setting the values on colorbar
%
% get the color limits
clim = caxis;
ylim(cbar,[clim(1) clim(2)]);
numpts = 24 ;    % Number of points to be displayed on colorbar
kssv = linspace(clim(1),clim(2),numpts);
set(cbar,'YtickMode','manual','YTick',kssv); % Set the tickmode to manual
for i = 1:numpts
    imep = num2str(kssv(i),'%+3.2E');
    vasu(i) = {imep} ;
end
set(cbar,'YTickLabel',vasu(1:numpts),'fontsize',9);


text:
function [OUT, lags] = CrossCorrelation(X,Y,lag,do_plot,Yname)
% =======================================================================
% Computes the cross-correlation between the vector X (Tx1) and each
% column of the matrix Y (TxN). If specified, it plots it
% =======================================================================
% [OUT lags] = CrossCorrelation(X,Y,lag,chart,row,col,)
% -----------------------------------------------------------------------
% INPUT
%   - X   = vector (Tx1) of interest [double]
%   - Y   = matrix (TxN) to correlate X with [double]
%   - lag = length of lags to consider [double]
% -----------------------------------------------------------------------
% OPTIONAL INPUT:
%   - do_plot = 1 for plot           [dflt 0]
%   - Yname = array (1xN) name of Y  [dflt Y#]
% -----------------------------------------------------------------------
% OUTPUT
%   - OUT  = matrix (2*lag+1,N) of cross-correlations, where every column
%           is a variable [double]
%   - lags = vector of lags used [double]
% -----------------------------------------------------------------------
% EXAMPLE
%   X = rand(50,1);
%   Y = rand(50,2);
%   [OUT lags] = CrossCorrelation(X,Y,5,1,{'Consumption','Investment'})
% =======================================================================
% VAR Toolbox 3.0
% Ambrogio Cesa-Bianchi
% ambrogiocesabianchi@gmail.com
% March 2012. Updated November 2020
% -----------------------------------------------------------------------

%% Check inputs
%==========================================================================
[a, b] = size(X);
[c, d] = size(Y);
if a ~= c
    disp('Error: vectors must have the same lenght')
    return
elseif b > 1
    disp('Error: X must be a Tx1 vector')
    return
elseif lag >= a-2
    disp('Error: too few observations. Reduce the number of lags')
    return
end
if exist('do_plot','var')==0
    do_plot = 0;
else
end

%% Compute the cross correlation with the max amount of observations
%==========================================================================
% Lead corr(Xt-1,Yt)
m = 1;
for jj=1:lag
    for ii=1:d
        OUT(m,ii) = corr(X(1:end-lag-1+jj),Y(1+lag+1-jj:end,ii));
    end
    m = m+1;
end

% Contemporaneous correlation
for ii=1:d
    OUT(m,ii) = corr(X,Y(:,ii));
end
m = m+1;

% Lag correlation corr(Xt+1,Yt)
for jj=1:lag
    for ii=1:d
        OUT(m,ii) = corr(X(1+jj:end),Y(1:end-jj,ii));
    end
    m = m+1;
end

lags = -lag:1:lag;

%% Plot an istogram with the cross correlations
%==========================================================================
if do_plot

    % If there are no Yname, create Yname (Y1, Y2,...)
    if ~exist('Yname','var')
        aux1 = 'Y';
        Yname(1,d) = {[]};
        for ii=1:d
            Yname(1,ii) = {[aux1 num2str(ii)]};
        end
        clear aux1
    end

    % Dimension of the matrix to plot
    dim=size(OUT);
    ntotlags = dim(1);
    ntotvars = dim(2);
    row = round(sqrt(ntotvars));
    col = ceil(sqrt(ntotvars));

    % Plot
    for jj=1:row*col
        if jj>ntotvars; break; end
        subplot(row,col,jj);
        stem(lags,OUT(:,jj),'*k','fill','LineWidth',1)
        grid off;
        title(Yname(jj));
        set(gca, 'XTick', -lag:lag); % set the number of ticks
    end
end

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox