text:
% --- Update waitbar
function DoUpdateBar(axes_handle,value,total_value)

% prevent matlab from stealing focus when updating bar
set(0,'CurrentFigure',get(axes_handle,'Parent'));
set(gcf,'CurrentAxes',axes_handle);

set(axes_handle,'Visible','on');
% axes(axes_handle);
cla(axes_handle);
axis(axes_handle,[0,total_value,0,1]);
patch([0,value,value,0],[0,0,1,1],'b','FaceAlpha',0.5);
text(0.5,0.5,[num2str((value/total_value)*100, '%02.0f') '%']);
% axis off;

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
clear; close all; clc;

% A = [2 2;2 7];
% A =[1,3;4,2]/4;
A=[2,1;1,2]; % shear
% angle = pi/2; A = [cos(angle) -sin(angle); sin(angle) cos(angle)]; %rotation
% A = [0, 1; 1, 0]; % permutation
% A = [1,0;0,0]; % projection
% vector = [1,2]'; A = vector*vector'; % projection on a vector

%% animation with DOTS
n_steps = 100;
[x,y] = ndgrid(-1:0.15:1);
xy_min = min(min(A*[x(:), y(:)]'))*1.5;
xy_max = max(max(A*[x(:), y(:)]'))*1.5;

dot_colors = jet(length(x(:)));

figure(2);
set(gcf,'color','w')
set(gcf,'position',[250 120 960 540])
scatter(x(:), y(:),30,dot_colors,'filled')
grid on;
xlim([xy_min, xy_max]); ylim([xy_min, xy_max]);
% pause;

step_mtx = (A-eye(2))/n_steps;

for i_steps = 1:n_steps

    new_xy = (eye(2)+step_mtx*i_steps)*[x(:), y(:)]';
    scatter(new_xy(1,:), new_xy(2,:),30,dot_colors,'filled')

    grid on;
    xlim([xy_min, xy_max]); ylim([xy_min, xy_max]);
    axis square

    pause(0.01);
end

%%
v = VideoWriter('permutation_transform.mp4','MPEG-4');
v.Quality = 100;

open(v);
writeVideo(v,F);
close(v);

author: angeloyeo

repo: gongdols

save_path: github-repos/MATLAB/angeloyeo-gongdols


text:
function [s, cfg] = ft_statfun_roc(cfg, dat, design)

% FT_STATFUN_ROC computes the area under the curve (AUC) of the Receiver Operator
% Characteristic (ROC). This is a measure of the separability of the data observed in
% two conditions. The AUC can be used for statistical testing whether the two
% conditions can be distinguished on the basis of the data.
%
% Use this function by calling one of the high-level statistics functions as
%   [stat] = ft_timelockstatistics(cfg, timelock1, timelock2, ...)
%   [stat] = ft_freqstatistics(cfg, freq1, freq2, ...)
%   [stat] = ft_sourcestatistics(cfg, source1, source2, ...)
% with the following configuration option
%   cfg.statistic = 'ft_statfun_roc'
%
% The experimental design is specified as:
%   cfg.ivar  = independent variable, row number of the design that contains the labels of the conditions to be compared (default=1)
%
% The labels for the independent variable should be specified as the number 1 and 2.
%
% Note that this statfun performs a one sided test in which condition "1" is assumed
% to be larger than condition "2". This function does not compute an analytic
% probability of condition "1" being larger than condition "2", but can be used in a
% randomization test, including clustering.
%
% A low-level example with 10 channel-time-frequency points and 1000 observations per
% condition goes like this:
%   dat1 = randn(10,1000) + 1;
%   dat2 = randn(10,1000);
%   design = [1*ones(1,1000) 2*ones(1,1000)];
%   stat = ft_statfun_roc([], [dat1 dat2], design);
%
% See also FT_TIMELOCKSTATISTICS, FT_FREQSTATISTICS or FT_SOURCESTATISTICS

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

% set the defaults
cfg.ivar = ft_getopt(cfg, 'ivar', 1);

% on-the-fly log transformation is not supported any more, please use FT_MATH instead
cfg = ft_checkconfig(cfg, 'forbidden', 'logtransform');
cfg = ft_checkconfig(cfg, 'forbidden', 'numbins');

% start with a quick test to see whether there appear to be NaNs
if any(isnan(dat(1,:)))
  % exclude trials that contain NaNs for all observed data points
  sel    = all(isnan(dat),1);
  dat    = dat(:,~sel);
  design = design(:,~sel);
end

% logical indexing is faster than using find(...)
selA = (design(cfg.ivar,:)==1);
selB = (design(cfg.ivar,:)==2);
% select the data in the two classes
datA = dat(:, selA);
datB = dat(:, selB);

nobs = size(dat,1);
na   = size(datA,2);
nb   = size(datB,2);
auc  = zeros(nobs, 1);

for k = 1:nobs
  % compute the area under the curve for each channel/time/frequency
  a = datA(k,:);
  b = datB(k,:);

  % to speed up the AUC, the critical value is determined by the actual
  % values in class B, which also ensures a regular sampling of the False Alarms
  b = sort(b);

  ca = zeros(nb+1,1);
  ib = zeros(nb+1,1);
  % cb = zeros(nb+1,1);
  % ia = zeros(nb+1,1);

  for i=1:nb
    % for the first approach below, the critval could also be choosen based on e.g. linspace(min,max,n)
    critval = b(i);

    % for each of the two distributions, determine the number of correct and incorrect assignments given the critical value
    % ca(i) = sum(a>=critval);
    % ib(i) = sum(b>=critval);
    % cb(i) = sum(b<critval);
    % ia(i) = sum(a<critval);

    % this is a much faster approach, which works due to using the sorted values in b as the critical values
    ca(i) = sum(a>=critval);   % correct assignments to class A
    ib(i) = nb-i+1;            % incorrect assignments to class B
  end

  % add the end point
  ca(end) = 0;
  ib(end) = 0;
  % cb(end) = nb;
  % ia(end) = na;

  hits = ca/na;
  fa   = ib/nb;

  % the numerical integration is faster if the points are sorted
  hits = fliplr(hits);
  fa   = fliplr(fa);

  if false
    % this part is optional and should only be used when exploring the data
    figure
    plot(fa, hits, '.-')
    xlabel('false positive');
    ylabel('true positive');
    title('ROC-curve');
  end

  % compute the area under the curve using numerical integration
  auc(k) = numint(fa, hits);

end

% return the area under the curve as the statistic of interest
s.stat = auc;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% NUMINT computes a numerical integral of a set of sampled points using
% linear interpolation. Alugh the algorithm works for irregularly sampled points
% along the x-axis, it will perform best for regularly sampled points
%
% Use as
%   z = numint(x, y)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function z = numint(x, y)
if ~all(diff(x)>=0)
  % ensure that the points are sorted along the x-axis
  [x, i] = sort(x);
  y = y(i);
end
n = length(x);
z = 0;
for i=1:(n-1)
  x0 = x(i);
  y0 = y(i);
  dx = x(i+1)-x(i);
  dy = y(i+1)-y(i);
  z = z + (y0 * dx) + (dy*dx/2);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function varargout=mvtest(varargin)
%EIG (overloaded)

switch class(varargin{1})

    case 'double'
        varargout{1} = sum(varargin{1}.^2,2);

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

        x = varargin{1};
        y = yalmip('definemulti',mfilename,x,[size(x,1) 1]);
        varargout{1} = y;

    case 'char'

        X = varargin{3};

        operator = struct('convexity','convex','monotonicity','none','definiteness','none','model','callback');

        varargout{1} = [];
        varargout{2} = operator;
        varargout{3} = X;

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

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [result]=randSimpleMatrix(size)
    result=rand(size,'double', 'gpuArray');
end

author: jiweil

repo: Hierarchical-Neural-Autoencoder


text:
%% tesSmooth_HC
% Below is a demonstration of the features of the |tesSmooth_HC| function

%%
clear; close all; clc;

%% Syntax
% |[P]=tesSmooth_HC(TES,V,IND_V,cPar);|

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function v = cvxcnst( rhs )
v = class( struct( 'rhs_', rhs ), 'cvxcnst' );

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function r8_roundb_test ( )

%*****************************************************************************80
%
%% R8_ROUNDB_TEST tests R8_ROUNDB.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2009
%
%  Author:
%
%    John Burkardt
%
  base = 3;
  x = pi;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8_ROUNDB_TEST\n' );
  fprintf ( 1, '  R8_ROUNDB rounds a number to a\n' );
  fprintf ( 1, '  specified number of base BASE digits.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Here, we will use BASE = %d\n', base );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Test effect on PI:\n' );
  fprintf ( 1, '  X = %f\n', x );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  NPLACE  XROUND\n' );
  fprintf ( 1, '\n' );

  for i = 0 : 20
    nplace = i;
    xround = r8_roundb ( base, nplace, x );
    fprintf ( 1, '  %8d  %f\n', i, xround );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Try with a negative base:\n' );
  x = 121.0;
  base = -3;
  nplace = 3;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Input quantity is X = %f\n', x );
  fprintf ( 1, '  to be rounded in base %d\n', base );

  for nplace = 1 : 5

    xround = r8_roundb ( base, nplace, x );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Output value to %d places is %f\n', nplace, xround );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ files ] = s1manifest_files( domnode )
%S1MANIFEST_FILES Extract relevant measurement and metadata files from a
%Sentinel manifest.safe file and group them together appropriately.
%
% Written by: Wade Schwartzkopf, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

%% Setup
% Overall wrapping node is in a namespace that makes things messier
% (than they already are) in a XPath call, so we go down one level.
domnode = domnode.getFirstChild;
xp=javax.xml.xpath.XPathFactory.newInstance.newXPath();

%% Find associated measurement and metadata files in this package
% First find all of the "Measurement Data Units".  This should provide each
% data object (measurements), together with its metadata and noise and
% calibration files.
mdu_xp_str = ['informationPackageMap/' ...
    '*[local-name()=''contentUnit'']/' ...
    '*[local-name()=''contentUnit''][@repID="s1Level1MeasurementSchema"]'];
num_measurment_data_units = str2double(xp.evaluate(...
    ['count(' mdu_xp_str ')'],...
    domnode));
files(num_measurment_data_units) = struct('data','','product','','noise','','calibration','');
for i = 1:num_measurment_data_units
    % Find data ("measurement") file itself
    data_object_ID = char(xp.evaluate(...
        [mdu_xp_str '[' num2str(i) ']/dataObjectPointer/@dataObjectID'],domnode));
    files(i).data = char(xp.evaluate([...
        'dataObjectSection/' ...
        'dataObject[@ID="' data_object_ID '"]/' ...
        'byteStream/' ...
        'fileLocation/' ...
        '@href'],domnode));
    % Find references to associated metadata files
    associated_DMDs = strsplit(strtrim(char(xp.evaluate(...
        [mdu_xp_str '[' num2str(i) ']/@dmdID'],domnode))));
    % The references are indirect.  They are equivalently pointers to
    % pointers. Not clear why it was done this way, but here we get the
    % real IDs for those files (used in the dataObjectSection).
    associated_IDs=cell(numel(associated_DMDs),1);
    for j = 1:numel(associated_DMDs)
        associated_IDs{j} = char(xp.evaluate([
            'metadataSection/' ...
            'metadataObject[@ID="' associated_DMDs{j} '"]/' ...
            'dataObjectPointer/' ...
            '@dataObjectID'],domnode));
    end
    % Construct a string that will look for any of these IDs
    ID_test_string = ['[@ID="' associated_IDs{1} '"'];
    for j = 2:numel(associated_IDs)
        ID_test_string = [ID_test_string ' or @ID="' associated_IDs{j} '"'];
    end
    ID_test_string = [ID_test_string ']'];
    files(i).product = char(xp.evaluate([...
        'dataObjectSection/' ...
        'dataObject[@repID="s1Level1ProductSchema"]' ... % Must be product type
        ID_test_string '/' ... % and an ID that matches one associated with the current measurement
        'byteStream/' ...
        'fileLocation/' ...
        '@href'],domnode));
    files(i).noise = char(xp.evaluate([...
        'dataObjectSection/' ...
        'dataObject[@repID="s1Level1NoiseSchema"]' ... % Must be noise type
        ID_test_string '/' ... % and an ID that matches one associated with the current measurement
        'byteStream/' ...
        'fileLocation/' ...
        '@href'],domnode));
    files(i).calibration = char(xp.evaluate([...
        'dataObjectSection/' ...
        'dataObject[@repID="s1Level1CalibrationSchema"]' ... % Must be calibration type
        ID_test_string '/' ... % and an ID that matches one associated with the current measurement
        'byteStream/' ...
        'fileLocation/' ...
        '@href'],domnode));
end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
function model = assignSubSystems(model,gpraFile)
% Assigns each reaction a subsystem in the model structure
%
% USAGE:
%
%    model = assignSubSystems(model, gpraFile)
%
% INPUTS:
%    model:       COBRA model structure
%    gpraFile:    SimPheny GPRA file
%
% OUTPUT:
%    model:       COBRA model structure with subsystem assignment
%
% .. Author: - Markus Herrgard

gpraModel = parseSimPhenyGPRA(gpraFile);

[isInList,rxnInd] = ismember(model.rxns,gpraModel.rxns);

for i = 1:length(model.rxns)
    if (isInList(i))
        model.subSystems{i} = gpraModel.subSystems{rxnInd(i)};
    else
        model.subSystems{i} = 'Exchange';
    end
end

model.subSystems = model.subSystems';

author: opencobra

repo: cobratoolbox


text:
function obj = spm2canlab_dataset(obj, subject, spm)
% Extract Event_Level data from subjects' SPM.mat files to add data to
% canlab_dataset object.
%
% :Usage:
% ::
%
%    obj = spm2canlab_dataset(obj, subject, spm)
%
%
% ..
%     Author and copyright information:
%
%     Copyright (C) 2015  Wani Woo
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
% ..
%
%
% :Inputs:
%
%   **obj:**
%        Canlab_dataset object (see canlab_dataset)
%
%   **subject:**
%        Subject list (it could be one subject [in a string
%        format], or it could be multiple subjects in cell array)
%
%   **spm:**
%        This could be loaded SPM (struct), or one path for one
%        subject's SPM.mat file (string), or multiple loaded SPM or
%        paths in cell array
%
% :Outputs:
%
%   **obj:**
%        Canlab_dataset object with new data
%
%
% :Examples:
% ::
%
%    subj = {'dpsp002','dpsp003'};
%    spm = {'dpsp002_SPM.mat', 'dpsp003_SPM.mat'};
%
%    D = canlab_dataset; % if D doesn't exist yet
%    D = spm2canlab_dataset(D, subj, spm);
%
%
% :See also:
%   canlab_dataset
%   spm_mat2batchinput
%

% ..
%    Programmers' notes:
%       After running this, please check D.Subj_Level.data and
%       D.Event_level.data to see if there are NaNs. NaNs could be there when
%       information cannot be extracted from the given SPM.mat files.
% ..

subj = parseSUBJ(subject); % Check 1: parsing subjects

% Check 2: parsing SPM input
SPM = parseSPM(spm);

% Adding data to the canlab_dataset obj
for subj_i = 1:numel(subj)

    % check if the subject exists in obj
    [isexist_sub, wh_subj] = exist_data(obj.Subj_Level.id, subj{subj_i});

    % if the subject is new, add the subject to obj.Subj_Level
    if ~isexist_sub

        obj.Subj_Level.id{wh_subj} = subj{subj_i};

        if ~isempty(obj.Subj_Level.names)

            str = ['D.Subj_Level.data for subject ' subj{subj_i} ' are now NaNs. Please enter the data.'];
            warning(str);

            obj.Subj_Level.data(wh_subj,:) = NaN(1,size(obj.Subj_Level.data,2));

        end
    end

    % get information from SPM
    out = spm2canlab_dataset_format(SPM{subj_i});
    names = out.names;

    for fields_i = 1:numel(names)

        % check the same information exists already
        [isexist_data, wh_col] = exist_data(obj.Event_Level.names, names{fields_i});

        % if the info is new, the existing subjects data for that column should be NaNs.
        if ~isexist_data && numel(obj.Subj_Level.id) ~= 1
            for i = 1:numel(obj.Event_Level.data)
                obj.Event_Level.data{i}(:,wh_col) = NaN(size(obj.Event_Level.data{i},1),1);
            end
            warning('For existing subjects, please check NaNs in D.Event_Level.data.');
        end

        % add the data
        obj.Event_Level.names{wh_col} = names{fields_i};
        obj.Event_Level.type{wh_col} = out.type{fields_i};
        obj.Event_Level.descrip{wh_col} = out.descrip{fields_i};
        obj.Event_Level.data{wh_subj}(:,wh_col) = out.data(:,fields_i);

        wh_cols(fields_i) = wh_col;

    end

    % if the subject is new and there were existing data, the data that the
    % new subject doesn't have should be NaNs
    if ~isexist_sub
        wh_old = ~ismember(1:size(obj.Event_Level.data{wh_subj},2),wh_cols);
        if any(wh_old)
            obj.Event_Level.data{wh_subj}(:,wh_old) = ...
                NaN(size(obj.Event_Level.data{wh_subj}(:,wh_old)));
            warning('For new subjects, please check NaNs in D.Event_Level.data.');
        end
    end

end

end

% ============ SUB-FUNCTIONS ============

function subj = parseSUBJ(subject)

if ischar(subject)
    subj{1} = subject;
else
    subj = subject;
end

end

function SPM = parseSPM(spm)

if ischar(spm)
    temp = load(spm);
    SPM{1} = temp.SPM;
else
    if iscell(spm)
        for i = 1:numel(spm)
            if ischar(spm{i})
                temp = load(spm{i});
                SPM{i} = temp.SPM;
            elseif isstruct(spm{i})
                SPM{i} = spm{i};
            else
                error('SPM is unknown format.');
            end
        end
    elseif isstruct(spm)
        SPM{1} = spm;
    else
        error('SPM is unknown format.');
    end
end

end

function [isexist, wh] = exist_data(obj_field, field)

isexist = any(strcmp(obj_field, field));

if isexist
    wh = find(strcmp(obj_field, field));
else
    wh = numel(obj_field) + 1;
end

end

% ---- MAJOR SUBFUNCTION ----

function out = spm2canlab_dataset_format(SPM)

% Extract Event_Level data from one subject's SPM.mat file to enter data to
% canlab_dataset object. This function is called by "spm2canlab_dataset".
%
% Usage:
% -------------------------------------------------------------------------
% out = spm2canlab_dataset_format(SPM)
%
% Inputs:
% -------------------------------------------------------------------------
% SPM            SPM saved in SPM.mat of the model directory.
%
% Outputs:
% -------------------------------------------------------------------------
% out.names      Event_Level names in cell array
%                (RunNumber, EventName, Onsets, Durations)
% out.descrip    Event_Level descript in cell array
% out.type       Event_Level variable types of data
% out.data       Event_Level data
%
% See also spm2canlab_dataset, canlab_dataset, spm_mat2batchinput
%
% Examples:
% -------------------------------------------------------------------------
% load SPM.mat;
% out = spm2dataset(SPM);
%
% % output:
% %   names: {'RunNumber'  'EventNames'  'Onsets'  'Durations'}
% % descrip: {'Run 1-8'  {12x1 cell}  'Onsets'  'Duration'}
% %    type: {'Categorical'  'Categorical'  'Continuous'  'Continuous'}
% %    data: [192x4 double]
%
% -------------------------------------------------------------------------
% Copyright (C) 2015  Wani Woo

% Programmers' notes:

% Parsing SPM to know the total number of Event types

[~, ~, ~, names] = spm_mat2batchinput(SPM);

% 1. Event_Level variable names
out.names = {'RunNumber','EventNames','Onsets','Durations'};

% 2. Event_Level descript
out.descrip{1} = ['Run 1-' num2str(numel(SPM.Sess))];
for i = 1:numel(unique(names))
    out.descrip{2}{i,1} = [num2str(i) ': ' deblank(names{i})];
end
out.descrip{3} = 'Onsets';
out.descrip{4} = 'Duration';

% 3. Event_Level variable type
out.type{1} = 'Categorical';
out.type{2} = 'Categorical';
out.type{3} = 'Continuous';
out.type{4} = 'Continuous';

% 4. Event_Level variable type
out.data = [];
for run_i = 1:numel(SPM.Sess)

    row_i = 0;
    for eventtype_i = 1:numel(SPM.Sess(run_i).U)

        EventName = find(strcmp(unique(names), SPM.Sess(run_i).U(eventtype_i).name));

        for event_i = 1:numel(SPM.Sess(run_i).U(eventtype_i).ons)

            row_i = row_i+1;
            Durations = SPM.Sess(run_i).U(eventtype_i).dur(event_i);
            Onsets = SPM.Sess(run_i).U(eventtype_i).ons(event_i);
            Ddata_run(row_i,:) = [run_i EventName Onsets Durations];

        end
    end

    % sorting events using onset information
    [~, idx] = sort(Ddata_run(:,3));
    out.data = [out.data; Ddata_run(idx,:)];

end

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
% RANSACFITLINE - fits line to 3D array of points using RANSAC
%
% Usage  [L, inliers] = ransacfitline(XYZ, t, feedback)
%
% This function uses the RANSAC algorithm to robustly fit a line
% to a set of 3D data points.
%
% Arguments:
%          XYZ - 3xNpts array of xyz coordinates to fit line to.
%          t   - The distance threshold between data point and the line
%                used to decide whether a point is an inlier or not.
%          feedback - Optional flag 0 or 1 to turn on RANSAC feedback
%                     information.
%
% Returns:.
%           V - Line obtained by a simple fitting on the points that
%               are considered inliers.  The line goes through the
%               calculated mean of the inlier points, and is parallel to
%               the principal eigenvector.  The line is scaled by the
%               square root of the largest eigenvalue.
%               This line is a n*2 matrix.  The first column is the
%               beginning point, the second column is the end point of the
%               line.
%           L - The two points in the data set that were found to
%               define a line having the most number of inliers.
%               The two columns of L defining the two points.
%           inliers - The indices of the points that were considered
%                     inliers to the fitted line.
%
% See also:  RANSAC, FITPLANE, RANSACFITPLANE

% Copyright (c) 2003-2006 Peter Kovesi and Felix Duvallet (CMU)
% School of Computer Science & Software Engineering
% The University of Western Australia
% http://www.csse.uwa.edu.au/
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software.
%
% The Software is provided "as is", without warranty of any kind.

% Aug  2006 - created ransacfitline from ransacfitplane
%             author: Felix Duvallet

function [V, L, inliers] = ransacfit3Dline(XYZ, t, feedback)

    if nargin == 2
	feedback = 0;
    end

    [rows, npts] = size(XYZ);

    if rows ~=3
        error('data is not 3D');
    end

    if npts < 2
        error('too few points to fit line');
    end

    s = 2;  % Minimum No of points needed to fit a line.

    fittingfn = @defineline;
    distfn    = @lineptdist;
    degenfn   = @isdegenerate;

    [L, inliers] = ransac(XYZ, fittingfn, distfn, degenfn, s, t, feedback);

    % Find the line going through the mean, parallel to the major
    % eigenvector
    V = fitline3d(XYZ(:, inliers));

%------------------------------------------------------------------------
% Function to define a line given 2 data points as required by
% RANSAC.

function L = defineline(X);
    L = X;

%------------------------------------------------------------------------
% Function to calculate distances between a line and an array of points.
% The line is defined by a 3x2 matrix, L.  The two columns of L defining
% two points that are the endpoints of the line.
%
% A line can be defined with two points as:
%        lambda*p1 + (1-lambda)*p2
% Then, the distance between the line and another point (p3) is:
%        norm( lambda*p1 + (1-lambda)*p2 - p3 )
% where
%                  (p2-p1).(p2-p3)
%        lambda =  ---------------
%                  (p1-p2).(p1-p2)
%
% lambda can be found by taking the derivative of:
%      (lambda*p1 + (1-lambda)*p2 - p3)*(lambda*p1 + (1-lambda)*p2 - p3)
% with respect to lambda and setting it equal to zero

function [inliers, L] = lineptdist(L, X, t)

    p1 = L(:,1);
    p2 = L(:,2);

    npts = length(X);
    d = zeros(npts, 1);

    for i = 1:npts
        p3 = X(:,i);

        lambda = dot((p2 - p1), (p2-p3)) / dot( (p1-p2), (p1-p2) );

        d(i) = norm(lambda*p1 + (1-lambda)*p2 - p3);
    end

    inliers = find(abs(d) < t);

%------------------------------------------------------------------------
% Function to determine whether a set of 2 points are in a degenerate
% configuration for fitting a line as required by RANSAC.
% In this case two points are degenerate if they are the same point
% or if they are exceedingly close together.

function r = isdegenerate(X)
    %find the norm of the difference of the two points
    % this will be 0 iff the two points are the same (the norm of their
    % difference is zero)
    r = norm(X(:,1) - X(:,2)) < eps;

author: DrGabor

repo: LiDAR

save_path: github-repos/MATLAB/DrGabor-LiDAR


text:
function L = ngaussNoiseLikelihood(noise, mu, varsigma, y)

% NGAUSSNOISELIKELIHOOD Likelihood of the data under the NGAUSS noise model.
% FORMAT
% DESC returns the likelihood of a data set under the  noiseless Gaussian noise model.
% ARG noise : the noise structure for which the likelihood is required.
% ARG mu : input mean locations for the likelihood.
% ARG varSigma : input variance locations for the likelihood.
% ARG y : target locations for the likelihood.
%
% SEEALSO : ngaussNoiseParamInit, ngaussNoiseLogLikelihood, noiseLikelihood
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005

% NOISE

L = gaussianNoiseLikelihood(noise, mu, varsigma, y);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function setPreferences()
%SETPREFERENCES Sets MATLAB preferences
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Author: Derek Macklin, macklin@stanford.edu
% Affilitation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 9/6/2011

%prevent MATLAB from using OpenGL to avoid segmentation faults
set(0, 'DefaultFigureRenderer', 'painters')
opengl('neverselect')

%set character encoding to UTF-8
if ~isdeployed
    slCharacterEncoding('UTF-8');
end

%default font
set(0, 'defaultAxesFontName', 'Arial')
set(0, 'defaultTextFontName', 'Arial')

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function [Fval,A,f,sig,sd] = ftestc(data,params,p,plt)
% computes the F-statistic for sine wave in locally-white noise (continuous data).
%
% [Fval,A,f,sig,sd] = ftestc(data,params,p,plt)
%
%  Inputs:
%       data        (data in [N,C] i.e. time x channels/trials or a single
%       vector) - required.
%       params      structure containing parameters - params has the
%       following fields: tapers, Fs, fpass, pad
%           tapers : precalculated tapers from dpss or in the one of the following
%                    forms:
%                    (1) A numeric vector [TW K] where TW is the
%                        time-bandwidth product and K is the number of
%                        tapers to be used (less than or equal to
%                        2TW-1).
%                    (2) A numeric vector [W T p] where W is the
%                        bandwidth, T is the duration of the data and p
%                        is an integer such that 2TW-p tapers are used. In
%                        this form there is no default i.e. to specify
%                        the bandwidth, you have to specify T and p as
%                        well. Note that the units of W and T have to be
%                        consistent: if W is in Hz, T must be in seconds
%                        and vice versa. Note that these units must also
%                        be consistent with the units of params.Fs: W can
%                        be in Hz if and only if params.Fs is in Hz.
%                        The default is to use form 1 with TW=3 and K=5
%
%	        Fs 	        (sampling frequency) -- optional. Defaults to 1.
%           fpass       (frequency band to be used in the calculation in the form
%                                   [fmin fmax])- optional.
%                                   Default all frequencies between 0 and Fs/2
%	        pad		    (padding factor for the FFT) - optional (can take values -1,0,1,2...).
%                    -1 corresponds to no padding, 0 corresponds to padding
%                    to the next highest power of 2 etc.
%			      	 e.g. For N = 500, if PAD = -1, we do not pad; if PAD = 0, we pad the FFT
%			      	 to 512 points, if pad=1, we pad to 1024 points etc.
%			      	 Defaults to 0.
%	    p		    (P-value to calculate error bars for) - optional.
%                           Defaults to 0.05/N where N is the number of samples which
%	                 corresponds to a false detect probability of approximately 0.05.
%       plt         (y/n for plot and no plot respectively)
%
%  Outputs:
%       Fval        (F-statistic in frequency x channels/trials form)
%  	    A		    (Line amplitude for X in frequency x channels/trials form)
%	    f		    (frequencies of evaluation)
%       sig         (F distribution (1-p)% confidence level)
%       sd          (standard deviation of the amplitude C)
if nargin < 1; error('Need data'); end;
if nargin < 2 || isempty(params); params=[]; end;
[tapers,pad,Fs,fpass,err,trialave,params]=getparams(params);
clear err trialave
data=change_row_to_column(data);
[N,C]=size(data);
if nargin<3 || isempty(p);p=0.05/N;end;
if nargin<4 || isempty(plt); plt='n';end;
tapers=dpsschk(tapers,N,Fs); % calculate the tapers
[N,K]=size(tapers);
nfft=max(2^(nextpow2(N)+pad),N);% number of points in fft
[f,findx]=getfgrid(Fs,nfft,fpass);% frequency grid to be returned
% errorchk = 0; % set error checking to default (no errors calculated)
% if nargout <= 3 % if called with 4 output arguments, activate error checking
%     errorchk = 0;
% else
%     errorchk = 1;
% end
Kodd=1:2:K;
Keven=2:2:K;
J=mtfftc(data,tapers,nfft,Fs);% tapered fft of data - f x K x C
Jp=J(findx,Kodd,:); % drop the even ffts and restrict fft to specified frequency grid - f x K x C
tapers=tapers(:,:,ones(1,C)); % add channel indices to the tapers - t x K x C
H0 = squeeze(sum(tapers(:,Kodd,:),1)); % calculate sum of tapers for even prolates - K x C
if C==1;H0=H0';end;
Nf=length(findx);% number of frequencies
H0 = H0(:,:,ones(1,Nf)); % add frequency indices to H0 - K x C x f
H0=permute(H0,[3 1 2]); % permute H0 to get dimensions to match those of Jp - f x K x C
H0sq=sum(H0.*H0,2);% sum of squares of H0^2 across taper indices - f x C
JpH0=sum(Jp.*squeeze(H0),2);% sum of the product of Jp and H0 across taper indices - f x C
A=squeeze(JpH0./H0sq); % amplitudes for all frequencies and channels
Kp=size(Jp,2); % number of even prolates
Ap=A(:,:,ones(1,Kp)); % add the taper index to C
Ap=permute(Ap,[1 3 2]); % permute indices to match those of H0
Jhat=Ap.*H0; % fitted value for the fft

num=(K-1).*(abs(A).^2).*squeeze(H0sq);%numerator for F-statistic
den=squeeze(sum(abs(Jp-Jhat).^2,2)+sum(abs(J(findx,Keven,:)).^2,2));% denominator for F-statistic
Fval=num./den; % F-statisitic
if nargout > 3
   sig=finv(1-p,2,2*K-2); % F-distribution based 1-p% point
   var=den./(K*squeeze(H0sq)); % variance of amplitude
   sd=sqrt(var);% standard deviation of amplitude
end;
if nargout==0 || strcmp(plt,'y');
   [S,f]=mtspectrumc(detrend(data),params);subplot(211); plot(f,10*log10(S));xlabel('frequency Hz'); ylabel('Spectrum dB');
   subplot(212);plot(f,Fval); line(get(gca,'xlim'),[sig sig],'Color','r');xlabel('frequency Hz');
   ylabel('F ratio');
end
A=A*Fs;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
% Test whether it is necessary to rotate VB PCA
function nc2010_publish_datasets

n = 200;
m = 50;
d = 10;

ncomps = 30;

randn('state', 1);
rand('state', 1);

% Generate multivariate normal data
disp('Using weak subspace')
eigs = (1 + [d:-1:1 zeros(1,m-d)]) .^ 2
datastring = 'weak';
plot_eigenvalues(eigs, datastring);

disp('Using strong subspace')
eigs = ([5*ones(1,d) 1*ones(1,m-d)]) .^ 2
datastring = 'strong';
plot_eigenvalues(eigs, datastring);

disp('Using no subspace')
eigs = (m:-1:1) .^ 2
datastring = 'no';
plot_eigenvalues(eigs, datastring);

function plot_eigenvalues(eigs, datastring)

filename = sprintf(['/home/jluttine/papers/neurocomputing2010/' ...
                    'fig_eigenvalues_dataset=%s'], datastring);

figure
plot(sqrt(eigs), 'k-');
set(gcf, 'units', 'centimeters', 'paperunits', 'centimeters');
pos = get(gcf, 'position');
set(gcf, 'position', [pos(1:2), 6,4]);
pos = get(gcf, 'paperposition');
set(gcf, 'paperposition', [pos(1:2),6,4])

ylim([0, max(sqrt(eigs))+1]);
xlim([1, length(eigs)]);

%ylabel('standard deviation')

print(gcf, '-depsc2', filename);

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function [isPhysioVisibleForSpmBatchEditor, pathSpm, pathPhysIO] = ...
    tapas_physio_check_spm_batch_editor_integration(isVerbose)
% Checks whether PhysIO-configuration file for matlabbatch is in subfolder
% of SPM/toolbox and returns a warning, if not
%
%   [isPhysioVisibleForSpmBatchEditor, pathSpm, pathPhysIO] = ...
%       tapas_physio_check_spm_batch_editor_integration()
%
% IN
%   isVerbose   if true, warning is issued
% OUT
%
% EXAMPLE
%   tapas_physio_check_spm_batch_editor_integration
%
%   See also

% Author: Lars Kasper
% Created: 2018-02-17
% Copyright (C) 2018 TNU, Institute for Biomedical Engineering,
%                    University of Zurich and ETH Zurich.
%
% This file is part of the TAPAS PhysIO Toolbox, which is released under the terms of the GNU General Public
% License (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.
%

if nargin < 1
    isVerbose = true;
end

[isSpmOnPath, pathSpm] = tapas_physio_check_spm();

isPhysioVisibleForSpmBatchEditor = isSpmOnPath; % minimum requirement for integration: SPM works!

% check for config matlabbatch file
if isSpmOnPath
    filePhysioCfgMatlabbatch = ...
        dir(fullfile(pathSpm, 'toolbox', '**/tapas_physio_cfg_matlabbatch.m'));

    isPhysioVisibleForSpmBatchEditor = ~isempty(filePhysioCfgMatlabbatch);

    % also important to set default modality of spm to fMRI and
    % initialize batch editor, if not done before
    spm('defaults', 'fmri');
    if ~exist('cfg_util')
        spm_jobman('initcfg');
    end
end

[isPhysioOnPath, pathPhysIO] = tapas_physio_check_path();

if ~isPhysioVisibleForSpmBatchEditor && isVerbose
    warning(['\n The PhysIO Toolbox code folder has not been copied (or linked)' ...
        ' to a subfolder of the SPM/toolbox folder. \n The Batch Editor will' ...
        ' not show PhysIO. \n To make PhysIO visible there, link its path' ...
        ' by typing the following in the Matlab command window (Linux/Mac) \n ' ...
        ' and restart SPM afterwards:\n' ...
        '       ln -s %s %s/toolbox/PhysIO'], pathPhysIO, pathSpm);
end

author: translationalneuromodeling

repo: tapas


text:
function runtop
% See Article 8.6
a=.01; h=4; b=2; n=6; m=3; k=20;
psi=45; theta=-45; phi=0;
for j=1:2:720

   psi=180+j; theta=-45+5*sin(.15*psi); phi=1.5*j;
   topdraw(a,h,b,8,m,18,psi,theta,phi);

   % topdraw(a,h,b,8,m,18,j,-45,1.5*j);
end
pause(1),%  close

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

function [x,y, z]=topdraw(a,h,b,n,m,k,psi,theta,phi)
% [x,y, z]=topdraw(a,h,b,n,m,k,psi,theta,phi)

% [X,Y,Z]=topsurf(n,m,k,a,h,b)
if nargin==0
  a=.01; h=4; b=2; n=6; m=3; k=20;
  psi=45; theta=-45; phi=0;
end
%s=sqrt(h^2+(a-b)^2); uax=cumsum([0,a,s,b]);
%u=cornrpts(uax/uax(4),nax+1);
%v=linspace(0,1,ncrc+1);

[dx,dy,dz]=topsurf(n,m,k,a,h,b);

x=dz; y=dx; z=dy;

m=eulerang(psi,theta,phi);
p=size(x); v=[x(:),y(:),z(:)]*m;
x=reshape(v(:,1),p); y=reshape(v(:,2),p);
z=reshape(v(:,3),p);
w=[-1 1 -1 1 -1 1]*sqrt(h^2+b^2);

%rotate3d on;
clf; surf(x,y,z),  view([-45,30])
axis equal; axis(w); axis on
xlabel('x axis'), ylabel('y axis')
zlabel('z axis')
title('NUTATING TOP PRECESSION')
colormap([127/255 1 212/255]); drawnow, shg

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

function m=eulerang(psi,theta,phi)
% m=eulerang(psi,theta,phi)
a=pi/180*[psi,theta,phi]; c=cos(a); s=sin(a);
m=[c(1)*c(2), s(1)*c(2), -s(2); -s(1)*c(3)+...
   c(1)*s(2)*s(3), c(1)*c(3)+s(1)*s(2)*s(3),...
   c(2)*s(3); s(1)*s(3)+c(1)*s(2)*c(3),...
  -c(1)*s(3)+s(1)*s(2)*c(3), c(2)*c(3)];

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

function [X,Y,Z]=topsurf(n,m,k,a,h,b)
% [X,Y,Z]=topsurf(n,m,k,a,h,b)
if nargin==0
n=8; m=4; k=20; a=.2; h=4; b=1;
end
tol=100*eps*(a+h+b);
a=a+(a==0)*tol; b=b+(b==0)*tol;

D=2*pi/n; u=cos(D/2); v=sin(D/2);
z=u+i*linspace(-v,v,m+1)'; z(m+1)=[];
z=z*exp(i*D*(0:n-1)); z=[z(:);u-i*v];
x=real(z); y=imag(z); N=length(x);
% plot(x,y,x,y,'.'), axis equal, shg, pause

qd=cumsum([0;a;sqrt(h^2+(a-b)^2);b]);
qd=qd/max(qd); q=cornrpts(qd,k);
K=length(q);
Z=interp1(qd,[0;0;h;h],q(:))*ones(1,N);
r=interp1(qd,[0;a;b;0],q); r=r(:);
X=r*x(:)'; Y=r*y(:)';

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

function v=cornrpts(u,N)
% v=cornrpts(u,N)
% This function generates a set of approximately
% N points between min(u) and max(u) including
% all points in u plus additional points evenly
% spaced in each successive interval.
% u   -  vector of points
% N   -  approximate number of output points
%        between min(u(:)) and max(u(:))
% v   -  vector of points in increasing order
u=sort(u(:))'; np=length(u); d=u(np)-u(1); v=u(1);
for j=1:np-1
  dj=u(j+1)-u(j); nj=max(1,fix(N*dj/d));
  v=[v,[u(j)+dj/nj*(1:nj)]];
end

author: Sable

repo: mcbench-benchmarks


text:
function varargout=bark2frq(varargin)
%
% For calling details please see v_bark2frq.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_bark2frq.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: bark2frq.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_bark2frq(varargin{:});
else
    v_bark2frq(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function value = r8_atan ( y, x )

%*****************************************************************************80
%
%% R8_ATAN computes the inverse tangent of the ratio Y / X.
%
%  Discussion:
%
%    R8_ATAN returns an angle whose tangent is ( Y / X ), a job which
%    the built in functions ATAN and ATAN2 already do.
%
%    However:
%
%    * R8_ATAN always returns a positive angle, between 0 and 2 PI,
%      while ATAN and ATAN2 return angles in the interval [-PI/2,+PI/2]
%      and [-PI,+PI] respectively;
%
%    * R8_ATAN accounts for the signs of X and Y, (as does ATAN2).  The ATAN
%     function by contrast always returns an angle in the first or fourth
%     quadrants.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real Y, X, two quantities which represent the tangent of
%    an angle.  If Y is not zero, then the tangent is (Y/X).
%
%    Output, real VALUE, an angle between 0 and 2 * PI, whose tangent is
%    (Y/X), and which lies in the appropriate quadrant so that the signs
%    of its cosine and sine match those of X and Y.
%

%
%  Special cases:
%
  if ( x == 0.0 )

    if ( 0.0 < y )
      value = pi / 2.0;
    elseif ( y < 0.0 )
      value = 3.0 * pi / 2.0;
    elseif ( y == 0.0 )
      value = 0.0;
    end

  elseif ( y == 0.0 )

    if ( 0.0 < x )
      value = 0.0;
    elseif ( x < 0.0 )
      value = pi;
    end
%
%  We assume that ATAN2 is correct when both arguments are positive.
%
  else

    abs_y = abs ( y );
    abs_x = abs ( x );

    theta_0 = atan2 ( abs_y, abs_x );

    if ( 0.0 < x && 0.0 < y )
      value = theta_0;
    elseif ( x < 0.0 && 0.0 < y )
      value = pi - theta_0;
    elseif ( x < 0.0 && y < 0.0 )
      value = pi + theta_0;
    elseif ( 0.0 < x && y < 0.0 )
      value = 2.0 * pi - theta_0;
    end

  end

  return
end


text:
% DEMO_BBCIONLINE_APPLY_LSL_STREAMING
% this is just for testing / demonstrating the streaming of EEG data using
% lsl

% set up dummy classifier
C= struct('b', 0);
C.w= randn(8*2, 1);  % 2 log-bandpower feature per channel

% setup the bbci variable to define the online processing chain
bbci= struct;
bbci.source.acquire_fcn= @bbci_acquire_lsl;

% define the dummy electrode setting
clab = {'AF5' 'AF3' 'AF1' 'AFz' 'AF2' 'AF4' 'AF6' ...
        'F5' 'F3', 'F1' 'Fz' 'F2' 'F4' 'F6' ...
        'FC7' 'FC5'};

bbci.signal.clab = clab;
% provide clab and markerstreamname to lsl acquire function
bbci.source.acquire_param = {'clab', clab, 'markerstreamname', 'MyMarkerStream'};

bbci.feature.proc= {@proc_variance, @proc_logarithm};
bbci.feature.ival= [-500 0];

bbci.classifier.C= C;

% specify file locations and logging
% bbci.log.output= 'screen&file';
% bbci.log.file= fullfile(BTB.DataDir, 'tmp\log');
bbci.source.record_signals = 0;
bbci.source.record_basename = fullfile(BTB.DataDir,'tmp\lsl_test');
bbci.quit_condition.running_time = 170;
bbci_apply(bbci);

author: bbci


text:
function [La,dLa,lambda0] = lagrange(U,s,b,more)
%LAGRANGE Plot the Lagrange function for Tikhonov regularization.
%
% [La,dLa,lambda0] = lagrange(U,s,b,more)
% [La,dLa,lambda0] = lagrange(U,sm,b,more)  ,  sm = [sigma,mu]
%
% Plots the Lagrange function
%    La(lambda) = || A x - b ||^2 + lambda^2*|| L x ||^2
% and its first derivative dLa = dLa/dlambda versus lambda.
% Here, x is the Tikhonov regularized solution.
%
% If nargin = 4, || A x - b || and || L x || are also plotted.
%
% Returns La, dLa, and the value lambda0 of lambda for which
% dLa has its minimum.

% Per Christian Hansen, IMM, Feb. 21, 2001.

% Set default number of points.
npoints = 200;

% Initialization.
[m,n] = size(U); [p,ps] = size(s);
beta = U'*b; beta2 = norm(b)^2 - norm(beta)^2;
if (ps==2)
  s = s(p:-1:1,1)./s(p:-1:1,2); beta = beta(p:-1:1);
end
xi = beta(1:p)./s;

% Compute the L-curve.
eta = zeros(npoints,1); rho = eta;
lambda(npoints,1) = s(p);
ratio = (s(1)/s(p))^(1/(npoints-1));
for i=npoints-1:-1:1, lambda(i) = ratio*lambda(i+1); end
for i=1:npoints
  f = fil_fac(s,lambda(i));
  eta(i) = norm(f.*xi);
  rho(i) = norm((1-f).*beta(1:p));
end
if (m > n && beta2 > 0), rho = sqrt(rho.^2 + beta2); end

% Compute the Lagrange function and its derivative.
La = rho.^2 + (lambda.^2).*(eta.^2);
dLa = 2*lambda.*(eta.^2);
[mindLa,mindLi] = min(dLa); lambda0 = lambda(mindLi);

% Plot the functions.
if (nargin==3)
  loglog(lambda,La,'-',lambda,dLa,'--',lambda0,mindLa,'o')
  legend('La','dLa/d\lambda')
else
  loglog(lambda,La,'-',lambda,dLa,'--',lambda,eta,':',lambda,rho,'-.',...
         lambda0,mindLa,'o')
  legend('La','dLa/d\lambda','|| L x ||_2','|| A x - b ||_2')
end
xlabel('\lambda')
title('Lagrange function La and its derivative')

repo: REST


text:
function [num_all_images, sensors, paths, info_count_text] = CheckImagesPath(path_data)
%CHECKIMAGEPATH Ensure the input path is right to find a Landsats 4-8, and
%Sentinel 2 image(s).
% path_data - the input path

    % searching deeps. 0 is default.
    [image_types_paths] = CheckImagePath(path_data,0);

    % all count info foe searched images at current folder.
    num_all_images = size(image_types_paths,1);
    num_L4_tm = 0;
    num_L5_tm = 0;
    num_L6_tm = 0;
    num_L7_tm_plus = 0;
    num_L8_oli_tirs = 0;
    num_S2A_msi = 0;
    num_S2B_msi = 0;

    % key info: the senor and path are needed as outputs, the first one is
    % to determine the default parameters for Fmask, and the second one is
    % to determine the path loading all data.
    sensors = [];
    paths = [];
    info_count_text = [];

    for i_all = 1:num_all_images
        % c: current.
        cimage_sensor = image_types_paths{i_all,1};
        cimage_num = image_types_paths{i_all,2};
        cimage_type = Convert2ImageType(cimage_sensor,cimage_num);
        clear cimage_num;
        cimage_path = image_types_paths{i_all,3};
        switch cimage_type
            case {'Landsat 4 TM'}
                num_L4_tm=num_L4_tm+1;
            case {'Landsat 5 TM'}
                num_L5_tm=num_L5_tm+1;
            case {'Landsat 6 TM'}
                num_L6_tm=num_L6_tm+1;
            case {'Landsat 7 ETM+'}
                num_L7_tm_plus=num_L7_tm_plus+1;
            case {'Landsat 8 OLI/TIRS'}
                num_L8_oli_tirs=num_L8_oli_tirs+1;
            case {'Sentinel 2A MSI'}
                % further check Sentinel 2 image folder, see there is
                % .SAFE.
                if contains(cimage_path,'.SAFE')
                    num_S2A_msi=num_S2A_msi+1;
                else
                    cimage_sensor = [];% no available Sentinel 2 data.
                end
            case {'Sentinel 2B MSI'}
                % see there is .SAFE.
                if contains(cimage_path,'.SAFE')
                    num_S2B_msi=num_S2B_msi+1;
                else
                    cimage_sensor = [];% no available Sentinel 2 data.
                end
        end
        if ~isempty(cimage_sensor)
            sensors = [sensors;{cimage_sensor}];
            paths = [paths;{cimage_path}];
        end
        clear cimage_sensor;
    end
    % renew num_all_images
    num_all_images = length(sensors);
    % used to notice user.
    text_line = 0;
    % multide
    if isequal(num_all_images,num_L4_tm)||...
            isequal(num_all_images,num_L5_tm)||...
            isequal(num_all_images,num_L6_tm)||...
            isequal(num_all_images,num_L7_tm_plus)||...
            isequal(num_all_images,num_L8_oli_tirs)||...
            isequal(num_all_images,num_S2A_msi)||...
            isequal(num_all_images,num_S2B_msi)
       % only for 1 type image
        if num_L4_tm > 0
        text_line = text_line+1;
        info_count_text{text_line} = sprintf('%s Landsat 4 TM images are found at ''%s''\n',...
            num2str(num_L4_tm),path_data);
        end
        if num_L5_tm > 0
            text_line = text_line+1;
            info_count_text{text_line} = sprintf('%s Landsat 5 TM images are found at ''%s''\n',...
                num2str(num_L5_tm),path_data);
        end
        if num_L6_tm > 0
            text_line = text_line+1;
            info_count_text{text_line} = sprintf('%s Landsat 6 TM images are found at ''%s''\n',...
            num2str(num_L6_tm),path_data);
        end
        if num_L7_tm_plus > 0
            text_line = text_line+1;
            info_count_text{text_line} = sprintf('%s Landsat 7 ETM+ images are found at ''%s''\n',...
                num2str(num_L7_tm_plus),path_data);
        end
        if num_L8_oli_tirs > 0
            text_line = text_line+1;
            info_count_text{text_line} = sprintf('%s Landsat 8 OLI/TIRS images are found at ''%s''\n',...
                 num2str(num_L8_oli_tirs),path_data);
        end
        if num_S2A_msi > 0
            text_line = text_line+1;
            info_count_text{text_line} = sprintf('%s Sentinel 2A MSI images are found at ''%s''\n',...
                num2str(num_S2A_msi),path_data);
        end
        if num_S2B_msi > 0
            text_line = text_line+1;
            info_count_text{text_line} = sprintf('%s Sentinel 2B MSI images are found at ''%s''\n',...
                num2str(num_S2B_msi),path_data);
        end
    else
         text_line = text_line+1;
         info_count_text{text_line} = sprintf(...
            'A total of %s images (as follows) are found at ''%s''\n',...
            num2str(num_all_images), path_data);
        if num_L4_tm > 0
        text_line = text_line+1;
        info_count_text{text_line} = sprintf('%s Landsat 4 TM images\n',...
            num2str(num_L4_tm));
        end
        if num_L5_tm > 0
            text_line = text_line+1;
            info_count_text{text_line} = sprintf('%s Landsat 5 TM images\n',...
                num2str(num_L5_tm));
        end
        if num_L6_tm > 0
            text_line = text_line+1;
            info_count_text{text_line} = sprintf('%s Landsat 6 TM images\n',...
            num2str(num_L6_tm));
        end
        if num_L7_tm_plus > 0
            text_line = text_line+1;
            info_count_text{text_line} = sprintf('%s Landsat 7 ETM+ images\n',...
                num2str(num_L7_tm_plus));
        end
        if num_L8_oli_tirs > 0
            text_line = text_line+1;
            info_count_text{text_line} = sprintf('%s Landsat 8 OLI/TIRS images\n',...
                 num2str(num_L8_oli_tirs));
        end
        if num_S2A_msi > 0
            text_line = text_line+1;
            info_count_text{text_line} = sprintf('%s Sentinel 2A MSI images\n',...
                num2str(num_S2A_msi));
        end
        if num_S2B_msi > 0
            text_line = text_line+1;
            info_count_text{text_line} = sprintf('%s Sentinel 2B MSI images\n',...
                num2str(num_S2B_msi));
        end
    end

    % if no data is found, give mention info.
    if isempty(sensors)
        text_line = text_line+1;
        info_count_text{text_line} = sprintf('%s available images are found at ''%s''\n',...
            '0',path_data);
        text_line = text_line+1;
        info_count_text{text_line} = sprintf('Please ensure the path is correct\n');
    end
%     fprintf(info_count);
end

function [image_types_paths] = CheckImagePath(path_data,subfolder_level)
%CHECKIMAGEPATH Ensure the input path is right to find a Landsats 4-8, and
%Sentinel 2 image(s).
% path_data - the input path
% subfolder_level - the level of subfolders that can be used to limit
% searching deeps. 0 is default.

    % If the searching deeps are more than 5, stop and return;
    if subfolder_level > 5
        image_types_paths = [];
        return;
    end
	image_types_paths = [];
    % first search the image(s) at current folder.image_types_paths
    [sensor,num_Lst,~,~] = LoadSensorType(path_data);
    if isempty(sensor)
        % if no available image at current folder,
        % and search its subfolders.
          subfolders = dir(path_data);
          for i_sub=1:length(subfolders)
              % filter out the file names starting with '.', that is not
              % right folder (system crashes).
              if strcmp(subfolders(i_sub).name(1),'.')
                  continue;
              end
              % go to search the images at each subfolder
              path_subfoler = fullfile(subfolders(i_sub).folder,...
                  subfolders(i_sub).name);
              [image_types_paths_sub] = CheckImagePath(path_subfoler,subfolder_level+1);
              if ~isempty(image_types_paths_sub)
                  image_types_paths =[image_types_paths;image_types_paths_sub];
              end
          end

    else
        % successfully searched a supported image.
        % and, return the sensor and the image path
        image_types_paths = [image_types_paths;{sensor,num_Lst,path_data}];
    end
end
function image_type = Convert2ImageType(sensor,num_Lst)
% CONVERT@IMAGETYPE Contruct image type from the sensor and num

    switch sensor
        case 'L_TM'
            image_type = ['Landsat ',num_Lst,' TM'];
        case 'L_ETM_PLUS'
            image_type = ['Landsat ',num_Lst,' ETM+'];
        case 'L_OLI_TIRS'
            image_type = ['Landsat ',num_Lst,' OLI/TIRS'];
        case 'S_MSI'
            image_type = ['Sentinel ',num_Lst,' MSI'];
        otherwise
            image_type=[];
%             error(['Errors occur when searching images at ''', path_data],'''.');
    end
end

author: GERSL

repo: Fmask

save_path: github-repos/MATLAB/GERSL-Fmask


text:
function [tsave, xsave] = quadSim(traj_obj, model_param, KK, t_sim, x0, Fmat, options)

[tsave, xsave] = ode45(@(t,s) quadDynamics(t, s, traj_obj, model_param, KK, Fmat), t_sim, x0, options);

desired_pos = zeros(length(tsave),3);
desired_vel = zeros(length(tsave),3);
desired_acc = zeros(length(tsave),3);

for i = 1:length(tsave)
    desired_s = desiredState(traj_obj, tsave(i));
    desired_pos(i,:) = desired_s.pos';
    desired_vel(i,:) = desired_s.vel';
    desired_acc(i,:) = desired_s.acc';
end

% position,attitude error plot
figure(2)
subplot(4,1,1)
plot(tsave,xsave(:,1)-desired_pos(:,1),'-b','LineWidth',1.0);title('Position error, x-x_d (m)');ylabel('x')
grid on
subplot(4,1,2)
plot(tsave,xsave(:,2)-desired_pos(:,2),'-b','LineWidth',1.0);ylabel('y')
grid on
subplot(4,1,3)
plot(tsave,xsave(:,3)-desired_pos(:,3),'-b','LineWidth',1.0);ylabel('z');%xlabel('time, sec')
grid on
% Orientation error plot (using error function)
yawd = desired_s.yaw;
yawd_dot = desired_s.yawdot;
yawd_2dot = desired_s.yawddot;
xcd = [cos(yawd) sin(yawd) 0]';
zbd = zeros(length(tsave),3);
ybd = zeros(length(tsave),3);
so3_error = zeros(length(tsave),1);
for i = 1:length(tsave)
    des_t = desired_acc(i,:) + [0 0 model_param.grav];
    zbd(i,:) = des_t/norm(des_t);
    ybd(i,:) = (hat_optr(zbd(i,:))*xcd/norm(hat_optr(zbd(i,:))*xcd))';
    xbd = hat_optr(ybd(i,:))*zbd(i,:)';
    Rd = [xbd ybd(i,:)' zbd(i,:)'];
    R = ROTZ(xsave(i,9))*ROTX(xsave(i,7))*ROTY(xsave(i,8));
    so3_error(i) = error_so3(R,Rd);
end
figure(2)
subplot(4,1,4)
plot(tsave,so3_error,'-b','LineWidth',1.0);title('Attitude error function, \Psi');xlabel('time, sec')
grid on
hold off
% Input plot
calc_acc = diff(xsave(:,4:6))./diff(tsave);
calc_ang_acc = diff(xsave(:,10:12))./diff(tsave);
t_vec = calc_acc;
t_vec(:,3) = t_vec(:,3)+model_param.grav;
norm_t = zeros(size(t_vec,1),1);
M = zeros(size(t_vec,1),3);
for i=1:length(norm_t)
    norm_t(i) = norm(t_vec(i,:));
    M(i,:) = model_param.I*calc_ang_acc(i,:)' + hat_optr(xsave(i,10:12))*model_param.I*xsave(i,10:12)';
end
u1 = model_param.mass*norm_t;

% plot for total thrust and moment
% figure(3)
% subplot(4,1,1)
% plot(tsave(1:end-1)+diff(tsave)/2,u1)
% subplot(4,1,2)
% plot(tsave(1:end-1)+diff(tsave)/2,M(:,1))
% subplot(4,1,3)
% plot(tsave(1:end-1)+diff(tsave)/2,M(:,2))
% subplot(4,1,4)
% plot(tsave(1:end-1)+diff(tsave)/2,M(:,3))
U = [u1 M];
L = model_param.arm_length;
c_tf = model_param.c_tf;
mapping_u = [1 1 1 1;0 L 0 -L;-L 0 L 0;c_tf -c_tf c_tf -c_tf];
motor_F = zeros(size(U));
for i=1:length(U)
    trans_F = mapping_u\U(i,:)';
    motor_F(i,:) = trans_F';
end
% plot for thrust of each rotor
figure(3)
subplot(4,1,1)
plot(tsave(1:end-1)+diff(tsave)/2,motor_F(:,1),'-b','LineWidth',1.0);title('Thrust of each rotor, (N)')
grid on
subplot(4,1,2)
plot(tsave(1:end-1)+diff(tsave)/2,motor_F(:,2),'-b','LineWidth',1.0)
grid on
subplot(4,1,3)
plot(tsave(1:end-1)+diff(tsave)/2,motor_F(:,3),'-b','LineWidth',1.0)
grid on
subplot(4,1,4)
plot(tsave(1:end-1)+diff(tsave)/2,motor_F(:,4),'-b','LineWidth',1.0);xlabel('time, sec')
grid on
hold off
end

author: yorgoon

repo: minimum-snap-geometric-control


text:
function [out] = rowvector(in)
% Transposes the input vector to a row vector if it is not
% already a row vector. Called by `verifyModel`.
%
% USAGE:
%
%    out = columnvector(in)
%
% INPUTS:
%    in:     Vector
%
% OUTPUTS:
%    out:    Column vector
%
% Please cite:
% `Sauls, J. T., & Buescher, J. M. (2014). Assimilating genome-scale
% metabolic reconstructions with modelBorgifier. Bioinformatics
% (Oxford, England), 30(7), 1036?8`. http://doi.org/10.1093/bioinformatics/btt747
%
% ..
%    Edit the above text to modify the response to help addMetInfo
%    Last Modified by GUIDE v2.5 06-Dec-2013 14:19:28
%    This file is published under Creative Commons BY-NC-SA.
%
%    Correspondance:
%    johntsauls@gmail.com
%
%    Developed at:
%    BRAIN Aktiengesellschaft
%    Microbial Production Technologies Unit
%    Quantitative Biology and Sequencing Platform
%    Darmstaeter Str. 34-36
%    64673 Zwingenberg, Germany
%    www.brain-biotech.de

dim = size(in) ;
if dim(1) > dim(2)
    out = in' ;
else
    out = in ;
end


text:
function ptch = pix2patch(I,pix,hsize,vsize)
%
% PIX2PATCH  Patch extraction
%   PATCH = PIX2PATCH(I,PIX,HSIZE,VSIZE)  extracts from the image I(MxN) the
%   patch of size HSIZExVSIZE centered at pixel PIX. If VSIZE is not
%   specified, a square patch HSIZExHSIZE is considered.
%
%   PATCH = PIX2PATCH(IMID,PIX,HSIZE,VSIZE)  extracts from the global image
%   Img{IMID} the patch of size HSIZExVSIZE centered at pixel PIX. IMID
%   is an integer. If VSIZE is not specified, a square patch HSIZExHSIZE is
%   considered
%
%   PATCH is a structure containing:
%     I    : matrix with all pixels - the patch itself
%     SI   : sum of all patch pixels;
%     SII  : sum of all squared pixels

global Img

if isscalar(I)  % ndims(I) == 1
    imId = I;
else
    Img{1} = I;
    imId  = 1;
end

if nargin < 4
    vsize = hsize;
end

% minimum necessary margins from central pixel to image edges
h2 = (hsize-1)/2;
v2 = (vsize-1)/2;

pix0 = round(pix); % in case the given pixel is not integer

if inSquare(pix,[1+h2 size(Img{imId},2)-h2 1+v2 size(Img{imId},1)-v2])

    ptch.I    = single(Img{imId}(pix0(2)-v2:pix0(2)+v2,pix0(1)-h2:pix0(1)+h2));
    ptch.SI   = sum(sum(ptch.I));
    ptch.SII  = sum(sum(ptch.I.*ptch.I));

else % pixel is too close to the edge: return blank patch

    ptch.I = single(zeros(vsize,hsize));
    [ptch.SI,ptch.SII] = deal(0);
    ptch.bias = [0;0];

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

%   # START GPL LICENSE

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

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

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
% VL_SAMPLINTHIST  Sample integral histogram
%   HISTS = VL_SAMPLINTHIST(INTHIST, BOXES) samples the integral
%   histogram INTHIST to obtain the histograms of the specified
%   BOXES.
%
%   INTHIST is a MxNxK array, where M x N are ``spatial'' dimensions,
%   and K is the number of histogram bins. INTHIST may be of class
%   UINT32 or DOUBLE.
%
%   Each box is a four dimensional vector [IMIN JMIN IMAX JMAX]' of
%   class UINT32 and correspond to the index set [IMIN, IMAX] x [JMIN,
%   JMAX]. To specify an empty box, let IMIN > IMAX.
%
%   HISTS stores one histogram per column (one for each box) and has K
%   rows, one for each histogram bin. HIST is of the same class of
%   INTHIST.
%
%   See also: VL_INTHIST(), VL_IMINTEGRAL(), VL_HELP().

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
function updateSSHypothesis( aid )
%UPDATESSHYPOTHESIS Summary of this function goes here
%   Detailed explanation goes here
%matlabpool(2);
[~,b] = system('hostname');
fprintf('Running on %s\n',b);

% diaryfilename = sprintf('../diary/diary%03d.txt', aid);
% diary(diaryfilename);
% diary on;

addpath(genpath('/n/fs/sun360/panoParser/SceneParsing'));
cd '/n/fs/sun360/panoParser/SceneParsing';
folderName = './annotation/dataset/bedroom/';
% folderName = './annotation/dataset/living_room/';

load([folderName '/IMGLIST.mat']);
load([folderName '/ANNO_ALL.mat']);
fprintf('Data loaded\n');
testfile2 = [folderName IMGLIST(aid).name '/bottomupcuboid_P1.mat'];
if exist(testfile2, 'file')
	fprintf('Job: %d: Previously done^_^\n', aid);
	return;
end

anno_valid = false( 1, length(ANNO_ALL));
for i = 1:length(anno_valid)
    anno_valid(i) = ~isempty(ANNO_ALL(i).ANNO3D) && ANNO_ALL(i).ANNO3D.b_singleroom;
end

if ~anno_valid(aid)
    return;
end

rotImg = imread([IMGLIST(aid).small_imgname(1:end-9) 'rotate.png']);
fprintf('Image loaded\n');

SS = getSelectiveSearch(rotImg);

if exist(testfile2, 'file')
	fprintf('Job: %d: Previously done^_^\n', aid);
	return;
end

AllCuboid = regionBasedHypothesisFromSS(SS, 1024, 2048);

if exist(testfile2, 'file')
	fprintf('Job: %d: Previously done^_^\n', aid);
	return;
end

rectangle = rectDetectionFlexFromSS(SS, 1024, 2048);

if exist(testfile2, 'file')
	fprintf('Job: %d: Previously done^_^\n', aid);
	return;
end

AllCuboid_NEW = pruneCuboid(rotImg, AllCuboid);
rectangle_NEW = pruneRectangle(rotImg, rectangle);

load([folderName IMGLIST(aid).name '/bottomupcuboid.mat']);
fprintf('Old data loaded\n');

regionCuboid.views = [regionCuboid.views; AllCuboid_NEW.views];
regionCuboid.xyzBox = [regionCuboid.xyzBox; AllCuboid_NEW.xyzBox];
regionCuboid.score = [regionCuboid.score; AllCuboid_NEW.score];
regionCuboid.count = regionCuboid.count + AllCuboid_NEW.count;

for i = 1:6
    rectangle(i).highPixelBox = [rectangle(i).highPixelBox; rectangle_NEW(i).highPixelBox];
%     rectangle(i).count = rectangle(i).count + rectangle_NEW(i).count;
    rectangle(i).xyzBox = [rectangle(i).xyzBox; rectangle_NEW(i).xyzBox];
    rectangle(i).score = [rectangle(i).score; rectangle_NEW(i).score];
end

save([folderName IMGLIST(aid).name '/bottomupcuboid_P1.mat'], ...
            'rectangle', 'rectCuboid2', 'rectCuboid3', 'regionCuboid');
fprintf('save data for bottomcuboid_P1.mat\n');

%% move updatePruneRectDetBySegmentation
testfile2 = [folderName IMGLIST(aid).name '/rectangleSegScore_P1.mat'];
[ rect_score ] = getSegScrOfRect( rotImg, rectangle );
save(testfile2, 'rect_score');
fprintf('save data for rectangleSegScore_P1.mat\n');

end

author: zouchuhang

repo: LayoutNet


text:
classdef TestConvexHull
    %TestConvexHull

    methods (Static)
        function test_input_1
            % cell array
            points = {[0 0], [1 1], [2 2], [3 3], [4 4]};
            hull = cv.convexHull(points);
        end

        function test_input_2
            % Nx2 matrix
            points = [0 0; 1 0; 2 2; 3 3; 3 4];
            hull = cv.convexHull(points);
        end

        function test_input_3
            % 1xNx2 matrix
            points = shiftdim([0 0; 1 0; 2 2; 3 3; 3 4], -1);
            hull = cv.convexHull(points);
        end

        function test_input_4
            % Nx1x2 matrix
            points = permute([0 0; 1 0; 2 2; 3 3; 3 4], [1 3 2]);
            hull = cv.convexHull(points);
        end

        function test_type_1
            points = single([0 0; 1 0; 2 2; 3 3; 3 4]);
            hull = cv.convexHull(points);
        end

        function test_type_2
            points = int32([0 0; 1 0; 2 2; 3 3; 3 4]);
            hull = cv.convexHull(points);
        end

        function test_option_1
            points = [0 0; 1 0; 2 2; 3 3; 3 4];
            hull = cv.convexHull(points, 'Clockwise',false);
            hull = cv.convexHull(points, 'Clockwise',true);
        end

        function test_output_1
            points = int32([0 0; 1 0; 2 2; 3 3; 3 4]);

            hull = cv.convexHull(points, 'ReturnPoints',true);
            validateattributes(hull, {'numeric'}, {'2d', 'ncols',2});
            assert(isa(hull, class(points)));

            hull = cv.convexHull(points, 'ReturnPoints',false);
            validateattributes(hull, {'int32'}, {'vector', 'nonnegative'});
        end

        function test_output_2
            points = num2cell(int32([0 0; 1 0; 2 2; 3 3; 3 4]),2);

            hull = cv.convexHull(points, 'ReturnPoints',true);
            validateattributes(hull, {'cell'}, {'vector'});
            assert(all(cellfun(@numel, hull) == 2));

            hull = cv.convexHull(points, 'ReturnPoints',false);
            validateattributes(hull, {'int32'}, {'vector', 'nonnegative'});
        end

        function test_error_argnum
            try
                cv.convexHull();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
% @authors:     Ahmad Humayun
% @contact:     ahumayun@cc.gatech.edu
% @affiliation: Georgia Institute of Technology
% @date:        Fall 2013 - Summer 2014

function [nonlambda_s, nonlambda_t, lambda_s, lambda_t, lambda_range, ...
          pairwise_edges, DISC_FACTOR, BIG_VALUE ] = ...
            preprocess_graphs(gp_obj, nonlambda_s, nonlambda_t, ...
                              lambda_s, lambda_t)
    l = -1;
    u = gp_obj.graph_sol_upper_bp;

    DISC_FACTOR = 1000; % original is 1000

    % i don't think it can be bigger (beyond certain value goes crazy)
    BIG_VALUE = 21475000000;

    % get the range of lambda values to use
    pmc_num_lambdas = gp_obj.segm_params.pmc_num_lambdas;

    if ~exist('overrride_lambdas', 'var')
        overrride_lambdas = [];
    end

    lambda_range = ...
        GraphProcessor.generate_param_lambdas(l, u, DISC_FACTOR, ...
                                              pmc_num_lambdas, ...
                                              overrride_lambdas);

    % capping lambda weight and offsets (parametric capacities) to BIG_VALUE
    nonlambda_s = nonlambda_s * DISC_FACTOR;
    nonlambda_t = nonlambda_t * DISC_FACTOR;
    nonlambda_s(nonlambda_s > BIG_VALUE) = BIG_VALUE;
    nonlambda_t(nonlambda_t > BIG_VALUE) = BIG_VALUE;

    % capping non-lambda weights (non-parametric capacities) to BIG_VALUE
    lambda_s = lambda_s * DISC_FACTOR;
    lambda_t = lambda_t * DISC_FACTOR;
    lambda_s(lambda_s > BIG_VALUE) = BIG_VALUE;
    lambda_t(lambda_t > BIG_VALUE) = BIG_VALUE;

    pairwise_vals = gp_obj.pairwise_vals;
    pairwise_vals = pairwise_vals * DISC_FACTOR;
    pairwise_vals(pairwise_vals > BIG_VALUE) = BIG_VALUE;

    % a < from-node > < to-node > < capacity >
    pairwise_edges = [gp_obj.pairwise_edges pairwise_vals];
end

author: Cloud-CV

repo: object-proposals


text:
function [cc] = oz2cc(oz)
% Convert volume from US liquid ounces to cubic centimeters.
% Chad Greene 2012
cc = oz*29.573529563;


text:
function [N,rho,T,P]=standardAtmosParam(dayOfYear,secondOfTheDay,latLonAlt,rhow)
%%STANDARDATMOSPARAM  Get basic parameters for atmospheric refractivity,
%                     density, pressure and temperature from the
%                     NRLMSISE-00 atmospheric model. The model is best
%                     suited for altitudes below 90km as the anomalous
%                     oxygen parameter of the NRLMSISE-00 model is not
%                     used by this function.
%
%INPUTS: dayOfYear The integer day of the year in the Gregorian calendar in
%                  universal coordinated time (UTC). Counting starts at 1.
%                  The resolution of the model is not sufficient for it to
%                  matter whether 365 or 366 is given as the day if it
%                  isn't/is a leap year.
%   secondOfTheDay The second of the day. This starts at zero. The
%                  resolution of the model is not high enough for leap
%                  seconds to matter, so values above 86400.0 are just
%                  clipped to 86400.0.
%        latLonAlt The location under consideration given in WGS-84
%                  ellipsoidal coordinates of latitude and longitude in
%                  radians and ellipsoidal height in meters.
%             rhow An optional parameter specifying the mass density of
%                  water vapor at the point in question in units of
%                  kilograms per cubic meter. If omitted, the air is
%                  assumed to be dry (rhow=0). The total density of the air
%                  is assumed to be the sum of the dry air density and
%                  rhow.
%
%OUTPUTS: N The refractivity of the atmosphere. In this model, N is always
%           real. N=10^6*(n-1) where n is the index of refraction. This is
%           generally valid for frequencies from L-band (1GHz) to 10 GHz
%           (the middle of X-band).
%       rho The atmospheric density at the point in question in units of
%           kilograms per cubic meter.
%         T The temperature at the point in question with units of degrees
%           Kelvin.
%         P The atmospheric pressure at the point in question in units of
%           Newtons per square meter (Pascals). It assumes that the gasses
%           can be treated as ideal gasses.
%
%The dry air density and temperature is obtained from the NRLMSISE-00
%atmospheric model using the default parameters for magnetic and solar
%parameters, which are generally valid below 90km, using the function
%NRLMSISE00GasTemp. The reslting standard atmospheric parameters are then
%passed to the atmosParam4GasTemp function.
%
%The atmospheric pressure is obtained using the Ideal Gas Law as described
%in [1]. The properties of atmospheric constituents used in
%atmosParam4GasTemp are more recent than those used in the NRLMSISE-00
%standard. Consequently, pressure obtained for the associated altitude is
%not completely numerically consistent with the function
%NRLMSISE00Alt4Pres.
%
%REFERENCES:
%[1] D. P. Drob, THE NRLMSISE-00 AND HWM-93 USERS GUIDE: Version 1.50,
%    Nov. 2003.
%
%December 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<4)
   rhow=0;
end

[gasTable,t]=NRLMSISE00GasTemp(dayOfYear,secondOfTheDay,latLonAlt);
T=t(2);%The temperature in Kelving at the point.

[N,rho,P]=atmosParam4GasTemp(gasTable,T,rhow);
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function pass = test_cumsum(pref)

if ( nargin == 0 )
    pref = chebfunpref;
end
tol = 100*pref.cheb3Prefs.chebfun3eps;

% Check cumsum on cube domain
[x, y, z] = cheb.xyz;

pass(1) = norm(cumsum(x) - .5*(x.^2-1)) < tol;
pass(2) = norm(cumsum(x, 1) - .5*(x.^2-1)) < tol;
pass(3) = norm(cumsum(x, 2) - x.*(y+1)) < tol;
pass(4) = norm(cumsum(x, 3) - x.*(z+1)) < tol;

pass(5) = norm(cumsum(y) - y.*(x+1)) < tol;
pass(6) = norm(cumsum(y, 1) - y.*(x+1)) < tol;
pass(7) = norm(cumsum(y, 2) - .5*(y.^2-1)) < tol;
pass(8) = norm(cumsum(y, 3) - y.*(z+1)) < tol;

pass(9) = norm(cumsum(z) - z.*(x+1)) < tol;
pass(10) = norm(cumsum(z, 1) - z.*(x+1)) < tol;
pass(11) = norm(cumsum(z, 2) - z.*(y+1)) < tol;
pass(12) = norm(cumsum(z, 3) - .5*(z.^2-1)) < tol;

% Check cumsum on rectangular box domain
x = chebfun3(@(x,y,z) x, [-1.1  2 -0.2 3 5 6]);
y = chebfun3(@(x,y,z) y, [-1.1  2 -0.2 3 5 6]);
z = chebfun3(@(x,y,z) z, [-1.1  2 -0.2 3 5 6]);

pass(13) = norm(cumsum(x) - .5*(x.^2-(-1.1)^2)) < tol;
pass(14) = norm(cumsum(x, 1) - .5*(x.^2-(-1.1)^2)) < tol;
pass(15) = norm(cumsum(x, 2) - x.*(y+0.2)) < tol;
pass(16) = norm(cumsum(x, 3) - x.*(z-5)) < tol;

pass(17) = norm(cumsum(y) - y.*(x+1.1)) < tol;
pass(18) = norm(cumsum(y, 1) - y.*(x+1.1)) < tol;
pass(19) = norm(cumsum(y, 2) - .5*(y.^2-(-0.2)^2)) < tol;
pass(20) = norm(cumsum(y, 3) - y.*(z-5)) < tol;

pass(21) = norm(cumsum(z) - z.*(x+1.1)) < tol;
pass(22) = norm(cumsum(z, 1) - z.*(x+1.1)) < tol;
pass(23) = norm(cumsum(z, 2) - z.*(y+0.2)) < tol;
pass(24) = norm(cumsum(z, 3) - .5*(z.^2-5^2)) < tol;

% Check that a triple cumsum is a cumsum3
f = sin((x-.1).*(y+.1).*(z+.1));
pass(25) = norm(cumsum(cumsum(cumsum(f), 2), 3) - cumsum3(f)) < tol;
pass(26) = norm(cumsum(cumsum(cumsum(f, 1), 2), 3) - cumsum3(f)) < tol;

% Look in one direction and make sure we get the right thing:
f = chebfun(@(x) exp(x));
f3 = chebfun3(@(x,y,z) exp(x));
g = cumsum(f);
g3 = cumsum(f3);
g3x = g3(:,0,.3);
pass(27) = norm(g3x-g) < tol;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function x = sindeg(degrees)
%sin of angle in degrees.
%JOD

rad = (degrees/180) * pi;

x = sin(rad);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function L = resymbol (L, A)						    %#ok
%RESYMBOL recomputes the symbolic Cholesky factorization of the matrix A.
%
%   Example:
%   L = resymbol (L, A)
%
%   Recompute the symbolic Cholesky factorization of the matrix A.  A must be
%   symmetric.  Only tril(A) is used.  Entries in L that are not in the Cholesky
%   factorization of A are removed from L.  L can be from an LL' or LDL'
%   factorization (lchol or ldlchol).  resymbol is useful after a series of
%   downdates via ldlupdate, since downdates do not remove any entries in L.
%   The numerical values of A are ignored; only its nonzero pattern is used.
%
% See also LCHOL, LDLUPDATE

%   Copyright 2006-2007, Timothy A. Davis
%   http://www.cise.ufl.edu/research/sparse

error ('resymbol not found') ;

author: fieldtrip

repo: fieldtrip


text:
function cr = slcorrectrate(scores, clabels, qlabels, op, varargin)
%SLCORRECTRATE Computes the correct rate of classification
%
% $ Syntax $
%   - cr = slcorrectrate(scores, clabels, qlabels, op, ...)
%
% $ Arguments $
%   - scores:           the scores to support the classification
%   - clabels:          the labels of classes
%   - qlabels:          the groundtruth of the labels of query samples
%   - op:               the option of the score
%   - cr:               the correct rate of the score-based classification
%
% $ Description $
%   - cr = slcorrectrate(scores, clabels, slabels, op, ...)
%     Computes the classification correct rate. Suppose we want to
%     classify n samples into m classes, then scores will be an m x n
%     matrix, with the entry at i-th row, j-th column representing the
%     score of the j-th sample in the i-th class.
%     For op, it can take either of the two values: 'high' and 'low'.
%     If op is 'high' then it means that the sample will be classified
%     to the class where it has the highest score value, vice versa
%     for 'low'.
%
% $ History $
%   - Created by Dahua Lin on Jun 10th, 2005
%   - Modified by Dahua Lin on May 1st, 2005
%     - To base on the sltoolbox v4.
%   - Modified by Dahua Lin on Aug 9th, 2006
%     - Extract slclassify as independent function and base on it
%   - Modified by Dahua Lin on Aug 16th, 2006
%     - Based on new slclassify to support multiple schemes
%

%% parse and verify input arguments

if nargin < 4
    raise_lackinput('slcorrectrate', 4);
end

%% Make decision

decisions = slclassify(scores, clabels, op, varargin{:});

%% Evaluate correct rate

qlabels = qlabels(:)';
cr = sum(decisions == qlabels) / length(qlabels);

author: lmthang


text:
function CVaR=ComputeCVaR(Units,Scenarios,Conf)

PnL=Scenarios*Units;
Sort_PnL=sort(PnL);

J=length(PnL);
Cut=round(J*(1-Conf));

CVaR=-mean(Sort_PnL(1:Cut));


text:
function lik = lik_qgp(varargin)
%LIK_QGP  Create a Quantile Gaussian Process likelihood (utility) structure
%
%  Description
%    LIK = LIK_QGP('PARAM1',VALUE1,'PARAM2,VALUE2,...)
%    creates a quantile gp likelihood structure in which the named
%    parameters have the specified values. Any unspecified
%    parameters are set to default values.
%
%    LIK = LIK_QGP(LIK,'PARAM1',VALUE1,'PARAM2,VALUE2,...)
%    modify a likelihood function structure with the named
%    parameters altered with the specified values.
%
%    Parameters for QGP likelihood function [default]
%      sigma2       - variance [0.1]
%      sigma2_prior - prior for sigma2 [prior_logunif]
%      quantile     - Quantile of interest [0.5]
%
%    Note! If the prior is 'prior_fixed' then the parameter in
%    question is considered fixed and it is not handled in
%    optimization, grid integration, MCMC etc.
%
%    The likelihood is defined as follows:
%                            __ n
%      p(y|f, sigma2, tau) = || i=1 tau*(1-tau)/sigma*exp(-(y-f)/sigma*
%                                 (tau - I(t <= f)))
%
%    where tau is the quantile of interest, sigma is the standard deviation
%    of the distribution and I(t <= f) = 1 if t <= f, 0 otherwise.
%
%    Note that because the form of the likelihood, second order derivatives
%    with respect to latent values are 0. Because this, EP should be used
%    instead of Laplace approximation.
%
%  See also
%    GP_SET, PRIOR_*, LIK_*
%
%   References
%     Boukouvalas et al. (2012). Direct Gaussian Process Quantile Regression
%     Using Expectation Propagation. Appearing in Proceedings of the 29th
%     International Conference on Machine Learning, Edinburg, Scotland, UK,
%     2012.
%

% Copyright (c) 2012 Ville Tolvanen

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

  ip=inputParser;
  ip.FunctionName = 'LIK_QGP';
  ip.addOptional('lik', [], @isstruct);
  ip.addParamValue('sigma2',0.1, @(x) isscalar(x) && x>0);
  ip.addParamValue('sigma2_prior',prior_logunif(), @(x) isstruct(x) || isempty(x));
  ip.addParamValue('quantile',0.5, @(x) isscalar(x) && x>0 && x<1);
  ip.parse(varargin{:});
  lik=ip.Results.lik;

  if isempty(lik)
    init=true;
    lik.type = 'QGP';
  else
    if ~isfield(lik,'type') || ~isequal(lik.type,'QGP')
      error('First argument does not seem to be a valid likelihood function structure')
    end
    init=false;
  end

  % Initialize parameters
  if init || ~ismember('sigma2',ip.UsingDefaults)
    lik.sigma2 = ip.Results.sigma2;
  end
  if init || ~ismember('quantile',ip.UsingDefaults)
    lik.quantile = ip.Results.quantile;
  end
  % Initialize prior structure
  if init
    lik.p=[];
  end
  if init || ~ismember('sigma2_prior',ip.UsingDefaults)
    lik.p.sigma2=ip.Results.sigma2_prior;
  end
  if init
    % Set the function handles to the subfunctions
    lik.fh.pak = @lik_qgp_pak;
    lik.fh.unpak = @lik_qgp_unpak;
    lik.fh.lp = @lik_qgp_lp;
    lik.fh.lpg = @lik_qgp_lpg;
    lik.fh.ll = @lik_qgp_ll;
    lik.fh.llg = @lik_qgp_llg;
    lik.fh.llg2 = @lik_qgp_llg2;
    lik.fh.llg3 = @lik_qgp_llg3;
    lik.fh.tiltedMoments = @lik_qgp_tiltedMoments;
    lik.fh.siteDeriv = @lik_qgp_siteDeriv;
    lik.fh.predy = @lik_qgp_predy;
    lik.fh.invlink = @lik_qgp_invlink;
    lik.fh.recappend = @lik_qgp_recappend;
  end

end

function [w s] = lik_qgp_pak(lik)
%LIK_QGP_PAK  Combine likelihood parameters into one vector.
%
%  Description
%    W = LIK_QGP_PAK(LIK) takes a likelihood structure LIK
%    and combines the parameters into a single row vector W.
%    This is a mandatory subfunction used for example in
%    energy and gradient computations.
%
%       w = [ log(lik.sigma2)
%             (hyperparameters of lik.magnSigma2)]'
%
%  See also
%    LIK_QGP_UNPAK

  w = []; s = {};
  if ~isempty(lik.p.sigma2)
    w = [w log(lik.sigma2)];
    s = [s; 'log(qgp.sigma2)'];
    % Hyperparameters of sigma2
    [wh sh] = lik.p.sigma2.fh.pak(lik.p.sigma2);
    w = [w wh];
    s = [s; sh];
  end

end

function [lik, w] = lik_qgp_unpak(lik, w)
%LIK_QGP_UNPAK  Extract likelihood parameters from the vector.
%
%  Description
%    W = LIK_QGP_UNPAK(W, LIK) takes a likelihood structure
%    LIK and extracts the parameters from the vector W to the LIK
%    structure. This is a mandatory subfunction used for example
%    in energy and gradient computations.
%
%    Assignment is inverse of
%       w = [ log(lik.sigma2)
%             (hyperparameters of lik.magnSigma2)]'
%
%  See also
%    LIK_QGP_PAK

  if ~isempty(lik.p.sigma2)
    lik.sigma2 = exp(w(1));
    w = w(2:end);

    % Hyperparameters of sigma2
    [p, w] = lik.p.sigma2.fh.unpak(lik.p.sigma2, w);
    lik.p.sigma2 = p;
  end
end

function lp = lik_qgp_lp(lik)
%LIK_QGP_LP  Evaluate the log prior of likelihood parameters
%
%  Description
%    LP = LIK_QGP_LP(LIK) takes a likelihood structure LIK and
%    returns log(p(th)), where th collects the parameters. This
%    subfunction is needed when there are likelihood parameters.
%
%  See also
%    LIK_QGP_PAK, LIK_QGP_UNPAK, LIK_QGP_G, GP_E

  lp = 0;

  if ~isempty(lik.p.sigma2)
    likp=lik.p;
    lp = likp.sigma2.fh.lp(lik.sigma2, likp.sigma2) + log(lik.sigma2);
  end
end

function lpg = lik_qgp_lpg(lik)
%LIK_QGP_LPG  Evaluate gradient of the log prior with respect
%                  to the parameters.
%
%  Description
%    LPG = LIK_QGP_LPG(LIK) takes a QGP likelihood
%    function structure LIK and returns LPG = d log (p(th))/dth,
%    where th is the vector of parameters. This subfunction is
%    needed when there are likelihood parameters.
%
%  See also
%    LIK_QGP_PAK, LIK_QGP_UNPAK, LIK_QGP_E, GP_G

  lpg = [];

  if ~isempty(lik.p.sigma2)
    likp=lik.p;

    lpgs = likp.sigma2.fh.lpg(lik.sigma2, likp.sigma2);
    lpg = lpgs(1).*lik.sigma2 + 1;
    if length(lpgs) > 1
      lpg = [lpg lpgs(2:end)];
    end
  end
end

function ll = lik_qgp_ll(lik, y, f, z)
%LIK_QGP_LL  Log likelihood
%
%  Description
%    LL = LIK_QGP_LL(LIK, Y, F, Z) takes a likelihood
%    structure LIK, observations Y and latent values F.
%    Returns the log likelihood, log p(y|f,z). This subfunction
%    is needed when using Laplace approximation or MCMC for
%    inference with non-Gaussian likelihoods. This subfunction
%    is also used in information criteria (DIC, WAIC) computations.
%
%  See also
%    LIK_QGP_LLG, LIK_QGP_LLG3, LIK_QGP_LLG2, GPLA_E

  tau=lik.quantile;
  sigma=sqrt(lik.sigma2);
  ll = sum(log(tau*(1-tau)/sigma) - (y-f)./sigma.*(tau-(y<=f)));
end

function llg = lik_qgp_llg(lik, y, f, param, z)
%LIK_QGP_LLG  Gradient of the log likelihood
%
%  Description
%    LLG = LIK_QGP_LLG(LIK, Y, F, PARAM) takes a likelihood
%    structure LIK, observations Y and latent values F. Returns
%    the gradient of the log likelihood with respect to PARAM.
%    At the moment PARAM can be 'param' or 'latent'. This subfunction
%    is needed when using Laplace approximation or MCMC for inference
%    with non-Gaussian likelihoods.
%
%  See also
%    LIK_QGP_LL, LIK_QGP_LLG2, LIK_QGP_LLG3, GPLA_E

  tau=lik.quantile;
  sigma2=sqrt(lik.sigma2);
  switch param
    case 'param'
      llg = sum(-1/(2.*sigma2) + (y-f)./(2.*sigma2^(3/2)).*(tau-(y<=f)));

      % correction for the log transformation
      llg = llg.*lik.sigma2;
    case 'latent'
      llg = (tau-(y<=f))/sqrt(sigma2);
  end
end

function llg2 = lik_qgp_llg2(lik, y, f, param, z)
%LIK_QGP_LLG2  Second gradients of the log likelihood
%
%  Description
%    LLG2 = LIK_QGP_LLG2(LIK, Y, F, PARAM) takes a likelihood
%    structure LIK, observations Y and latent values F. Returns
%    the Hessian of the log likelihood with respect to PARAM.
%    At the moment PARAM can be 'param' or 'latent'. LLG2 is
%    a vector with diagonal elements of the Hessian matrix
%    (off diagonals are zero). This subfunction is needed
%    when using Laplace approximation or EP for inference
%    with non-Gaussian likelihoods.
%
%  See also
%    LIK_QGP_LL, LIK_QGP_LLG, LIK_QGP_LLG3, GPLA_E

  tau=lik.quantile;
  sigma2=lik.sigma2;
  switch param
    case 'param'
      llg2 = sum(1/(2*sigma2^2) - 3.*(tau-(y<=f)).*(y-f)./(4.*sigma2^(5/2)));

      % correction due to the log transformation
      llg2 = llg2.*lik.sigma2;
    case 'latent'
      llg2 = zeros(size(f));
    case 'latent+param'
      llg2 = -(tau-(y<=f))./(2*sigma2^(3/2));

      % correction due to the log transformation
      llg2 = llg2.*lik.disper;
  end
end

function llg3 = lik_qgp_llg3(lik, y, f, param, z)
%LIK_QGP_LLG3  Third gradients of the log likelihood
%
%  Description
%    LLG3 = LIK_QGP_LLG3(LIK, Y, F, PARAM) takes a likelihood
%    structure LIK, observations Y and latent values F and
%    returns the third gradients of the log likelihood with
%    respect to PARAM. At the moment PARAM can be 'param' or
%    'latent'. LLG3 is a vector with third gradients. This
%    subfunction is needed when using Laplace approximation for
%    inference with non-Gaussian likelihoods.
%
%  See also
%    LIK_QGP_LL, LIK_QGP_LLG, LIK_QGP_LLG2, GPLA_E, GPLA_G

  tau=lik.quantile;
  sigma2=lik.sigma2;
  switch param
    case 'param'
      llg3 = sum(-1/sigma2^3 + 15.*(tau-(y<=f)).*(y-f)./(8.*sigma2^(7/2)));
    case 'latent'
      llg3 = 0;
    case 'latent2+param'
      llg3 = 0;

      % correction due to the log transformation
      llg3 = llg3.*lik.sigma2;
  end
end

function [logM_0, m_1, sigm2hati1] = lik_qgp_tiltedMoments(lik, y, i1, sigm2_i, myy_i, z)
%LIK_QGP_TILTEDMOMENTS  Returns the marginal moments for EP algorithm
%
%  Description
%    [M_0, M_1, M2] = LIK_QGP_TILTEDMOMENTS(LIK, Y, I, S2,
%    MYY, Z) takes a likelihood structure LIK, observations
%    Y, index I and cavity variance S2 and mean MYY. Returns
%    the zeroth moment M_0, mean M_1 and variance M_2 of the
%    posterior marginal (see Rasmussen and Williams (2006):
%    Gaussian processes for Machine Learning, page 55). This
%    subfunction is needed when using EP for inference with
%    non-Gaussian likelihoods.
%
%  See also
%    GPEP_E

  yy = y(i1);
  sigma2 = lik.sigma2;
  tau=lik.quantile;
  logM_0=zeros(size(yy));
  m_1=zeros(size(yy));
  sigm2hati1=zeros(size(yy));

  for i=1:length(i1)
    % get a function handle of an unnormalized tilted distribution
    % (likelihood * cavity = Quantile-GP * Gaussian)
    % and useful integration limits
    [tf,minf,maxf]=init_qgp_norm(yy(i),myy_i(i),sigm2_i(i),sigma2,tau);

    % Integrate with quadrature
    RTOL = 1.e-6;
    ATOL = 1.e-10;
    [m_0, m_1(i), m_2] = quad_moments(tf, minf, maxf, RTOL, ATOL);
    sigm2hati1(i) = m_2 - m_1(i).^2;

    % If the second central moment is less than cavity variance
    % integrate more precisely. Theoretically for log-concave
    % likelihood should be sigm2hati1 < sigm2_i.
    if sigm2hati1(i) >= sigm2_i(i)
      ATOL = ATOL.^2;
      RTOL = RTOL.^2;
      [m_0, m_1(i), m_2] = quad_moments(tf, minf, maxf, RTOL, ATOL);
      sigm2hati1(i) = m_2 - m_1(i).^2;
      if sigm2hati1(i) >= sigm2_i(i)
        error('lik_qgp_tilted_moments: sigm2hati1 >= sigm2_i');
      end
    end
    logM_0(i) = log(m_0);
  end
end

function [g_i] = lik_qgp_siteDeriv(lik, y, i1, sigm2_i, myy_i, z)
%LIK_QGP_SITEDERIV  Evaluate the expectation of the gradient
%                      of the log likelihood term with respect
%                      to the likelihood parameters for EP
%
%  Description [M_0, M_1, M2] =
%    LIK_QGP_SITEDERIV(LIK, Y, I, S2, MYY, Z) takes a
%    likelihood structure LIK, observations Y, index I
%    and cavity variance S2 and mean MYY. Returns E_f
%    [d log p(y_i|f_i) /d a], where a is the likelihood
%    parameter and the expectation is over the marginal posterior.
%    This term is needed when evaluating the gradients of
%    the marginal likelihood estimate Z_EP with respect to
%    the likelihood parameters (see Seeger (2008):
%    Expectation propagation for exponential families).This
%    subfunction is needed when using EP for inference with
%    non-Gaussian likelihoods and there are likelihood parameters.
%
%  See also
%    GPEP_G

  yy = y(i1);
  sigma2=lik.sigma2;
  tau=lik.quantile;

  % get a function handle of an unnormalized tilted distribution
  % (likelihood * cavity = Quantile-GP * Gaussian)
  % and useful integration limits
  [tf,minf,maxf]=init_qgp_norm(yy,myy_i,sigm2_i,sigma2,tau);
  % additionally get function handle for the derivative
  td = @deriv;

  % Integrate with quadgk
  [m_0, fhncnt] = quadgk(tf, minf, maxf);
  [g_i, fhncnt] = quadgk(@(f) td(f).*tf(f)./m_0, minf, maxf);
  g_i = g_i.*sigma2;

  function g = deriv(f)

    g = -1/(2.*sigma2) + (yy-f)./(2.*sigma2^(3/2)).*(tau-(yy<=f));

  end
end

function [lpy, Ey, Vary] = lik_qgp_predy(lik, Ef, Varf, yt, zt)
%LIK_QGP_PREDY  Returns the predictive mean, variance and density of y
%
%  Description
%    LPY = LIK_QGP_PREDY(LIK, EF, VARF YT, ZT)
%    Returns logarithm of the predictive density PY of YT, that is
%        p(yt | zt) = \int p(yt | f, zt) p(f|y) df.
%    This subfunction is needed when computing posterior predictive
%    distributions for future observations.
%
%    [LPY, EY, VARY] = LIK_QGP_PREDY(LIK, EF, VARF) takes a
%    likelihood structure LIK, posterior mean EF and posterior
%    Variance VARF of the latent variable and returns the
%    posterior predictive mean EY and variance VARY of the
%    observations related to the latent variables. This
%    subfunction is needed when computing posterior predictive
%    distributions for future observations.
%

%
%  See also
%    GPLA_PRED, GPEP_PRED, GPMC_PRED

  sigma2=lik.sigma2;
  tau=lik.quantile;

  Ey=[];
  Vary=[];

  % Evaluate the posterior predictive densities of the given observations
  lpy = zeros(length(yt),1);
  for i1=1:length(yt)
    % get a function handle of the likelihood times posterior
    % (likelihood * posterior = Quantile-GP * Gaussian)
    % and useful integration limits
    [pdf,minf,maxf]=init_qgp_norm(...
      yt(i1),Ef(i1),Varf(i1),sigma2, tau);
    % integrate over the f to get posterior predictive distribution
    lpy(i1) = log(quadgk(pdf, minf, maxf));
  end
end

function [df,minf,maxf] = init_qgp_norm(yy,myy_i,sigm2_i,sigma2,tau)
%INIT_QGP_NORM
%
%  Description
%    Return function handle to a function evaluating
%    Quantile-GP * Gaussian which is used for evaluating
%    (likelihood * cavity) or (likelihood * posterior) Return
%    also useful limits for integration. This is private function
%    for lik_qgp. This subfunction is needed by subfunctions
%    tiltedMoments, siteDeriv and predy.
%
%  See also
%    LIK_QGP_TILTEDMOMENTS, LIK_QGP_SITEDERIV,
%    LIK_QGP_PREDY

  sigma=sqrt(sigma2);
% avoid repetitive evaluation of constant part
  ldconst = log(tau*(1-tau)/sigma) ...
            - log(sigm2_i)/2 - log(2*pi)/2;
  % Create function handle for the function to be integrated
  df = @qgp_norm;
  % use log to avoid underflow, and derivates for faster search
  ld = @log_qgp_norm;
  ldg = @log_qgp_norm_g;
%   ldg2 = @log_qgp_norm_g2;

  % Set the limits for integration
  % Quantile-GP likelihood is log-concave so the qgp_norm
  % function is unimodal, which makes things easier
  if yy==0
    % with yy==0, the mode of the likelihood is not defined
    % use the mode of the Gaussian (cavity or posterior) as a first guess
    modef = myy_i;
  else
    % use precision weighted mean of the Gaussian approximation
    % of the Quantile-GP likelihood and Gaussian
    modef = (myy_i/sigm2_i + yy/sigma2)/(1/sigm2_i + 1/sigma2);
  end
  % find the mode of the integrand using Newton iterations
  % few iterations is enough, since the first guess in the right direction
  niter=8;       % number of Newton iterations

  minf=modef-6*sigm2_i;
  while ldg(minf) < 0
    minf=minf-2*sigm2_i;
  end
  maxf=modef+6*sigm2_i;
  while ldg(maxf) > 0
    maxf=maxf+2*sigm2_i;
  end
  for ni=1:niter
%     h=ldg2(modef);
    modef=0.5*(minf+maxf);
    if ldg(modef) < 0
      maxf=modef;
    else
      minf=modef;
    end
  end
  % integrand limits based on Gaussian approximation at mode
  minf=modef-6*sqrt(sigm2_i);
  maxf=modef+6*sqrt(sigm2_i);
  modeld=ld(modef);
  iter=0;
  % check that density at end points is low enough
  lddiff=20; % min difference in log-density between mode and end-points
  minld=ld(minf);
  step=1;
  while minld>(modeld-lddiff)
    minf=minf-step*sqrt(sigm2_i);
    minld=ld(minf);
    iter=iter+1;
    step=step*2;
    if iter>100
      error(['lik_qgp -> init_qgp_norm: ' ...
             'integration interval minimun not found ' ...
             'even after looking hard!'])
    end
  end
  maxld=ld(maxf);
  step=1;
  while maxld>(modeld-lddiff)
    maxf=maxf+step*sqrt(sigm2_i);
    maxld=ld(maxf);
    iter=iter+1;
    step=step*2;
    if iter>100
      error(['lik_qgp -> init_qgp_norm: ' ...
             'integration interval maximun not found ' ...
             'even after looking hard!'])
    end
  end

  function integrand = qgp_norm(f)
  % Quantile-GP * Gaussian
    integrand = exp(ldconst ...
                    -(yy-f)./sqrt(sigma2).*(tau-(yy<=f)) ...
                    -0.5*(f-myy_i).^2./sigm2_i);
  end

  function log_int = log_qgp_norm(f)
  % log(Quantile-GP * Gaussian)
  % log_qgp_norm is used to avoid underflow when searching
  % integration interval
    log_int = ldconst...
              -(yy-f)./sqrt(sigma2).*(tau-(yy<=f)) ...
              -0.5*(f-myy_i).^2./sigm2_i;
  end

  function g = log_qgp_norm_g(f)
  % d/df log(Quantile-GP * Gaussian)
  % derivative of log_qgp_norm
    g = (tau-(yy<=f))/sqrt(sigma2) ...
        + (myy_i - f)./sigm2_i;
  end

end

function mu = lik_qgp_invlink(lik, f, z)
%LIK_QGP_INVLINK  Returns values of inverse link function
%
%  Description
%    MU = LIK_QGP_INVLINK(LIK, F) takes a likelihood structure LIK and
%    latent values F and returns the values MU of inverse link function.
%    This subfunction is needed when using function gp_predprctmu.
%
%     See also
%     LIK_QGP_LL, LIK_QGP_PREDY

  mu = f;
end

function reclik = lik_qgp_recappend(reclik, ri, lik)
%RECAPPEND  Append the parameters to the record
%
%  Description
%    RECLIK = LIK_QGP_RECAPPEND(RECLIK, RI, LIK) takes a
%    likelihood record structure RECLIK, record index RI and
%    likelihood structure LIK with the current MCMC samples of
%    the parameters. Returns RECLIK which contains all the old
%    samples and the current samples from LIK.  This subfunction
%    is needed when using MCMC sampling (gp_mc).
%
%  See also
%    GP_MC

  if nargin == 2
    % Initialize the record
    reclik.type = 'Quantile-GP';

    % Initialize parameter
    reclik.sigma2 = [];

    % Set the function handles
    reclik.fh.pak = @lik_qgp_pak;
    reclik.fh.unpak = @lik_qgp_unpak;
    reclik.fh.lp = @lik_qgp_lp;
    reclik.fh.lpg = @lik_qgp_lpg;
    reclik.fh.ll = @lik_qgp_ll;
    reclik.fh.llg = @lik_qgp_llg;
    reclik.fh.llg2 = @lik_qgp_llg2;
    reclik.fh.llg3 = @lik_qgp_llg3;
    reclik.fh.tiltedMoments = @lik_qgp_tiltedMoments;
    reclik.fh.predy = @lik_qgp_predy;
    reclik.fh.invlink = @lik_qgp_invlink;
    reclik.fh.recappend = @lik_qgp_recappend;
    reclik.p=[];
    reclik.p.sigma2=[];
    if ~isempty(ri.p.sigma2)
      reclik.p.sigma2 = ri.p.sigma2;
    end
  else

    % Append to the record
    reclik.sigma2(ri,:)=lik.sigma2;
    if ~isempty(lik.p)
      reclik.p.sigma2 = lik.p.sigma2.fh.recappend(reclik.p.sigma2, ri, lik.p.sigma2);
    end
  end
end