text:
%SerialLink.MANIPLTY Manipulability measure
%
% M = R.maniplty(Q, OPTIONS) is the manipulability index (scalar) for the
% robot at the joint configuration Q (1xN) where N is the number of robot
% joints.  It indicates dexterity, that is, how isotropic the robot's
% motion is with respect to the 6 degrees of Cartesian motion. The measure
% is high when the manipulator is capable of equal motion in all directions
% and low when the manipulator is close to a singularity.
%
% If Q is a matrix (MxN) then M (Mx1) is a vector of  manipulability
% indices for each joint configuration specified by a row of Q.
%
% [M,CI] = R.maniplty(Q, OPTIONS) as above, but for the case of the Asada
% measure returns the Cartesian inertia matrix CI.
%
% R.maniplty(Q) displays the translational and rotational manipulability.
%
% Two measures can be computed:
% - Yoshikawa's manipulability measure is based on the shape of the velocity
%   ellipsoid and depends only on kinematic parameters (default).
% - Asada's manipulability measure is based on the shape of the acceleration
%   ellipsoid which in turn is a function of the Cartesian inertia matrix and
%   the dynamic parameters.  The scalar measure computed here is the ratio of
%   the smallest/largest ellipsoid axis.  Ideally the ellipsoid would be
%   spherical, giving a ratio of 1, but in practice will be less than 1.
%
% Options::
% 'trans'       manipulability for transational motion only (default)
% 'rot'         manipulability for rotational motion only
% 'all'         manipulability for all motions
% 'dof',D       D is a vector (1x6) with non-zero elements if the
%               corresponding DOF is to be included for manipulability
% 'yoshikawa'   use Yoshikawa algorithm (default)
% 'asada'       use Asada algorithm
%
% Notes::
% - The 'all' option includes rotational and translational dexterity, but
%   this involves adding different units.  It can be more useful to look at the
%   translational and rotational manipulability separately.
% - Examples in the RVC book (1st edition) can be replicated by using the 'all' option
%
% References::
%
% - Analysis and control of robot manipulators with redundancy,
%   T. Yoshikawa,
%   Robotics Research: The First International Symposium (M. Brady and R. Paul, eds.),
%   pp. 735-747, The MIT press, 1984.
% - A geometrical representation of manipulator dynamics and its application to
%   arm design,
%   H. Asada,
%   Journal of Dynamic Systems, Measurement, and Control,
%   vol. 105, p. 131, 1983.
% - Robotics, Vision & Control, P. Corke, Springer 2011.
%
% See also SerialLink.inertia, SerialLink.jacob0.

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

%TODO
% return the ellipsoid?

function [w,mx] = maniplty(robot, q, varargin)

    opt.method = {'yoshikawa', 'asada'};
    opt.axes = {'all', 'trans', 'rot'};
    opt.dof = [];

    opt = tb_optparse(opt, varargin);

    if nargout == 0
        opt.axes = 'trans';
        mt = maniplty(robot, q, 'setopt', opt);
        opt.axes = 'rot';
        mr = maniplty(robot, q, 'setopt', opt);
        for i=1:numrows(mt)
        fprintf('Manipulability: translation %g, rotation %g\n', mt(i), mr(i));
        end
        return;
    end

    if isempty(opt.dof)
        switch opt.axes
            case 'trans'
                dof = [1 1 1 0 0 0];
            case 'rot'
                dof = [0 0 0 1 1 1];
            case 'all'
                dof = [1 1 1 1 1 1];
        end
    else
        dof = opt.dof;
    end

    opt.dof = logical(dof);

    if strcmp(opt.method, 'yoshikawa')
        w = zeros(numrows(q),1);
        for i=1:numrows(q)
            w(i) = yoshi(robot, q(i,:), opt);
        end
    elseif strcmp(opt.method, 'asada')
        w = zeros(numrows(q),1);
        if nargout > 1
            dof = sum(opt.dof);
            MX = zeros(dof,dof,numrows(q));
            for i=1:numrows(q)
                [ww,mm] = asada(robot, q(i,:), opt);
                w(i) = ww;
                MX(:,:,i) = mm;
            end
        else
            for i=1:numrows(q)
                w(i) = asada(robot, q(i,:), opt);
            end
        end
    end

    if nargout > 1
        mx = MX;
    end

function m = yoshi(robot, q, opt)
    J = robot.jacob0(q);

    J = J(opt.dof,:);
    m2 = det(J * J');
    m2 = max(0, m2);    % clip it to positive
    m = sqrt(m2);

function [m, mx] = asada(robot, q, opt)
    J = robot.jacob0(q);

    if rank(J) < 6
        warning('robot is in degenerate configuration')
        m = 0;
        return;
    end

    Ji = pinv(J);
    M = robot.inertia(q);
    Mx = Ji' * M * Ji;
    d = find(opt.dof);
    Mx = Mx(d,d);
    e = eig(Mx);
    m = min(e) / max(e);

    if nargout > 1
        mx = Mx;
    end

author: petercorke

repo: robotics-toolbox-matlab


text:
hydro = struct();

hydro = readCAPYTAINE(hydro,'coer_comp_full.nc');
hydro = radiationIRF(hydro,10,[],[],[],[]);
hydro = radiationIRFSS(hydro,[],[]);
hydro = excitationIRF(hydro,20,[],[],[],[]);
writeBEMIOH5(hydro)
plotBEMIO(hydro)

author: WEC-Sim

repo: WEC-Sim

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


text:
function C = spm_dctmtx(N,K,n,f)
% Creates basis functions for Discrete Cosine Transform.
% FORMAT C = spm_dctmtx(N,K,n)
%     OR C = spm_dctmtx(N,K)
%     OR D = spm_dctmtx(N,K,n,'diff')
%     OR D = spm_dctmtx(N,K,'diff')
% N - dimension
% K - order
% n - optional points to sample
%____________________________________________________________________________
% spm_dctmtx creates a matrix for the first few basis functions of a one
% dimensional discrete cosine transform.
% With the 'diff' argument, spm_dctmtx produces the derivatives of the
% DCT.
%
% See:    Fundamentals of Digital Image Processing (p 150-154).
%         Anil K. Jain 1989.
%____________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id$

d = 0;

if nargin == 1, K = N; end;

if any(nargin == [1 2]),
    n = (0:(N-1))';
elseif nargin == 3,
    if strcmp(n,'diff'),
        d = 1;
        n = (0:(N-1))';
    elseif strcmp(n,'diff2'),
        d = 2;
        n = (0:(N-1))';
    else
        n = n(:);
    end
elseif nargin == 4,
    n = n(:);
    if strcmp(f,'diff'),
        d = 1;
    elseif strcmp(f,'diff2'),
        d = 2;
    else
        error('Incorrect Usage');
    end
else
    error('Incorrect Usage');
end

C = zeros(size(n,1),K);

if d == 0,
    C(:,1)=ones(size(n,1),1)/sqrt(N);
    for k=2:K
        C(:,k) = sqrt(2/N)*cos(pi*(2*n+1)*(k-1)/(2*N));
    end
elseif d == 1,
    for k=2:K
        C(:,k) = -2^(1/2)*(1/N)^(1/2)*sin(1/2*pi*(2*n*k-2*n+k-1)/N)*pi*(k-1)/N;
    end
elseif d == 2,
    for k=2:K,
        C(:,k) = -2^(1/2)*(1/N)^(1/2)*cos(1/2*pi*(2*n+1)*(k-1)/N)*pi^2*(k-1)^2/N^2;
    end;
else
    error('Can''t do this');
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [sessionInfo] = bz_getSessionInfo(basePath,varargin)
%[sessionInfo] = bz_getSessionInfo(basePath) loads the sessionInfo metadata
%for the recording in basePath. basePath should be in the format:
%       /whateverPath/baseName/
%           a file  basePath/baseName.sessionInfo.mat
%           or      basePath/baseName.xml
%           should exist.
%If no baseName.sessionInfo.mat exists, loads from the xml.
%
%INPUT
%   basePath            directory: '/whatevetPath/baseName/'
%   (options)
%       'saveMat'       (default: prompt)
%       'noPrompts'     (default: false) prevents prompts about
%                       saving/adding metadata
%       'editGUI'       (default: false) opens a GUI to edit select
%                       sessionInfo fields (beta, please add/improve!)
%
%OUTPUT
%   sessionInfo         metadata structure
%
%2017 DLevenstein and DTingley
%% inputs and defaults
p = inputParser;
addParameter(p,'noPrompts',false,@islogical);
addParameter(p,'editGUI',false,@islogical);
addParameter(p,'saveMat',false,@islogical);
parse(p,varargin{:})
noPrompts = p.Results.noPrompts;
editGUI = p.Results.editGUI;
saveMat = p.Results.saveMat;

if ~exist('basePath','var')
    basePath = pwd;
end
baseName = bz_BasenameFromBasepath(basePath);
filename = fullfile(basePath,[baseName,'.sessionInfo.mat']);

%% Load the stuff
%d = dir('*sessionInfo*'); %all files with sessioninfo in the name
if exist(filename,'file')
    sessionInfostruct = load(filename);
    %Checks that there is a single structure in the sessionInfo file
    varsInFile = fieldnames(sessionInfostruct);
    if numel(varsInFile)==1
        sessionInfo = sessionInfostruct.(varsInFile{1});
    else
        warning('Your .sessionInfo.mat has multiple variables/structures in it... wtf.')
        sessionInfo = sessionInfostruct;
    end
    SIexist = true;  %Marks that session info exists as expected
else
   warning(['could not find file ',baseName,'.sessionInfo.mat ',...
       'running LoadParameters instead..'])
   sessionInfo = LoadParameters(basePath);
   SIexist = false;
end

%% Check sessionInfo using bz_isSessionInfo and update if necessary
bz_isSessionInfo(sessionInfo);

%Here: prompt user to add any missing sessionInfo fields and save
if editGUI
    [ sessionInfo ] = bz_sessionInfoGUI(sessionInfo);
    SIexist = false;
% elseif ~isfield(sessionInfo,'region') && ~noPrompts
%     regionadd = questdlg(['Your sessionInfo is missing regions, ',...
%         'would you like to add them?'],'Add Regions?','Yes');
%     switch regionadd
%         case 'Yes'
%             [sessionInfo] = bz_sessionInfoGUI(sessionInfo,'Regions');
%             SIexist = false;
%         case 'Cancel'
%             return
%     end
end

%Should check that sessionInfo.session.name and sesioninfo.session.path
%match basePath....  if not, prompt the user to save with the correct
%information.

%% Save sessionInfo file
%Prompt user to save basePath/baseName.sessionInfo.mat
%if loaded from xml or changed
if ~noPrompts && ~SIexist %Inform the user that they should save a file for later
    savebutton = questdlg(['Would you like to save your sessionInfo in ',...
        filename '?'],'Save sessionInfo?','Yes');
    switch savebutton
        case 'Yes'
            saveMat = true;
        case 'No'
            saveMat = false;
        case 'Cancel'
            return
    end
end

if saveMat
    disp(['saving ',baseName,'.sessionInfo.mat'])
    save(filename,'sessionInfo');
end
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [candidates, priority] = run_rantalankila(im, num_candidates, varargin)
  if ~isdeployed
    old_path = path;
    addpath(fullfile('..', 'Rantalankila'));
    addpath(fullfile('..', 'Rantalankila', 'GCMex'));
    addpath(fullfile('..', 'Rantalankila', 'features'));
    addpath(fullfile('..', 'Rantalankila', 'vlfeat-0.9.18', 'toolbox', 'mex', 'mexa64'));

  end

  opts = get_opts();

  interp_data = load('/path/to/all/rantalankila_num_candidates_interpolation_data.mat');
  if num_candidates <= min(interp_data.num_candidates)
    opts.gc_branches = 0;
  elseif num_candidates >= max(interp_data.num_candidates)
    opts.gc_branches = 500;
  else
    opts.gc_branches = interp1(interp_data.num_candidates, interp_data.gc_branches, num_candidates);
    assert(~isnan(opts.gc_branches));
  end

  argl = length(varargin);
  if mod(argl,2) ~= 0 % string-value pairs
      error('Give pairs of extra arguments.');
  end
  for argi = 0:((argl/2) - 1)
      arg_str = varargin(2*argi + 1);
      arg_val = varargin(2*argi + 2);
      assert(isfield(opts, arg_str{1}));
      opts.(arg_str{1}) = arg_val{1};
  end

  opts.seg_method = 'felz';
  try
    [region_parts{1}, orig_sp{1}] = spagglom(im, opts);
  catch
    fprintf('rantalankila run on felsenzwalb segmentation failed\n');
  end
  opts.seg_method = 'slic';
  [region_parts{2}, orig_sp{2}] = spagglom(im, opts);

  candidates = sp_to_candidates(im, region_parts, orig_sp);
  priority = [];

  if ~isdeployed
    path(old_path);
  end
end

function bboxes = sp_to_candidates(im, region_parts, orig_sp)
%   candidates = [];
%   for sus = 1:length(orig_sp)
%     ys = double(orig_sp{sus}.pixels(:,1));
%     xs = double(orig_sp{sus}.pixels(:,2));
%     candidates = cat(1, candidates, [min(xs), min(ys), max(xs), max(ys)]);
%   end

  [h, w, ~] = size(im);

    % For each sp of orig_sp, find its bounding box
    for region_set = 1:length(region_parts)
        sp_edges{region_set} = zeros(length(orig_sp{region_set}), 4);
        bboxes{region_set} = zeros(length(region_parts{region_set}), 4);

        % find edges of Ri
        for ses = 1:length(orig_sp{region_set})
            Ri = false(h, w); % note ',', not '*' as below
            Ri(orig_sp{region_set}{ses}.spind) = 1;
            Xi = sum(Ri,1);
            Yi = sum(Ri,2);
            sp_edges{region_set}(ses,:) = [find(Xi,1,'first'), find(Yi,1,'first'), find(Xi,1,'last'), find(Yi,1,'last')];
        end

        % Using the above bounding boxes, solve bounding box of each region
        % proposal
        for j = 1:length(region_parts{region_set})
            sus = region_parts{region_set}{j};
            if ~isempty(sus)
                tg = sp_edges{region_set}(sus,:);
                bboxes{region_set}(j,:) = [min(tg(:,1)), min(tg(:,2)), max(tg(:,3)), max(tg(:,4))]; % bounding box is the most extreme values of individual bounding boxes
            else
              error('what is this?');
                bboxes = [];
            end
        end
    end % for each region set

    bboxes = cat(1, bboxes{:});
end

function [opts] = get_opts()
% Sets options for spagglom.m that are propagated through the functions. Many of the
% options are zero/one flags. The default settings were used to produce the
% results presented in the paper.

%% Superpixelation stage
opts.seg_method = 'slic'; % 'felz', 'slic' or 'qshift'. Two first recommended

% Felzenswalb superpixelation parameters
opts.felz_k = 50; % default 50
opts.felz_sigma = 0.8; % default 0.8
opts.felz_min_area = 150; % default 150

% SLIC superpixelation parameters
opts.slic_regularizer = 800; % default 800 (uses different definition than SLIC authors)
opts.slic_region_size = 20; % default 20

%% Features
opts.diagonal_connections = 0; % default 0. Whether only diagonally connected pixels are considered connected.
opts.dsift_step = 2; % default 2. calculate dsfit only every n step (quadratic speedup in n).

% Histogram features to use
opts.feature_dsift_bow = 1;    % 1 % denseSIFT bag-of-words
opts.feature_color_bow = 1;    % 2 % color bag-of-words in various color spaces
opts.feature_rgb_raw = 0;      % 3 % raw rgb histograms
opts.feature_grad_texture = 0; % 4 % feature used by van de Sande. Has two implementations, see code
opts.feature_lbp = 0;          % 6 % Local binary patterns
opts.feature_size = 1;         % 8 % size of combined superpixel
opts.features = 1:6; % This will be used in similarity.m to change histogram distances feature-wise
opts.features = opts.features(logical([opts.feature_dsift_bow, opts.feature_color_bow, opts.feature_rgb_raw, opts.feature_grad_texture, opts.feature_lbp ,opts.feature_size]));

opts.feature_weights = [1,1,2]; % default [1,1,2] with above dsift_bow, color_bow, size enabled and rest disabled

opts.collect_merged_regions = 1; % default 1. Every time a pair is merged during the greedy pairing algorithm, the new pair is saved as a region

opts.gc_branches = 15; % default 15. Number of graphcut branches.

opts.start_phase2 = 0.8; % 0.8 default. Score at which to change features and/or start branching

% Load precalculated data
opts.load_color_dict = 1; % default 1. 0 means dictionary will be created from the image (very slow). This option is used for both, rgb and lab features.
opts.load_dsift_dict = 1; % default 1. same as above but for dsift
opts.load_dsift_words = 0; % default 0. load precalculated dsift words using precalculated dict. This option overrides opts.load_dsift_dict.
opts.load_init_segs = 0; % default 0. Load initial felz or slic rgb segmentations with "default" parameters (see conditionals in code at spagglom.m)

if opts.load_dsift_dict
    load('/path/to/Rantalankila/dicts/dsift_dict_k500');
    opts.dsift_dict = dsift_dict;
    clear dsift_dict;
end
end

author: hosang

repo: detection-proposals


text:
function [x_snv] = snv(x)
% Standard Normal Variate
%
% [x_snv] = snv(x)
%
% input:
% x (samples x variables) data to preprocess
%
% output:
% x_snv (samples x variables) preprocessed data
%
% By Cleiton A. Nunes
% UFLA,MG,Brazil

[m,n]=size(x);
rmean=mean(x,2);
dr=x-repmat(rmean,1,n);
x_snv=dr./repmat(sqrt(sum(dr.^2,2)/(n-1)),1,n);

author: Sable

repo: mcbench-benchmarks


text:
function removal(seg, segnum, between, label, near, centroids, url, grad, texthist)

    disp 'Removal'
    im = imread(url);

    epsilon=0.01;
    l_label = label;

    between = zeros([segnum, segnum]);
    for i = 1:size(between, 1)
        for j = 1:size(between, 2)
            %between(i, j) = sum((grad(i,:)-grad(j,:)).^2./(grad(i,:)+grad(j,:)+epsilon)) + sum((texthist(i,:)-texthist(j,:)).^2./(texthist(i,:)+texthist(j,:)+epsilon));
            distance = [centroids(i,1), centroids(i,2); centroids(j,1), centroids(j,2)];
            distance = sqrt((distance(1,1)-distance(2,1))^2 + (distance(1,2)-distance(2,2))^2) / max(size(im, 1), size(im, 2));
            between(i, j) = sum(abs(grad(i,:) - grad(j,:))) + sum(abs(texthist(i,:) - texthist(j,:))) + distance;

            if i == j
                between(i,j) = 100;
            end
        end
    end

    nim = im;
    [gx gy] = gradient(double(seg));
    eim = (gx.^2+gy.^2)>1e-10;
    t = nim(:,:,1); t(eim)=0; nim(:,:,1)=t;
    t = nim(:,:,2); t(eim)=0; nim(:,:,2)=t;
    t = nim(:,:,3); t(eim)=0; nim(:,:,3)=t;
%      imshow(nim);
%      hold on;

    c_im = double(rgb2hsv(im));
    c1_im = c_im(:,:,1);
    c2_im = c_im(:,:,2);
    c3_im = c_im(:,:,3);
    lab = calHsvHist(c_im, seg, segnum);

    for i = 1:size(label,2)
        if label(i) == 0
            j = near(i);
            num = 0;
            while label(j) ~= 255
                [value, j] = min(between(i,:));
                between(i,j) = 100;
                num = num + 1;

            end
            near(i) = j;
            if num ~= 21
                plot([centroids(i,1) centroids(j,1)], [centroids(i,2) centroids(j,2)], 'r');
                %text(centroids(i,1),centroids(i,2),[ num2str(centroids(i,1)) ,  num2str(centroids(i,2))])
            end
        end
    end

    for i = 1:size(label, 2)
        if label(i) == 0 && label(near(i)) == 255
            j = near(i);
            temp3 = reshape(lab(j,101:150),[50,1]);
            c3_im(seg==i) = hist_match(c3_im(seg==i),temp3);
            temp2 = reshape(lab(j,51:100),[50,1]);
            c2_im(seg==i) = hist_match(c2_im(seg==i),temp2);
            temp1 = reshape(lab(j,1:50),[50,1]);
            c1_im(seg==i) = hist_match(c1_im(seg==i),temp1);
            %c2_im(seg==i) = c2_im(seg==i) + (median(c2_im(seg==j)) - median(c2_im(seg==i)));
            %c1_im(seg==i) = c1_im(seg==i) + (median(c1_im(seg==j)) - median(c1_im(seg==i)));
        end
    end
    c_im(:,:,1) = c1_im;
    c_im(:,:,2) = c2_im;
    c_im(:,:,3) = c3_im;
    %imshow(hsv2rgb(c_im));

%end

%%
test_im = hsv2rgb(c_im);
circle = zeros(size(c_im, 1), size(c_im, 2));
for i = 8:size(c_im, 1)-8
    for j = 8:size(c_im, 2)-8
        if label(seg(i,j)) ~= label(seg(i-1,j)) || label(seg(i,j)) ~= label(seg(i+1,j)) || label(seg(i,j)) ~= label(seg(i,j-1)) || label(seg(i,j)) ~= label(seg(i,j+1))
            circle(i-1:i+1,j-1:j+1) = 1;

        end
    end
end
h = fspecial('gaussian', 15, 15);
pattern = imfilter(test_im, h);

label = l_label;

for i = 1:segnum
    if label(i) == 0

        for ch = 1:3
            fig = test_im(:,:,ch);

            for x = 20:size(fig, 1) - 20
                for y = 20:size(fig,2) - 20

                    if seg(x,y) == i && seg(x-4, y) ~= i && circle(x,y) == 1
                        %if abs(test_im(x,y,3) - avg(i,3)) > 0.05
                            fig(x,y) = pattern(x,y,ch);
                            fig(x-4,y) = pattern(x,y,ch);
                        %end
                        %
                    elseif seg(x,y) == i && seg(x+4, y) ~= i && circle(x,y) == 1
                        %if abs(test_im(x,y,3) - avg(i,3)) > 0.05
                            fig(x,y) = pattern(x,y,ch);
                            fig(x+4,y) = pattern(x,y,ch);
                        %end
                        %
                    elseif seg(x,y) == i && seg(x, y+4) ~= i && circle(x,y) == 1
                        %if abs(test_im(x,y,3) - avg(i,3)) > 0.05
                            fig(x,y) = pattern(x,y,ch);
                            fig(x,y+4) = pattern(x,y,ch);
                        %end
                        %
                    elseif seg(x,y) == i && seg(x, y-4) ~= i && circle(x,y) == 1
                        %if abs(test_im(x,y,3) - avg(i,3)) > 0.05
                            fig(x,y) = pattern(x,y,ch);
                            fig(x,y-4) = pattern(x,y,ch);
                        %end
                        %
                    end
                end
            end
            test_im(:,:,ch) = fig;
        end
    end

end
    imshow(test_im);
end

author: kittenish

repo: Image-Shadow-Detection-and-Removal


text:
classdef TP4 < PROBLEM
% <multi> <real> <large/none> <robust>
% Test problem for robust multi-objective optimization
% delta --- 0.05 --- Maximum disturbance degree
% H     ---   50 --- Number of disturbances

%------------------------------- Reference --------------------------------
% A. Gaspar-Cunha, J. Ferreira, and G. Recio, Evolutionary robustness
% analysis for multi-objective optimization: benchmark problems, Structural
% and Multidisciplinary Optimization, 2014, 49: 771-793.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    properties
        delta;      % Maximum disturbance degree
        H;          % Number of disturbances
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            [obj.delta,obj.H] = obj.ParameterSet(0.05,50);
            obj.M = 2;
            if isempty(obj.D); obj.D = 10; end
            obj.lower    = zeros(1,obj.D);
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            PopObj(:,1) = (exp(PopDec(:,1))-1)/(exp(1)-1);
            g = 1 + 10*mean(PopDec(:,2:end),2);
            h = sin(4*pi*PopDec(:,1))/15 - PopDec(:,1) + 1;
            PopObj(:,2) = h.*g;
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(~,N)
            x = linspace(0,1,N)';
            R(:,1) = (exp(x)-1)/(exp(1)-1);
            R(:,2) = sin(4*pi*x)/15 - x + 1;
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            R = obj.GetOptimum(100);
        end
        %% Calculate the metric value
        function score = CalMetric(obj,metName,Population)
            switch metName
                case {'Mean_IGD','Mean_HV','Worst_IGD','Worst_HV'}
                    score = feval(metName,Population,obj);
                otherwise
                    score = feval(metName,Population,obj.optimum);
            end
        end
        %% Perturb solutions multiple times
        function PopX = Perturb(obj,PopDec,N)
            if nargin < 3; N = obj.H; end
            Delta = repmat(obj.delta.*(obj.upper-obj.lower),N*size(PopDec,1),1);
            w     = UniformPoint(N,obj.D,'Latin');
            Dec   = 2*Delta.*w(reshape(repmat(1:end,size(PopDec,1),1),1,[]),:) + repmat(PopDec,N,1) - Delta;
            Dec   = obj.CalDec(Dec);
            PopX  = SOLUTION(Dec,obj.CalObj(Dec),obj.CalCon(Dec));
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function data = read_data()

fprintf('Reading data ...\n');
fid = fopen('TrainingData.csv');
fgetl(fid);

data = zeros(37821,95);
days = {'"Saturday"','"Sunday"','"Monday"','"Tuesday"','"Wednesday"','"Thursday"','"Friday"'};
row_cnt = 0;

while ~feof(fid)
    row_cnt = row_cnt + 1;
    line = fgetl(fid);
    C = strread(line,'%s','delimiter',',');
    for i=1:95
        if i==5
            data(row_cnt,5) = find(strcmp(days,C{5}));
        else
            if strcmp(C{i},'NA')
                data(row_cnt,i) = -1000000;
            else
                data(row_cnt,i) = str2num(C{i});
            end
        end
    end
end

fprintf('Read in %d rows\n', row_cnt);

author: benhamner


text:
function outsig = noise(siglen,varargin)
% NOISE  Stochastic noise generator
%   Usage: outsig = noise(siglen,nsigs,type);
%
%   Input parameters:
%       siglen    : Length of the noise (samples)
%       nsigs     : Number of signals (default is 1)
%       type      : type of noise. See below.
%
%   Output parameters:
%       outsig    : $siglen \times nsigs$ signal vector
%
%   `noise(siglen,nsigs)` generates *nsigs* channels containing white noise
%   of the given type with the length of *siglen*. The signals are arranged as
%   columns in the output. If only *siglen* is given, a column vector is
%   returned.
%
%   `noise` takes the following optional parameters:
%
%     'white'  Generate white (gaussian) noise. This is the default.
%
%     'pink'   Generate pink noise.
%
%     'brown'  Generate brown noise.
%
%     'red'    This is the same as brown noise.
%
%   By default, the noise is normalized to have a unit energy, but this can
%   be changed by passing a flag to |setnorm|.
%
%   Examples:
%   ---------
%
%   White noise in the time-frequency domain:::
%
%     sgram(noise(5000,'white'),'dynrange',70);
%
%   Pink noise in the time-frequency domain:::
%
%     sgram(noise(5000,'pink'),'dynrange',70);
%
%   Brown/red noise in the time-frequency domain:::
%
%     sgram(noise(5000,'brown'),'dynrange',70);
%
%   See also: setnorm

%   AUTHOR: Hagen Wierstorf and Peter L. Søndergaard.

% ------ Checking of input parameter -------------------------------------

if nargin<1
  error('%s: Too few input parameters.',upper(mfilename));
end;

if ~isnumeric(siglen) || ~isscalar(siglen) || siglen<=0
    error('%s: siglen has to be a positive scalar.',upper(mfilename));
end

definput.import={'setnorm'};
definput.importdefaults={'2'};
definput.flags.noisetype={'white','pink','brown','red'};
definput.keyvals.nsigs=1;
[flags,kv,nsigs]  = ltfatarghelper({'nsigs'},definput,varargin);

if flags.do_white
  outsig=randn(siglen,nsigs);
end;

if flags.do_brown || flags.do_red
  outsig=cumsum(randn(siglen,nsigs));
end;

if flags.do_pink
  % --- Handle trivial condition

  if siglen==1
    outsig=ones(1,nsigs);
    return;
  end;

  % ------ Computation -----------------------------------------------------
  fmax = floor(siglen/2)-1;
  f = (2:(fmax+1)).';
  % 1/f amplitude factor
  a = 1./sqrt(f);
  % Random phase
  p = randn(fmax,nsigs) + i*randn(fmax,nsigs);
  sig = bsxfun(@times,a,p);

  outsig = ifftreal([ones(1,nsigs); sig; 1/(fmax+2)*ones(1,nsigs)],siglen);

end;

outsig=setnorm(outsig,flags.norm);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function varAtDistX = LVx(structNum,varCalcDist,planC)
% function varAtDistX = LVx(structNum,varCalcDis,planC)
%
% This function calculates the variance of absolute difference between
% voxels within structure (structNum) that are x cm (varCalcDist) apart.
%
% APA, 10/03/2013

indexS = planC{end};

% Get scan associated with this structNum
scanNum = getStructureAssociatedScan(structNum, planC);

% Get x,y,z vals for this structure voxels
[rasterSegments, planC, isError] = getRasterSegments(structNum, planC);
[mask3M, uniqueSlices] = rasterToMask(rasterSegments,scanNum,planC);
[i,j,k] = find3d(mask3M);
[xValsScan, yValsScan, zValsScan] = getScanXYZVals(planC{indexS.scan}(scanNum));
structXvals = xValsScan(j);
structYvals = yValsScan(i);
structZvals = zValsScan(uniqueSlices(k));
scanStructV = getScanAt(scanNum, structXvals, structYvals, structZvals, planC);

% Randomly sample 1000 points
numSamples = 1000;
randIndV = randi(length(structXvals),numSamples,1);
samplePtsXv = structXvals(randIndV);
samplePtsYv = structYvals(randIndV);
samplePtsZv = structZvals(randIndV);
sampleScanV = getScanAt(scanNum, samplePtsXv, samplePtsYv, samplePtsZv, planC);

% calculate distance between sample points and the structure voxels
distM = sepsq([samplePtsXv(:) samplePtsYv(:) samplePtsZv(:)]', [structXvals(:) structYvals(:) structZvals(:)]');
distM = distM.^0.5;

% Get scan resolution
dx = abs(xValsScan(1) - xValsScan(2));
dy = abs(yValsScan(1) - yValsScan(2));
dz = abs(zValsScan(1) - zValsScan(2));
res = max([dx dy dz])/2;
%indM = distM >= varCalcDist - res & distM <= varCalcDist + res;
indM = distM <= varCalcDist + res;

if ~any(indM(:))
    varAtDistX = NaN;
    warning(['Cannot sample points that are at ',num2str(varCalcDist),' cm'])
    return;
end

sampleScanM = repmat(sampleScanV',[1 length(scanStructV)]);
scanStructM = repmat(scanStructV,[length(sampleScanV) 1]);
absScanDiffAllPtsM = abs(scanStructM - sampleScanM);
absScanDiffSampledPtsV = absScanDiffAllPtsM(indM);

% Calculate variance for the samples
varAtDistX = var(absScanDiffSampledPtsV);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


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

%%
clear; close all; clc;

%%
% PLOT SETTINGS
fontSize=15;
faceAlpha1=0.5;

%% Creating a triangulated mesh of a box

%%
% Specifying dimensions and number of elements for each direction
boxDim=[4 5 6]; %Width in each direction
pointSpacing=1; %Desired point spacing

%%
% Using |triBox| to build the patch model

[F,V,faceBoundaryMarker]=triBox(boxDim,pointSpacing);

%%
% Visualisation
cFigure; hold on;
title('Box triangular faces and normals','FontSize',fontSize);

gpatch(F,V,faceBoundaryMarker,'k',1);
patchNormPlot(F,V);

axisGeom(gca,fontSize);
colormap(gjet(6)); icolorbar;
drawnow;

%%
%
% <<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:
%% Mesh Generator for Tecplot
% Generating a 2D mesh for Tecplot for an axisymmetric diffuser.
% by Manuel Diaz, 2012.09.15.

clear all; close all;
%% Parameters
L0 = 0; L1 = 10 + L0; L2 = 6 + L1 + L0;
H0 = 0; H1 = 0.5;
theta = deg2rad(5); % theta = 5 degress converted to radians
H2 = H1 + L2*tan(theta);

%% Number of elements of edges (ne)
ne_x1 = 25; ne_x2 = 40; ne_y = 40;

%% End Points for the rectangular zones:
Z1=[ L0,H1 ; L1,H1 ; L1,H0 ; L0,H0 ]; % Zone1
Z2=[ L1,H1 ; L2,H2 ; L2,H0 ; L1,H0 ]; % Zone2

%% Zone 1:
% There are 4 edge lines, each having several edge points.
% We use the Function "Pinterpol" to create the edge points for zone 1.

E1 = Pinterpol(Z1(1,1),Z1(1,2),Z1(2,1),Z1(2,2),ne_x1);
E2 = Pinterpol(Z1(4,1),Z1(4,2),Z1(3,1),Z1(3,2),ne_x1);

% Now generating the mesh for Zone 1
M1 = zeros(ne_x1+1,ne_y+1,2);
for i = 1:(ne_x1+1);
    V = Pinterpol(E1(i,1),E1(i,2),E2(i,1),E2(i,2),ne_y);
    M1(i,1:(ne_y+1),1:2) = V; % mesh1
end

%% Zone 2:
% Similarly for zone 2
E1 = Pinterpol(Z2(1,1),Z2(1,2),Z2(2,1),Z2(2,2),ne_x2);
E2 = Pinterpol(Z2(4,1),Z2(4,2),Z2(3,1),Z2(3,2),ne_x2);

% generating the mesh for Zone 2
M2 = zeros(ne_x2+1,ne_y+1,2);
for i = 1:(ne_x2+1);
    V = Pinterpol(E1(i,1),E1(i,2),E2(i,1),E2(i,2),ne_y);
    M2(i,1:(ne_y+1),1:2) = V; % mesh2
end

%% To Tecplot
% Writing Output data for Tecplot format:

file = fopen('mesh1.tec','w');
% h1 gets the handel for the file "mesh1.tec".
% 'w' specifies that it will be written.
% similarly 'r' is for reading and 'a' for appending.

fprintf(file, 'TITLE = "Mesh from Matlab"\n');
fprintf(file, 'VARIABLES = "X" "Y"\n');
fprintf(file, 'ZONE T = "Inlet zone"\n');
fprintf(file, 'I = %d, J = %d, K = 1, F = POINT\n\n', ne_x1+1,ne_y+1);

for j = 1:(ne_y+1)
    for i = 1:(ne_x1+1)
        fprintf(file, '%f\t%f\n', M1(i,j,1),M1(i,j,2));
    end
end

fprintf(file, 'ZONE T="Diffuser zone"\n');
fprintf(file, 'I = %d, J=%d, k=1, F=POINT\n\n', ne_x2+1,ne_y+1);

for j = 1:(ne_y+1)
    for i = 1:(ne_x2+1)
        fprintf(file, '%f\t%f\n', M2(i,j,1),M2(i,j,2));
    end
end

fclose(file);

repo: Aero-matlab


text:
function Stock_PCA()
clc; clear all; clearvars
addpath('functions_IO');

%
% Download all necessary data
%
finviz_dir = '/Users/ansonwong/Desktop/MLF_data/finviz'; % finviz file will be finviz.csv
chart_dir = '/Users/ansonwong/Desktop/MLF_data/charts'; % charts will be XXXX.txt
%datestr_use = datestr(now,'mmm_dd_yyyy');

need_download_finviz = 0; % downloads finviz.csv and places it in finviz_dir
need_download_chart = 0; % downloads all charts into chart_dir based on finviz.csv
need_convert_finviz = 0; % create a finviz.mat that is based off of finviz.csv data (numbers/NaN only)

% Download finviz and chart data (as well as converting finviz data to a
% convenient format)
if(need_download_finviz == 1)
  download_finviz_data(finviz_dir);
end
if(need_download_chart == 1)
  download_chart_data(chart_dir,finviz_dir)
end
if(need_convert_finviz == 1)
  data2num_finviz(finviz_dir)
end

%
% Analyze data now
%

% Load all data
% Now we have 'FHEAD', 'FHEAD_names', 'name_cell', 'data_label', 'data'
finviz2_filename = strcat(finviz_dir,'/finviz.mat');
M = load(finviz2_filename,'-mat');

% Set new data variables
X = M.data'; % n=7091 x d=62
XS = M.name_cell; % n=7091 x 4
FH = M.FHEAD; % structure array with indices for X
FH_name = M.FHEAD_names; % cell array of string n=7091 x 1
clear M
[n,d]=size(X);

%{
    'MarketCap'
    'PE'
    'PEForward'
    'PEG'
    'PS'
    'PB'
    'PCash'
    'PFreeCash'
    'DividendYield'
    'PayoutRatio'
    'EPS'
    'EPSGrowthThisYear'
    'EPSGrowthNextYear'
    'EPSGrowthPast5Years'
    'EPSGrowthNext5Years'
    'SalesGrowth5Years'
    'EPSGrowthLastQuarter'
    'SalesGrowthLastQuarter'
    'SharesOutstanding'
    'Float'
    'InsiderOwnership'
    'InsiderTransactions'
    'InstitutionalOwnership'
    'InstitutionalTransactions'
    'FloatShort'
    'ShortRatio'
    'ReturnOnAssets'
    'ReturnOnEquity'
    'ReturnOnInvestment'
    'CurrentRatio'
    'QuickRatio'
    'LTDebtEquity'
    'TotalDebtEquity'
    'GrossMargin'
    'OperatingMargin'
    'ProfitMargin'
    'PerformanceWeek'
    'PerformanceMonth'
    'PerformanceQuarter'
    'PerformanceHalfYear'
    'PerformanceYear'
    'PerformanceYTD'
    'Beta'
    'AverageTrueRange'
    'VolatilityWeek'
    'VolatilityMonth'
    'SMA20day'
    'SMA50day'
    'SMA200day'
    'High50day'
    'Low50day'
    'High52week'
    'Low52week'
    'RSI14'
    'ChangeFromOpen'
    'Gap'
    'AnalystRecommendation'
    'AverageVolume'
    'RelativeVolume'
    'Price'
    'ChangePercentage'
    'Volume'
    'TargetPrice'
%}

%{
corrcoef(X(:,FH.ReturnOnAssets),X(:,FH.ReturnOnAssets))

%}
% Find the NaN value fractions
if(0)
  fprintf('[n=%d,d=%d]\n',n,d);
  FID = fopen('stats_pca.txt','w');
  for j=1:d
    n_nonnan = sum(~isnan(X(:,j)));
    fprintf('%d) %s (%.0f%% good)\n',j,FH_name{j},100*n_nonnan/n);
    fprintf(FID,'%d) %s (%.0f%% good)\n',j,FH_name{j},100*n_nonnan/n);
  end
  fclose(FID);
end

%
% Find a subset of X
%
%index_Xsubset = find(X(:,FH.MarketCap)>100);
index_Xsubset = find(X(:,FH.AverageVolume)>1000);
fprintf('index_subset = %d\n',length(index_Xsubset));
X = X(index_Xsubset,:);

%plot_2var(1,FH.ReturnOnAssets,FH.ChangePercentage,X,FH_name);

% Plot between 2 variables
% 1) Gap, ChangePercentage
% 2)
%plot_2var(1,FH.ReturnOnAssets,FH.ChangePercentage,X,FH_name);
%plot_2var(1,FH.ReturnOnAssets,FH.ChangePercentage,X,FH_name);
%plot_2var(1,FH.ReturnOnEquity,FH.ChangePercentage,X,FH_name);
%plot_2var(1,FH.ReturnOnEquity,FH.ChangePercentage,X,FH_name);

%
% Use features and reduce to non-NaN, X_nonnan
%
j_list = ...
  [FH.FloatShort,...
   FH.ChangePercentage];

size(X)
[X,indices_nan] = collect_nonnan(X,j_list);
size(X)

rho = corr(X,'type','Pearson')
figure(2)
clf(2)
plot(X(:,1),X(:,2),'.');

%[X,indices_nan] = collect_nonnan(X,1:d);

%
% Normalize Xtilde_nonnan = X_nonnan
%
X = standardizeCols(X);

%
% Perform PCA
%
[coeff,score,latent,tsquared,explained] = pca(X);

coeff

for i=1:length(explained)
  fprintf('PC %d = %.1f%%\n',i,explained(i));
end

%{
indices_nan = zeros(n,1);
for j=length(j_list),
  indices_nan = indices_nan + isnan(X(:,j_list(j)));
end
indices_nan = indices_nan > 0;
X(indices_nan,:) = []; % remove tickers without marketcap info

size(X)
sum(~isnan(X))
%}

clearvars
rmpath('functions_IO');
end

% Given a j_list, reduces the list so that we have no nans
function [X_new,indices_nan] = collect_nonnan(X,j_list)
[n,d]=size(X);

% Find ticker indices that are nan in j_list
indices_nan = zeros(n,1);
for j=1:length(j_list),
  indices_nan = indices_nan + isnan(X(:,j_list(j)));
end
indices_nan = indices_nan > 0;
indices_nonnan = ~indices_nan;

% Make new X
X_new = X(indices_nonnan,:);
vec_setempty = 1:d;
vec_setempty(j_list) = [];
X_new(:,vec_setempty) = [];

% Make sure every value of X_new is non-nan
if( any( diff(sum(~isnan(X_new))) ) )
  error('collect_nonnan did not work')
end

end

function plot_2var(fignum,jx,jy,X,FH_name)
index_nonnan = find(~isnan(X(:,jx)) & ~isnan(X(:,jy)));
xplot = X(index_nonnan,jx); yplot = X(index_nonnan,jy);
figure(fignum)
clf(fignum)
plot(xplot,yplot,'.')
hold on
xlabel(FH_name{jx});
ylabel(FH_name{jy});
hold off
end

author: gudbrandtandberg


text:
% Copyright 2019 Jonas Koenemann, Moritz Diehl, University of Freiburg
% Redistribution is permitted under the 3-Clause BSD License terms. Please
% ensure the above copyright notice is visible in any derived work.
%

function [nFails] = core(testExamples,saveLog,suffix)
  % runTests()
  % runTests(testExamples)
  % runTests(1,1,'computername')
  ocl.utils.checkStartup()

  if nargin < 1
    testExamples = false;
  end
  if nargin <=1
    saveLog = false;
    suffix = '';
  else
    suffix = ['_',suffix];
  end

  testDir = getenv('OPENOCL_TEST');
  oclDir = getenv('OPENOCL_PATH');

  % go to main dir and get current git hash
  cd(oclDir);
  [status,version] = system('git rev-parse HEAD');
  if status~=0
    ocl.utils.info('Could not get git hash, log filename will not include the hash. This is not a problem.')
    version='nogithash';
  end

  if isempty(testDir)
    ocl.utils.error('Test directory not set. Run StartupOCL again.')
  end

  ocl.utils.setTestRun(true)

  tests{1}.name = 'Variable';         tests{end}.file = 'ocl.tests.variable';
  tests{end+1}.name = 'TreeVariable'; tests{end}.file = 'ocl.tests.treeVariable';
  tests{end+1}.name = 'VarStructure'; tests{end}.file = 'ocl.tests.structure';
  tests{end+1}.name = 'ocl.Model';    tests{end}.file = 'ocl.tests.model';
  tests{end+1}.name = 'ocl.Stage';       tests{end}.file = 'ocl.tests.stage';
  tests{end+1}.name = 'ocl.Integrator';   tests{end}.file = 'ocl.tests.integrator';

  if testExamples
    tests{end+1}.name = 'Examples';      tests{end}.file = 'ocl.tests.examples';
    tests{end+1}.name = 'Features';      tests{end}.file = 'ocl.tests.features';
  end

  NTests = length(tests);

  % turn off figures for testing examples
  if testExamples
    close all
    set(0,'DefaultFigureVisible','off');
  end

  %% run all tests
  testResults = cell(1,NTests);
  for k=1:NTests
    test = tests{k};
    testResults{k} = runTest(test.name,str2func(test.file));
  end

  nFails = 0;
  for k=1:NTests
    if testResults{k}.passed == false
      nFails = nFails + 1;
    end
  end

  %% save results
  fileName = [datestr(now,'yyyy-mm-dd_HHMMSS'), '_', version(1:7), suffix,  '.txt'];
  filePath = fullfile(testDir,fileName);
  resultsFile = fopen(filePath,'w');
  fprintf(resultsFile,'Test on %s\nVersion: %s\n\n',datestr(now),version);

  for k=1:NTests
    printResults(testResults{k});
  end

  if nFails >0
    outputString = sprintf('%i Tests failed.\n\n',nFails);
    fprintf(resultsFile,outputString);fprintf(2,outputString);
  else
    outputString = 'All tests passed!\n\n';
    fprintf(resultsFile,outputString);fprintf(outputString);
  end

  fclose(resultsFile);

  if saveLog
    copyfile(filePath,'Log');
  end

  if testExamples
    set(0,'DefaultFigureVisible','on');
  end

  ocl.utils.setTestRun(false);

  function testResult = runTest(testName,scriptHandle)
    testResult = struct;
    testResult.name = testName;
    testResult.passed = true;
    testResult.runtime = 0;
    testResult.exception = '';
    try
      testTic = tic;
      testResult.outputs = cell(nargout(scriptHandle),1);
      [testResult.outputs{:}] = scriptHandle();
      testResult.runtime = toc(testTic);
    catch exception
      testResult.passed = false;
      testResult.exception = exception;
    end
  end

  function printResults(testResult)
    if testResult.passed
    	outputString = sprintf('%s Tests passed\n',testResult.name);
      fprintf(resultsFile,outputString);fprintf(outputString);
      outputString = sprintf('It took %.4f seconds.\n\n',testResult.runtime);
      fprintf(resultsFile,outputString);fprintf(outputString);

      for i=1:length(testResult.outputs)
        out = testResult.outputs{i};
        names = fieldnames(out);
        sum = 0;
        for j=1:length(names)
          val = out.(names{j});
          sum = sum+val;
          outputString = sprintf('Test %i, %s: %.4f seconds.\n', i, names{j}, val);
          fprintf(resultsFile,outputString);fprintf(outputString);
        end
        outputString = sprintf('Test %i, sum: %.4f seconds.\n\n', i, sum);
        fprintf(resultsFile,outputString);fprintf(outputString);
      end

    else
      outputString = sprintf('%s Tests failed\n',testResult.name);
      fprintf(resultsFile,outputString);fprintf(outputString);

      if ocl.utils.isOctave()
        disp(getReportOctave(testResult.exception))  %% OCTAVE
      else
        disp(getReport(testResult.exception))  %% MATLAB
      end
    end
  end
end

author: OpenOCL

repo: OpenOCL

save_path: github-repos/MATLAB/OpenOCL-OpenOCL


text:
function res = backtest(rex,Nh);
% PURPOSE: Backtesting of temporal disaggregation.
% -----------------------------------------------------------------------
% SYNTAX: res = backtest(rex,Nh);
% -----------------------------------------------------------------------
% OUTPUT: res: a structure with...
%         rex  : reference structure
%         yh   : recursive estimations from FROM T=N+Nh T=N
%         yrev : revision as %, compared with final (full info) estimation
% -----------------------------------------------------------------------
% INPUT: rex: a structure generated by chowlin(), fernandez(), litterman()
%             or ssc()
%        Nh: 1x1 --> number of low-freq. to compute revisions
% -----------------------------------------------------------------------
% LIBRARY: chowlin, fernandez, litterman, ssc
% -----------------------------------------------------------------------
% NOTE: Revision of estimates when adding low-frequency data from T=N-Nh to
% T=N

% written by:
%  Enrique M. Quilis
%  Macroeconomic Research Department
%  Ministry of Economy and Competitiveness
%  <enrique.quilis@mineco.es>

% Version 2.1 [May 2008]

t0 = clock;

% ------------------------------------------------------------
% Loading the structure

Y 		= rex.Y;
x 		= rex.x(:,2:end); %Excluding intercept
ta 	    = rex.ta;
sc 	    = rex.sc;
type 	= rex.type;
opC     = rex.opC;

% ------------------------------------------------------------
% Size of the problem

N = rex.N;     % Size of low-frequency input
n = rex.n;     % Size of high-frequency input
p = rex.p - 1; % Number of regressors (without intercept)

% ------------------------------------------------------------
% BACKTEST: REVISIONS OF ESTIMATES, AS Y EXPANDS (x all sample)
% ------------------------------------------------------------

x = x(1:sc*N,:);
yh = NaN * ones(n,Nh+1);
for h=Nh:-1:0
    Yh = Y(1:end-h);
    % Calling temporal disaggregation function
    switch rex.meth
        case {'Chow-Lin'}
            rl   = rex.rl;
            res1 = chowlin(Yh,x,ta,sc,type,opC,rl);
        case {'Fernandez'}
            res1 = fernandez(Yh,x,ta,sc,opC);
        case {'Litterman'}
            rl   = rex.rl;
            res1 = litterman(Yh,x,ta,sc,type,opC,rl);
        case {'Santos Silva-Cardoso'}
            rl   = rex.rl;
            res1 = ssc(Yh,x,ta,sc,type,opC,rl);
        otherwise
            error ('*** BACKTESTING IS NOT AVAILABLE FOR THIS METHOD ***');
    end
    % Generating yh
    yh(:,Nh-h+1) = res1.y;
end

% Computation of revisons as percentages
yrev = (yh - kron(ones(1,Nh+1),yh(:,end))) ./ yh;

% -----------------------------------------------------------
% Loading the structure
% -----------------------------------------------------------

% Reference structure
res.rex = rex;
% Sequence of estimates
res.yh = yh;
% Revisions of stimates
res.yrev = yrev;
% Elapsed time
res.et = etime(clock,t0);


text:
function [evtonsets,X,bestp,meff,origvec] = ga2(gensize,numgen,conditions,mspec,varargin)
%[evtonsets,X,bestp,meff,origvec] = ga2(gensize,numgen,conditions,mspec,[nooverlap])
%
% everything in conditions should be in number of frames (TRs), not seconds
% model_setup_ui.m automatically creates the correct format
%
% optional last argument is a 1 or 0.  1 Constrains the events so that
% events of different trials do not overlap.
%
% recommended: test conditions setup before running with
% [X,pv] = construct_model(mspec,conditions,[]);
% also [X,pv] = paramvec2onsets(mspec,conditions,[]);
%
% OR
%
% [X,pv] = construct_model(mspec,conditions,[],[1 14]);
% ...to constrain overlaps among trials to a min of 14 TRs apart,
% (variable trial onsets will only be properly constrained by the 2nd element).
% OR add in mspec.toverlap = [1 10], e.g.
%
% Tor Wager
% April 19, 2002
% Modified April 2, 2003: changed xover, convergence, constrained
%   so that min time between trial onsets is maintained
%   also added non-identity intrinsic autocorrelation
%   and capability for d-optimality
%
% Modified May 14, 2003: changed autocorrelation function to colored noise
%
% convergence criteria:
%   all efficiency scores are the same, OR no improvement in 100 generations
%
% defaults:
%   exponential autocorrelation function based on empirical estimates
%   saturation (clipping) of the predicted response at 2.5 times the original height
%       in order to approximate nonlinear effects
%   A-optimal design efficiency computation
%
%   Also possible:
%   Combo A-D optimal design, equally weighted
%       for A-optimal design, mspec.dopt = 0;
%       for D-optimality, specify mspec.dopt = 1;
%   But I'm still testing the D-optimality stuff.
%
% to see a test of A vs D optimality, try:
% Vi = eye(mspec.numframes);
% for i = 1:1000,
%    [X,paramvec] = construct_model(mspec,conditions,[]);
%    a(i) = calcefficiency(ones(1,size(X,2) - 1),[],pinv(X),Vi,0);
%    d(i) = 1 ./ (inv(det(X'*X)) .^ (1./size(X,2)));
% end
% d = d - mean(d);
% figure; plot(a,d,'x'); xlabel('A-efficiency'),ylabel('D-efficiency')
% a1 = calcefficiency(ones(1,size(X,2) - 1),[],pinv(X),Vi,0);
% d1 = 1 ./ (inv(det(X'*X)) .^ (1./size(X,2)));
% hold on;plot(a1,d1,'rs','MarkerFaceColor','r')

t0 = clock;

% --------------------------------------------------------------------
%
% * Setup and defaults
%
% --------------------------------------------------------------------

if length(varargin) > 0, nooverlap = varargin{1};,else, nooverlap = 0;,end

for j = 1:gensize, newvec{j} = [];,end

% intrinsic autocorrelation
% --------------------------------------------------------------------
%Vi = eye(mspec.numframes);      % use this for scanner white noise

% new autocorr, from VNL experiment (n = 10)
% this is invariant to the TR as well.
f = inline('5.6034 * exp(-.93 * x)','x');
t1 = f(1:mspec.TR:10); t1 = t1 ./ t1(1);

Vi = getv('make',t1,ceil(mspec.numframes));

% nonlinearity and fitness function
% --------------------------------------------------------------------

% nonlinear clipping threshold
satthresh = 2.5;

% display results vs. random?
showres = 1;

% converge if no change in this many generations
genconverge = 10;

% optimality - default is A
if ~isfield(mspec,'dopt'), mspec.dopt = 0;, end
if mspec.dopt == 2
    disp('Computing both A- and D-optimality'),
elseif mspec.dopt == 1
    disp('Computing D-optimality'),
else
    disp('Computing A-optimality'),
end

if nooverlap & length(nooverlap) == 1
    % we need a list of all trial onset conditions to make sure we don't create
    % illegal lists in xover
    % older, more cumbersome way based on each condition
    trialonsetinfo = get_tonsetinfo(conditions);
elseif nooverlap & length(nooverlap) > 1
    % newer way: minimum time constraint between all trial onsets
    trialonsetinfo = nooverlap(2);
else
    % no constraints
    trialonsetinfo = -Inf;
end

% model seeding - m-sequences, 10 %
% works for designs with one subpart per trial - don't know about
% others.  will have to insert more paramvec entries
% --------------------------------------------------------------------
nmseq = round(gensize ./ 10);
%[mseqs,newvec(1:nmseq)] = genMseq(mspec.numframes,[],[],nmseq,length(conditions));

[tmp,mseqs] = genMseq(mspec.numframes,[],[],nmseq,length(conditions));  % 10 sequences

for i = 1:length(mseqs)
    % get paramvec, which constrains length to correct values
    [tmp,pvm] = construct_model(mspec,conditions,mseqs{i});
    newvec{i} = pvm;
end

% set up contrasts across conditions
% --------------------------------------------------------------------
    numframes = [];,    % ne4eded for ga3power at the end, and for contrasts
    for i = 1:length(conditions),
        for j=1:length(conditions(i).subcond), numframes(end+1) = conditions(i).subcond(j).hrfest;,
        end
    end

if isfield(mspec,'contrasts')
    % executes the code to apply contrast, apply_dx_contrast.m, in
    % construct_model
    disp(['Found contrast to collapse columns: ']), mspec.contrasts
end

f1 = figure; set(gcf,'Color','w'); hold on; title('Max efficiency'),xlabel('Generation')

% matrix of multiplier values for each param in first cell to add noise
tmp = cat(1,newvec{:}); for i = 1:length(tmp), tmp{i} = diff(tmp{i}); end , tmp = cat(2, tmp{:});
noisestd = .2 .* range(tmp);
disp(['Noise std if no change is: ' num2str(noisestd)])

% --------------------------------------------------------------------
%
% * Start generations
%
% --------------------------------------------------------------------

for i = 1:numgen

    t1 = clock;

    for j = 1:gensize

        % --------------------------------------------------------------------
        % * make models
        % --------------------------------------------------------------------
        % if empty (1st generation), then it constructs a model.
        [X,paramvec{j}] = construct_model(mspec,conditions,newvec{j},nooverlap);
        X = modelSaturation(X,satthresh);

        % tmp for Emily
        %X(:,15*4+1:15*5) = X(:,15*4+1:15*5) + X(:,15*6+1:15*7); X(:,15*6+1:15*7) = [];

        XX{j} = X;

        if i == i & j == 1
            cweights = ones(1,size(X,2) - 1);
        end

        if mspec.dopt == 0
            eff(j) = calcEfficiency(cweights,[],pinv(X),Vi,0);
        elseif mspec.dopt == 1
            eff(j) = det(X'*X) .^ (1./size(X,2)); % 1 ./ (inv(det(X'*X)) .^ (1./size(X,2)));  % D-optimality
        elseif mspec.dopt == 2  % for A- D-optimality
            eff2(j) = calcEfficiency(cweights,[],pinv(X),Vi,0);
            eff(j) = det(X'*X) .^ (1./size(X,2));  % D-optimality
        end

    end

    if mspec.dopt == 2  % for A- D-optimality
        zeff = (eff - mean(eff)) ./ std(eff);
        zeff2 = (eff2 - mean(eff2)) ./ std(eff2);
        zeff = zeff + zeff2;
        meffi = max(zeff);
        meffind = find(zeff == meffi); meffind = meffind(1);
        meff(i) = eff(meffind);
        meff2(i) = eff2(meffind);
        t2 = clock;
        fprintf(1,'Gen %3.0f: D-eff is %3.4f, A-eff is %3.4f, model %3.0f, time = %3.2f, elapsed = %3.0f\n',i,meff(i),meff2(i),meffind,etime(t2,t1),etime(t2,t0));
        figure(f1); subplot 121; plot(1:i,meff,'r','LineWidth',2); title('D-efficiency')
        subplot 122; plot(1:i,meff2,'b','LineWidth',2); title('A-efficiency')
        drawnow

    else
        [meffi, meffind] = max(eff);
        meffind = meffind(1);
        meff(i) = meffi;
        t2 = clock;
        fprintf(1,'Gen %3.0f: best eff is %3.4f, model %3.0f, time = %3.2f, elapsed = %3.0f\n',i,meffi,meffind,etime(t2,t1),etime(t2,t0));
        figure(f1); plot(1:i,meff,'LineWidth',2); drawnow
    end

    % --------------------------------------------------------------------
    % add noise if no improvement at all
    % rounding and illegal values should be handled by construct_model
    % but need to be integers >= 0 for HRF power below
    % in next iteration
    % --------------------------------------------------------------------
    if meffind == 1

        for j = 1:length(paramvec)
            if j == meffind
                % leave the best one alone!
            else
                for k = 1:length(paramvec{j})
                    paramvec{j}{k} = paramvec{j}{k} + round( noisestd .* randn(size(paramvec{j}{k})) );

                    paramvec{j}{k}(paramvec{j}{k} < 0) = 0;
                end
            end
        end

    end

    % --------------------------------------------------------------------
    % * save best and crossover
    % --------------------------------------------------------------------
    bestp = paramvec{meffind};

    % test power of bestp
    clear delta
    for reg = 1:length(bestp)
        delta{reg} = zeros(mspec.numframes,1);
        delta{reg}(bestp{reg}+1) = 1;
        delta{reg} = delta{reg}(1:mspec.numframes);
    end
    P = hrf_power(mspec.TR,conditions(1).subcond.hrfest,delta,eye(length(bestp)));
    bestipow(i) = max(P.ipow);
    figure(f1); hold on; plot(bestipow,'r','LineWidth',2),drawnow
    %figure;subplot(1,2,1);imagesc(P.X2); subplot(1,2,2);imagesc(XX{meffind}); colormap gray

    % --------------------------------------------------------------------
    % * Check for convergence
    % --------------------------------------------------------------------
    if i > genconverge
        last_gens = meff(i-genconverge+1:i);
        unique_in_last_gens = length(unique(last_gens));
        % if the above is one, there has been no change in last n
        % generations
    else
        unique_in_last_gens = Inf;
    end

    if (all(eff == median(eff))) || unique_in_last_gens == 1
        isconverged = 1;
        disp(['System converged at generation ' num2str(i)])

        break
    end

    % save 1st generation in origvec output
    if i == 1, origvec = newvec; end

    [newvec] = xover(paramvec,eff,trialonsetinfo);
        %newvec{ceil(rand * length(newvec))} = bestp;   % this causes lots of the same list to be inserted randomly
     newvec{1} = bestp;

end

% --------------------------------------------------------------------
% * End of GA
% --------------------------------------------------------------------

X = construct_model(mspec,conditions,bestp,nooverlap);

% --------------------------------------------------------------------
% * Figure of final model matrix
% --------------------------------------------------------------------
try
    figure('Color','w'); imagesc(X); colormap(gray)
    bestcondno = cond(X);
    tmp = corrcoef(X); tmp=tmp(:); tmp(isnan(tmp) | tmp == 1) = []; maxcorr = max(tmp);
    tmp = sprintf('Best model: Cond. # = %3.2f, max pairwise corr = %3.2f',bestcondno,maxcorr);
    title(tmp),
catch
    disp('Problem with figures')
end

% this plots lines for onsets of first event type
%hold on; plot(repmat([.5 1],length(evtonsets{1}),1)',[1+(evtonsets{1} ./ mspec.TR) 1+(evtonsets{1} ./ mspec.TR)]','r','LineWidth',2)

% --------------------------------------------------------------------
% * Text list of onset times
% --------------------------------------------------------------------
try
    disp(' '); disp('GA2.m output')
    [evtonsets,paramvec] = paramvec2onsets(mspec,conditions,bestp,nooverlap);
catch
    disp('Error building evtonsets!  Incorrect model specification in model_setup_ui.m')
    evtonsets = [];
end

fname = ['ga2_output' datestr(now,31)]; fname(fname==' ')='_';fname(end-5)='h';fname(end-2)='m';fname=fname(1:end-2);
try
    eval(['save ' fname ' mspec conditions X bestp nooverlap evtonsets paramvec'])
catch
    warning('Not all vars created?  Saving everything.')
    eval(['save ' fname])
end

if showres

    % --------------------------------------------------------------------
    % * Figure of final model matrix
    % --------------------------------------------------------------------
    try
        [out,P] = ga3power(bestp,mspec,conditions,eye(length(numframes)));
    catch
        disp('Could not execute power plot.')
    end

    % --------------------------------------------------------------------
    % * Figure of A vs D optimality
    % --------------------------------------------------------------------
try
    figure('Color','w');hist(diff(evtonsets{1})); xlabel('Time between trial onsets'),ylabel('Frequency'); set(gca,'FontSize',18)

    if ~exist('Vi') == 1
        disp('No Vi detected: Using white noise')
        Vi = eye(mspec.numframes);      % use this for scanner white noise
    end

    disp('Comparing GA results (red square) to 1000 random designs (blue x)')
     for i = 1:1000,
        [X,paramvec] = construct_model(mspec,conditions,[],nooverlap);
        a(i) = calcefficiency(ones(1,size(X,2) - 1),[],pinv(X),Vi,0);
        d(i) = det(X'*X) .^ (1./size(X,2));     % 1 ./ (inv(det(X'*X)) .^ (1./size(X,2)));
     end
    figure('Color','w'); plot(a,d,'x'); set(gca,'FontSize',18); xlabel('A-efficiency'),ylabel('D-efficiency')
    [X,paramvec] = construct_model(mspec,conditions,bestp,nooverlap);
    a1 = calcefficiency(ones(1,size(X,2) - 1),[],pinv(X),Vi,0);
    d1 = 1 ./ (inv(det(X'*X)) .^ (1./size(X,2)));
    hold on;plot(a1,d1,'rs','MarkerFaceColor','r')
    title('GA result vs. random designs')
catch
    disp('Problem with A-vs-D optimality figure')
end

end

return

% --------------------------------------------------------------------
% * get best half, crossover to fill in lists
% --------------------------------------------------------------------
function newvec = xover(paramvec,eff,t)

% save best one - now done outside xover
%b = find(eff == max(eff)); b = b(1); b = paramvec{b};

% add a little random noise to efficiency - 1% of var of efficiency
eff = eff + randn(1,length(eff)) .* .01 * var(eff);

w = find(eff > median(eff));

% we can only do crossover if not all the designs are the same
if isempty(w)
    warning('Extremely homogenous sample!')
    % add a little random noise to efficiency - 1% of var of efficiency
    eff = eff + randn(1,length(eff)) .* .01 * mean(eff);
    w = find(eff > median(eff));
end

% check all paramvec values - works if all trial onsets are fixed
%for i=1:length(eff),ttmp=cat(1,paramvec{1}{:});ttmp=ttmp(:);,tmax(i)=max(ttmp);tmin(i)=min(ttmp);,end,[tmin; tmax]

% save best half
newvec = paramvec(w);
last = size(newvec,2);  % save size for adding random variation to existing

% fill in 2nd half with crossovers of first half
while length(newvec) < length(paramvec)

    w = ceil(rand(1,2) * length(newvec));
    babyv = [];

    for i = 1:length(newvec{w(1)})

        babyv{i} = rcomb(newvec{w(1)}{i},newvec{w(2)}{i});

        if length(t) > 1    % if we have full trialonsetinfo
        % older, more cumbersome way based on each condition
        % construct a set of legal vectors - cannot be less time between events than min onset time

        it = 1;
        while any(diff(babyv{i}) < t(1,i))

            tmp = [Inf diff(babyv{i})]; tmp = t(1,i) - tmp; tmp(tmp < 0) = 0;
            babyv{i} = babyv{i} + ceil(tmp);

            it = it + 1;
            if it > 20,
                %warning(['Too many iterations to get a legal vector for ' num2str(i)]),
                babyv{i} = rcomb(newvec{w(1)}{i},newvec{w(2)}{i});
            end

            if it > 100
                warning(['Too many iterations to get a legal vector for cond ' num2str(i)]),
                w = ceil(rand(1,2) * length(newvec));
                babyv{i} = rcomb(newvec{w(1)}{i},newvec{w(2)}{i});
                keyboard
            end
        end

        else
            % newer way constrains min time only
            % t is constrained in construct_model, so leave out here
        end

    end

    newvec{end+1} = babyv;

end

% add random variation to existing best half
for i = 1:length(paramvec{1})   % i indexes events

    for j = 1:last  % j indexes organisms

        p = sign(randn(1,length(newvec{j}{i})));

        if length(t) > 1    % if old way

        if t(3,i) == 1

            p(1) = 0;
            d = diff(newvec{j}{i});   % get diffs to make sure we don't move something illegally
            d = [Inf d] - t(1,i);
            p(find(d <= 1 & p < 0)) = 0; % set shift value to 0 for these points
            newvec{j}{i} = newvec{j}{i} + p;

            while any(diff(newvec{j}{i}) < t(1,i))

                tmp = [Inf diff(newvec{j}{i})]; tmp = t(1,i) - tmp; tmp(tmp < 0) = 0;
                newvec{j}{i} = newvec{j}{i} + ceil(tmp);

            end

        elseif t(3,i) < 1

            tmp2 = newvec{j}{i} + p;
            %disp('orig')
            %tmp2
            % ensure no values go outside proscribed range
            tmp2(tmp2 < t(1,i)) = t(1,i);
            tmp2(tmp2 > t(2,i)) = t(2,i);
            %disp('final')
            %tmp2

            newvec{j}{i} = tmp2;
        else
            warning('This should not happen.  t(3,:) should be 1 or -Inf')
        end

        else   % not old way, no constraints here (constraints in construct_model)
            newvec{j}{i} = newvec{j}{i} + p;
            newvec{j}{i}(newvec{j}{i} < 0) = 0;
        end

    end

    % warning: may not really ensure that we have stimuli too close
end

%newvec{1} = b;      % re-insert best one

return

% --------------------------------------------------------------------
% * crossover
% --------------------------------------------------------------------
function c = rcomb(a,b)
% combines 2 vectors of numbers at a random crossover point

w = ceil(rand * (length(a) - 1));
c(1:w) = a(1:w);
c = [c b(w+1:end)];

return

function [t,sc] = get_tonsetinfo(conditions)
% get list of onset times to match paramvec
% this contains info needed to see if one trial is impinging on the previous one, based on the length
% so it saves a vector, t, of the minimum times in TRs that each paramvec can have
% this value is compared, for each onset vector in paramvec, to the difference among onset times (trial length)
% to make sure no trial length is less than the minimum.
%
% t has first row = min, 2nd = max, 3rd = trial onsets (1) or subpart (0)

ind = 1;
for i = 1:length(conditions)
    if size(conditions(i).onsets,2) == 2   % variable range of trial onsets
        tmp = sum(conditions(i).onsets);
        t(1,ind) = tmp(1);
        t(2,ind) = tmp(2);
        t(3,ind) = 1;
        ind = ind + 1;
    end

    for j = 1:length(conditions(i).subcond)
        if length(conditions(i).subcond(j).onsets) > 1 & any(diff(conditions(i).subcond(j).onsets))
            t(1:2,ind) = [min(conditions(i).subcond(j).onsets); max(conditions(i).subcond(j).onsets)];
            t(3,ind) = 0;
            % need this for subcond, i think   %conditions(i).subcond(j).onsets(1);
            ind = ind + 1;
        end
    end
end

t(t == 0) = -Inf;

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [t,q,R,Rt] = splitFrame(F)

% SPLITFRAME  Split frame information into useful matrices and vectors.
%   [T,Q,R,Rt] = SPLITFRAME(F), for a frame F, returns the translation
%   vector T, quaternion Q, rotation matrix R and its transpose Rt. The
%   frame F can be either a 7-vector F=[T;Q] or a structure containing, at
%   least, the fields F.t, F.q, F.R and F.Rt.
%
%   See also FRAME, UPDATEFRAME.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

if isa(F,'struct')
    t  = F.t;
    q  = F.q;
    R  = F.R;
    Rt = F.Rt;

else % F is a 7-vector
    t  = F(1:3);
    q  = F(4:7);
    R  = q2R(q);
    Rt = R';

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:
function image = motionCompMeanImage(view, scans, frames)

%    function image = motionCompMeanImage(view, [scans], [frames])
%
% gb 01/22/05
%
% Finds the image closest to the mean image across frames for all scans
% Default : scans = current scan
%           frames = all frames

% Initilizes arguments and variables

global dataTYPES
curType = viewGet(view, 'curdatatype');

if ieNotDefined('scans')
    scans = 1:length(dataTYPES(curType).scanParams);
end
if ieNotDefined('frames')
    frames = 1:length(dataTYPES(curType).scanParams(scans(1)).nFrames);
end

% Computes the mean image
meanImage = motionCompComputeMean(view,scans,frames);

% Finds the image closest to the mean image
[scan,frame,mn] = motionCompNearestImage(view, scans, frames, meanImage);

% Loads this image
frame = frames(frame);
image = motionCompLoadImages(view,scan,frame);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
classdef OptimizationVariablesNearBoundsValidator < AbstractLaunchVehicleDataValidator
    %NoOptimizationVariablesValidator Summary of this class goes here
    %   Detailed explanation goes here

    properties
        lvdData LvdData
    end

    methods
        function obj = OptimizationVariablesNearBoundsValidator(lvdData)
            obj.lvdData = lvdData;
        end

        function [errors, warnings] = validate(obj)
            errors = LaunchVehicleDataValidationError.empty(0,1);
            warnings = LaunchVehicleDataValidationWarning.empty(0,1);

            varSet = obj.lvdData.optimizer.vars;
            varSet.sortVarsByEvtNum();
            [x, ~, varNameStrs, ~] = varSet.getTotalScaledXVector();
            [lb, ub] = varSet.getTotalScaledBndsVector();

            warnInd = [];
            for(i=1:length(x)) %#ok<*NO4LP>
                normX = (x(i) - lb(i)) / (ub(i) - lb(i));

                if(normX >= 0.99 || normX <= 0.01 || (ub(i) == lb(i) && x(i) == lb(i)))
                    warnInd(end+1) = i; %#ok<AGROW>
                end
            end

            if(isempty(warnInd))
                return;
            end

            for(i=1:length(warnInd))
                varNames{i} = varNameStrs{warnInd(i)}; %#ok<AGROW>
            end

            varNamesSubStr = strjoin(varNames,'\n');

            str = sprintf('Variables are near optimization bounds on some events:\n%s', varNamesSubStr);
            warnings(end+1) = LaunchVehicleDataValidationWarning(str);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [adfreq, n, ad] = plx_ad_span_v(filename, channel, startCount, endCount)
% plx_ad_span_v(filename, channel): read a span of a/d data from a .plx or .pl2 file
%
% [adfreq, n, ad] = plx_ad_span_v(filename, channel, startCount, endCount)
%
% INPUT:
%   filename - if empty string, will use File Open dialog
%   startCount - index of first sample to fetch
%   endCount - index of last sample to fetch
%   channel - 0 - based channel number or channel name
%
% OUTPUT:
%   adfreq - digitization frequency for this channel
%   n - total number of data points
%   ad - array of a/d values converted to mV

adfreq = 0;
n = 0;
ad = -1;

if nargin ~= 4
    error 'expected 4 input arguments';
end

[ filename, isPl2 ] = internalPL2ResolveFilenamePlx( filename );
if isPl2 == 1
    % return data from .pl2 file
    pl2Channel = channel;
    % PL2 uses 1-based channel numbers
    if ischar(channel) == 0
        pl2Channel = channel+1;
    end
    pl2ad = PL2AdSpan(filename, pl2Channel, startCount, endCount);
    if numel(pl2ad.Values) > 0
        adfreq = pl2ad.ADFreq;
        n = numel(pl2ad.Values);
        ad = pl2ad.Values;
    end
    return
end

% return plx data
channelNumber = plx_ad_resolve_channel(filename, channel);
if channelNumber == -1
    fprintf('\n plx_ad_span_v: no header for the specified A/D channel.');
    fprintf('\n    use plx_ad_info(filename) to print the list of a/d channels\n');
    return
end

[adfreq, n, ad] = mexPlex(18, filename, channelNumber, startCount, endCount);

end


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

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function H = PlotSwathe(bar,swathe,col,transp,nticks,fo,frequency,swatheonly)
% =========================================================================
% Plot a line with a shaded swathe
% =========================================================================
% H = PlotSwathe(bar,swathe,col,transp,nticks,fo,frequency,swatheonly)
% -------------------------------------------------------------------------
% INPUT
%   - bar: line to plot
%   - swathe: if a vector, draws symmetric errorbars. Otherwise draws
%     	asymmetric error bars, with first columnn being the upper bar and
%       second column being the lower bar
%--------------------------------------------------------------------------
% OPTIONAL INPUT
%   - col: color of the line/patch, can be rgb number (e.g. [1 0 0]
%       or string (e.g. 'blue')
%   - transp: 'transparent' for transparency, [empty] otherwise.
%       Default is [empty]
%   - nticks: number ticks on the horizontal axis
%   - fo: first observation (convention: 1987.00 = 1987Q1)
%   - frequency: quarterly ('q') or annual ('y') frequency. Default: 'q'
%   - swatheonly: set to 1 to plot only the swathe. Default: 0
%--------------------------------------------------------------------------
% OUPUT
%   - H.bar: handle for the median
%   - H.patch: handle for the patch
%   - H.swathe: H.swathe(1) handle for upp , H.swathe(2) handle for low
%   - H.ax: hande for the axis
%   - H.nobs: number of observations
% =========================================================================
% EXAMPLE
% x = 3*sin(1:50);
% PlotSwathe(x,1,'dark blue','transparent',2,1970)
% hold on
% y = cos(1:50);
% PlotSwathe(y,1,'light red','transparent',2,1970)
% =========================================================================
% Ambrogio Cesa Bianchi, March 2015
% ambrogio.cesabianchi@gmail.com

%% INITIALIZE VARIABLES
if ~exist('col','var')
    colrgb = rgb('dark grey');
else
    if isnumeric(col) % in case user provides already the rgb code
        colrgb = col;
    else % otherwise use rgb to get the code
        colrgb = rgb(col);
    end
end

if ~exist('transp','var')
    transp = '';
end

if ~exist('fo','var')
    do_dates=0;
else
    do_dates=1;
    if ~exist('frequency','var')
        frequency = 'q';
    end
end

if ~exist('swatheonly','var')
    swatheonly = 0;
end

% Number of observations
nobs = length(bar); H.nobs = nobs;

% Make sure data is row-vector
if size(bar,1)>1
    bar =bar';
end

% Check if swathe is vector and set error bands
if isvector(swathe)
    % Check if swathe has the right dimension [1 x T]
    if size(swathe,1)>1
        swathe = swathe';
        low = bar - swathe;
        upp = bar + swathe;
    else
        low = bar - swathe;
        upp = bar + swathe;
    end
else
    % Check if i has the right dimension [2 x T]
    if size(swathe,1)>2
        swathe = swathe';
        upp = swathe(1,:);
        low = swathe(2,:);
    else
        upp = swathe(1,:);
        low = swathe(2,:);
    end
end

% Create the x axis
xaxis = 1:nobs;

%% PLOT
% Initialize patch
edgeColor = colrgb+(1-colrgb)*0.55;
patchSaturation = 0.15; %How de-saturated or transp to make the patch
if strcmp(transp,'transparent')==1
    faceAlpha = patchSaturation;
    patchColor = colrgb;
    set(gcf,'renderer','openGL')
else
    faceAlpha = 1;
    patchColor = colrgb+(1-colrgb)*(1-patchSaturation);
    set(gcf,'renderer','painters')
end

% Plot the 'bar' line
if swatheonly==0
    H.bar = plot(xaxis,bar,'LineWidth',2,'Color',colrgb);
end

%Add the error-bar plot elements
holdStatus = ishold;
if ~holdStatus, hold on,  end

% Make the cordinats for the patch
yP = [low,fliplr(upp)];
xP = [xaxis,fliplr(xaxis)];

% Remove any nans otherwise patch won't work
xP(isnan(yP))=[];
yP(isnan(yP))=[];
H.patch=patch(xP,yP,1,'facecolor',patchColor,'edgecolor','none','facealpha',faceAlpha);

%Make nice edges around the patch.
H.swathe(1)=plot(xaxis,low,'-','color',edgeColor);
H.swathe(2)=plot(xaxis,upp,'-','color',edgeColor);
if swatheonly==0
    delete(H.bar)
end
if swatheonly==0
    H.bar=plot(xaxis,bar,'LineWidth',2,'Color',colrgb);
end
if ~holdStatus, hold off, end

% Make smart x labels
if do_dates==1
    DatesPlot(fo,nobs,nticks,frequency)
end
H.ax = gca;

author: ambropo

repo: VAR-Toolbox

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


text:
function channels = PW(N, beta)
m = log2(N);
channels = zeros(N, 1);
for i = 0 : N - 1
    bin_seq = dec2bin(i, m);
    sum = 0;
    for j = 1 : m
        if bin_seq(j) == '1'
            sum = sum + beta^(m - j);
        end
    end
    channels(i + 1) = sum;
end
end

author: YuYongRun

repo: PolarCodeDecodersInMatlab


text:
function [L, Nc] = discretize(S, N)
%DISCRETIZE   Discretize a SPINOP.
%   [L, NC] = DISCRETIZE(S, N) uses a Fourier spectral method in coefficient
%   space to discretize the SPINOP S with N grid points. L is the linear part, a
%   NxN diagonal matrix stored as a Nx1 vector, and NC is the differentiation
%   term of the nonlinear part (and hence is linear).
%
% Remark: DISCRETIZE will fail to discretize SPINOP objects which are not of the
%         right form. See HELP/SPINOP.
%
% See also SPINOP.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

%% Set-up:

% Get the domain DOM, the linear part FUNCL, the nonlinear part FUNCN, and the
% number of variables NVARS:
dom = S.domain;
funcL = S.lin;
funcNc = S.nonlinearPartCoeffs;
nVars = nargin(funcL);

% Get the variables of the workspace:
func = functions(funcL);
wrk = func.workspace{1};
names = fieldnames(wrk);
if ( isempty(names) == 0 )
    lengthNames = size(names, 1);
    for k = 1:lengthNames
        eval(sprintf('%s = wrk.(names{k});', names{k}));
    end
end

% Create a CHEBOPPREF object with TRIGSPEC discretization:
pref = cheboppref();
pref.discretization = @trigspec;

%% Discretize the linear part:

% Add the dependent variable 'x':
strL = func2str(funcL);
strL = strrep(strL, '@(', '@(x,');
funcL = eval(strL);

% USE CHEBOP, LINOP and MATRIX with TRIGSPEC discretization:
matL = matrix(linop(chebop(funcL, dom)), N, pref);
L = [];
for l = 1:nVars
    idx = (l-1)*N + 1;
    temp = full(diag(matL(idx:idx+N-1,idx:idx+N-1)));
    if ( mod(N, 2) == 0 )
        temp = fftshift(temp);
    else
        temp = ifftshift(temp);
    end
    L = [L; temp]; %#ok<*AGROW>
end

%% Disretize the differentiation term of the nonlinear part:

% If NC if a FUNCTION HANDLE, discretize it with TRIGSPEC:
if ( isa(funcNc, 'function_handle') == 1 )

    % Add the dependent variable 'x':
    strNc = func2str(funcNc);
    strNc = strrep(strNc, '@(', '@(x,');
    funcNc = str2func(strNc);

    % USE CHEBOP, LINOP and MATRIX with TRIGSPEC discretization:
    matNc = matrix(linop(chebop(funcNc, dom)), N, pref);
    Nc = [];
    for l = 1:nVars
        idx = (l-1)*N + 1;
        temp = full(diag(matNc(idx:idx+N-1,idx:idx+N-1)));
        if ( mod(N, 2) == 0 )
            temp = fftshift(temp);
        else
            temp = ifftshift(temp);
        end
        Nc = [Nc; temp]; %#ok<*AGROW>
    end

% NC might be 1 and in that case, nothing to do:
elseif ( isa(funcNc, 'double') == 1 )
    Nc = funcNc;
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function configCMatrix()
path = fileparts(mfilename('fullpath'));
addpath(fullfile(path, '..'));
config;

addpath(fullfile(path, 'include'));

author: opencobra

repo: cobratoolbox


text:
function [ n_data, x, fx ] = bessel_j0_int_values ( n_data )

%*****************************************************************************80
%
%% BESSEL_J0_INT_VALUES returns some values of the Bessel J0 integral.
%
%  Discussion:
%
%    The function is defined by:
%
%      J0_INT(x) = Integral ( 0 <= t <= x ) J0(t) dt
%
%    The data was reported by McLeod.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 August 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Allan McLeod,
%    Algorithm 757, MISCFUN: A software package to compute uncommon
%      special functions,
%    ACM Transactions on Mathematical Software,
%    Volume 22, Number 3, September 1996, pages 288-301.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real X, the argument of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 20;

  fx_vec = [ ...
      0.97656242238978822427E-03, ...
      0.39062450329491108875E-02, ...
     -0.62479657927917933620E-01, ...
      0.12483733492120479139E+00, ...
     -0.48968050664604505505E+00, ...
      0.91973041008976023931E+00, ...
     -0.14257702931970265690E+01, ...
      0.10247341594606064818E+01, ...
     -0.12107468348304501655E+01, ...
      0.11008652032736190799E+01, ...
     -0.10060334829904124192E+01, ...
      0.81330572662485953519E+00, ...
     -0.10583788214211277585E+01, ...
      0.87101492116545875169E+00, ...
     -0.88424908882547488420E+00, ...
      0.11257761503599914603E+01, ...
     -0.90141212258183461184E+00, ...
      0.91441344369647797803E+00, ...
     -0.94482281938334394886E+00, ...
      0.92266255696016607257E+00 ];

  x_vec = [ ...
       0.0009765625E+00, ...
       0.0039062500E+00, ...
      -0.0625000000E+00, ...
       0.1250000000E+00, ...
      -0.5000000000E+00, ...
       1.0000000000E+00, ...
      -2.0000000000E+00, ...
       4.0000000000E+00, ...
      -8.0000000000E+00, ...
      16.0000000000E+00, ...
     -16.5000000000E+00, ...
      18.0000000000E+00, ...
     -20.0000000000E+00, ...
      25.0000000000E+00, ...
     -30.0000000000E+00, ...
      40.0000000000E+00, ...
     -50.0000000000E+00, ...
      75.0000000000E+00, ...
     -80.0000000000E+00, ...
     100.0000000000E+00 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    x = 0.0;
    fx = 0.0;
  else
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [nev] = neuralynx_read_nev(filename, varargin)
% READ_NEURALYNX_NEV reads the event information from the *.nev file in a Neuralynx dataset directory
% Use as
%   nev = read_neuralynx_hdr(datadir, ...)
%   nev = read_neuralynx_hdr(eventfile, ...)
%
% Optional input arguments should be specified in key-value pairs and may include
%   implementation  should be 1, 2 or 3 (default = 3)
%   value           number or list of numbers
%   mintimestamp    number
%   maxtimestamp    number
%   minnumber       number
%   maxnumber       number
%
% The output structure contains all events and timestamps.

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

% get the optional input arguments
% flt_value        = ft_getopt(varargin, 'value');
% flt_mintimestamp = ft_getopt(varargin, 'mintimestamp');
% flt_maxtimestamp = ft_getopt(varargin, 'maxtimestamp');
% flt_minnumber    = ft_getopt(varargin, 'minnumber');
% flt_maxnumber    = ft_getopt(varargin, 'maxnumber');
% implementation   = ft_getopt(varargin, 'implementation', 3);
flt_value        = [];
flt_mintimestamp = [];
flt_maxtimestamp = [];
flt_minnumber    = [];
flt_maxnumber    = [];
implementation   = 3;

% if ft_filetype(filename, 'neuralynx_ds')
%   % replace the directory name by the filename
%   if     exist(fullfile(filename, 'Events.Nev'))
%     filename = fullfile(filename, 'Events.Nev');
%   elseif exist(fullfile(filename, 'Events.nev'))
%     filename = fullfile(filename, 'Events.nev');
%   elseif exist(fullfile(filename, 'events.Nev'))
%     filename = fullfile(filename, 'events.Nev');
%   elseif exist(fullfile(filename, 'events.nev'))
%     filename = fullfile(filename, 'events.nev');
%   end
% end

% The file starts with a 16*1024 bytes header in ascii, followed by a
% number of records (c.f. trials).
%
% The format of an event record is
%   int16 PktId
%   int16 PktDataSize
%   int64 TimeStamp
%   int16 EventId
%   int16 TTLValue
%   int16 CRC
%   int32 Dummy
%   int32 Extra[0]
%   ...
%   int32 Extra[7]
%   char EventString[0]
%   ...
%   char EventString[127]
% PktId is usually 0x1002.
% PktDataSize is random data.
% Dummy is random data.
% CRC may contain random data.
% Extra is user-defined data.
% TTLValue is the value sent to the computer on a parallel input port.

fid = fopen(filename, 'rb', 'ieee-le');

headersize = 16384;
offset     = headersize;
if ~isempty(flt_minnumber)
  offset = offset + (flt_minnumber-1)*184;
end

nev = [];

if implementation==1
  if ~isempty(flt_maxnumber)
    warning('filtering on maximum number not yet implemneted');
  end
  % this is the slow way of reading it
  % it also does not allow filtering
  fseek(fid, offset, 'bof');
  while ~feof(fid)
    nev(end+1).PktStart     = fread(fid, 1, 'int16');
    nev(end  ).PktId        = fread(fid, 1, 'int16');
    nev(end  ).PktDataSize  = fread(fid, 1, 'int16');
    nev(end  ).TimeStamp    = fread(fid, 1, 'uint64');
    nev(end  ).EventId      = fread(fid, 1, 'int16');
    nev(end  ).TTLValue     = fread(fid, 1, 'uint16');
    nev(end  ).CRC          = fread(fid, 1, 'int16');
    nev(end  ).Dummy        = fread(fid, 1, 'int32');
    nev(end  ).Extra        = fread(fid, 8, 'int32');
    nev(end  ).EventString  = fread(fid, 128, 'char');
  end
end

if implementation==2
  if ~isempty(flt_maxnumber)
    warning('filtering on maximum number not yet implemneted');
  end
  % this is a faster way of reading it and it is still using the automatic type conversion from Matlab
  fp = offset;
  fseek(fid, fp+ 0, 'bof'); PktStart       = fread(fid, inf, 'uint16', 184-2);
  num = length(PktStart);
  fseek(fid, fp+ 2, 'bof'); PktId          = fread(fid, num, 'uint16', 184-2);
  fseek(fid, fp+ 4, 'bof'); PktDataSize    = fread(fid, num, 'uint16', 184-2);
  fseek(fid, fp+ 6, 'bof'); TimeStamp      = fread(fid, num, 'uint64=>uint64', 184-8);
  fseek(fid, fp+14, 'bof'); EventId        = fread(fid, num, 'uint16', 184-2);
  fseek(fid, fp+16, 'bof'); TTLValue       = fread(fid, num, 'uint16', 184-2);
  fseek(fid, fp+18, 'bof'); CRC            = fread(fid, num, 'uint16', 184-2);
  fseek(fid, fp+22, 'bof'); Dummy          = fread(fid, num, 'int32' , 184-4);
  % read each of the individual extra int32 values and concatenate them
  fseek(fid, fp+24+0*4, 'bof'); Extra1     = fread(fid, num, 'int32', 184-4);
  fseek(fid, fp+24+1*4, 'bof'); Extra2     = fread(fid, num, 'int32', 184-4);
  fseek(fid, fp+24+2*4, 'bof'); Extra3     = fread(fid, num, 'int32', 184-4);
  fseek(fid, fp+24+3*4, 'bof'); Extra4     = fread(fid, num, 'int32', 184-4);
  fseek(fid, fp+24+4*4, 'bof'); Extra5     = fread(fid, num, 'int32', 184-4);
  fseek(fid, fp+24+5*4, 'bof'); Extra6     = fread(fid, num, 'int32', 184-4);
  fseek(fid, fp+24+6*4, 'bof'); Extra7     = fread(fid, num, 'int32', 184-4);
  fseek(fid, fp+24+7*4, 'bof'); Extra8     = fread(fid, num, 'int32', 184-4);
  Extra = [Extra1 Extra2 Extra3 Extra4 Extra5 Extra6 Extra7 Extra8];
  % read the complete data excluding header as char and cut out the piece with the EventString content
  fseek(fid, fp, 'bof'); EventString        = fread(fid, [184 num], 'char');
  EventString = char(EventString(57:184,:)');
end

if implementation==3
  % this is an even faster way of reading it
  if isempty(flt_minnumber)
    flt_minnumber = 1;
  end
  if isempty(flt_maxnumber)
    flt_maxnumber = inf;
  end
  if fseek(fid, offset, 'bof')~=0
    error(ferror(fid));
  end
  buf = fread(fid, (flt_maxnumber-flt_minnumber+1)*184, 'uint8=>uint8');
  [PktStart, PktId , PktDataSize , TimeStamp , EventId , TTLValue , CRC , Dummy , Extra1 , Extra2 , Extra3 , Extra4 , Extra5 , Extra6 , Extra7 , Extra8 , EventString] = ...
    cstructdecode(buf, 'int16', 'int16', 'int16', 'uint64', 'int16', 'uint16', 'int16', 'int32', 'int32', 'int32', 'int32', 'int32', 'int32', 'int32', 'int32', 'int32', 'char128');
  % the cstructdecode function does not respect the byte order for other-endian numbers
  if 0 % bigendian
    PktStart     = swapbytes(PktStart);
    PktId        = swapbytes(PktId);
    PktDataSize  = swapbytes(PktDataSize);
    TimeStamp    = swapbytes(TimeStamp);
    EventId      = swapbytes(EventId);
    TTLValue     = swapbytes(TTLValue);
  % CRC          = swapbytes(CRC);
  % Dummy        = swapbytes(Dummy);
  % Extra1       = swapbytes(Extra1);
  % Extra2       = swapbytes(Extra2);
  % Extra3       = swapbytes(Extra3);
  % Extra4       = swapbytes(Extra4);
  % Extra5       = swapbytes(Extra5);
  % Extra6       = swapbytes(Extra6);
  % Extra7       = swapbytes(Extra7);
  % Extra8       = swapbytes(Extra8);
  end
end

fclose(fid);

if implementation~=1
  % make a selection of events
  sel = true(size(TTLValue));
  if ~isempty(flt_mintimestamp)
    sel = sel & (TimeStamp>=flt_mintimestamp);
  end
  if ~isempty(flt_maxtimestamp)
    sel = sel & (TimeStamp<=flt_maxtimestamp);
  end
  if ~isempty(flt_value)
    tmp = false(size(sel));
    for i=1:length(flt_value)
      tmp = tmp | (TTLValue==flt_value(i));
    end
    sel = sel & tmp;
  end

  % restructure the data into a struct array, by first making cell-arrays out of it ...
  numsel = sum(sel);
  PktStart      = mat2cell(PktStart(sel)     , ones(1,numsel), 1);
  PktId         = mat2cell(PktId(sel)        , ones(1,numsel), 1);
  PktDataSize   = mat2cell(PktDataSize(sel)  , ones(1,numsel), 1);
  TimeStamp     = mat2cell(TimeStamp(sel)    , ones(1,numsel), 1);
  EventId       = mat2cell(EventId(sel)      , ones(1,numsel), 1);
  TTLValue      = mat2cell(TTLValue(sel)     , ones(1,numsel), 1);
  % CRC           = mat2cell(CRC(sel)          , ones(1,numsel), 1);
  % Dummy         = mat2cell(Dummy(sel)        , ones(1,numsel), 1);
  % Extra1        = mat2cell(Extra1(sel)       , ones(1,numsel), 1);
  % Extra2        = mat2cell(Extra2(sel)       , ones(1,numsel), 1);
  % Extra3        = mat2cell(Extra3(sel)       , ones(1,numsel), 1);
  % Extra4        = mat2cell(Extra4(sel)       , ones(1,numsel), 1);
  % Extra5        = mat2cell(Extra5(sel)       , ones(1,numsel), 1);
  % Extra6        = mat2cell(Extra6(sel)       , ones(1,numsel), 1);
  % Extra7        = mat2cell(Extra7(sel)       , ones(1,numsel), 1);
  % Extra8        = mat2cell(Extra8(sel)       , ones(1,numsel), 1);
  EventString   = mat2cell(EventString(sel,:), ones(1,numsel), 128);
  EventNumber   = mat2cell(find(sel) + flt_minnumber - 1, ones(1,numsel), 1); % this helps in the external filtering for BCI
  % ... and then convert the cell-arrays into a single structure
  nev = struct(...
    'PktStart'    , PktStart     , ...
    'PktId'       , PktId        , ...
    'PktDataSize' , PktDataSize  , ...
    'TimeStamp'   , TimeStamp    , ...
    'EventId'     , EventId      , ...
    'TTLValue'    , TTLValue     , ...
    'EventString' , EventString  , ...
    'EventNumber' , EventNumber);
end

% remove null values and convert to strings
for i=1:length(nev)
  nev(i).EventString = nev(i).EventString(find(nev(i).EventString));
  nev(i).EventString = char(nev(i).EventString(:)');
end

author: brainstorm-tools


text:
function tx = updateFastBoundaryClassifierFeatures(X, s1, s2, eid)
% tx = updateFastBoundaryClassifierFeatures(X, s1, s2, eid)
%
% Computes a simple set of features based on two regions s1 and s2, each
% containing one or more segments from X and separated by edglets given by
% eid.

ng = 5; % five geometric classes

tx = zeros([1 33], 'single');

f = 0;

%% Edge features
tx(f+1) = sum(X.edge.pb(eid).*X.edge.length(eid)) / sum(X.edge.length(eid));

area1 = X.region.area(s1);
area2 = X.region.area(s2);

f = f+1;

%% Region features

% area
sumarea1 = sum(area1);
sumarea2 = sum(area2);
tx(f+(1:2)) = [min([sumarea1 sumarea2], [], 2) max([sumarea1 sumarea2], [], 2)];

% color
meanColor1 = sum(X.region.colorMean(s1, :).*[area1 area1 area1], 1) / sumarea1;
meanColor2 = sum(X.region.colorMean(s2, :).*[area2 area2 area2], 1) / sumarea2;
tx(f+3) = sqrt(sum((meanColor1-meanColor2).^2, 2));

% position
left1 = min(X.region.x(s1, 1));  right1 = max(X.region.x(s1, 3));
left2 = min(X.region.x(s2, 1));  right2 = max(X.region.x(s2, 3));
bot1 = min(X.region.y(s1, 1));  top1 = max(X.region.y(s1, 3));
bot2 = min(X.region.y(s2, 1));  top2 = max(X.region.y(s2, 3));

tx(f+4) = top1-top2;
tx(f+5) = bot1-bot2;
tx(f+6) = top1-bot2;
tx(f+7) = bot1-top2;
tx(f+8) = top1-bot1;
tx(f+9) = top2-bot2;
tx(f+10) = left1-left2;
tx(f+11) = right1-right2;
tx(f+12) = right1-left1;
tx(f+13) = right2-left2;

% x alignment
x1 = [left1 right1];
x2 = [left2 right2];
tx(f+14) = (min([x1(:, 2) x2(:, 2)], [], 2)-max([x1(:, 1) x2(:, 1)], [], 2)) ./ ...
    (max([x1(:, 2) x2(:, 2)], [], 2)-min([x1(:, 1) x2(:, 1)], [], 2));

% y overlap
y1 = [bot1 top1];
y2 = [bot2 top2];
tx(f+15) = (min([y1(:, 2) y2(:, 2)], [], 2)-max([y1(:, 1) y2(:, 1)], [], 2)) ./ ...
    (max([y1(:, 2) y2(:, 2)], [], 2)-min([y1(:, 1) y2(:, 1)], [], 2));

f = f + 15;

%% 3D Geometry features

% geometric context features
gc = X.region.geomContext;
gc1 = sum(gc(s1, :).*repmat(area1, [1 size(gc, 2)]), 1) / sumarea1;
gc2 = sum(gc(s2, :).*repmat(area2, [1 size(gc, 2)]), 1) / sumarea2;

tx(f+(1:ng)) = gc1;
tx(f+ng+(1:ng)) = gc2;
tx(f+2*ng+(1:ng)) = gc1-gc2;
tx(f+3*ng+1) = sum(abs(gc1-gc2), 2)/2;

[maxval1, maxlab1] = max([gc1(:, 1) sum(gc1(:, 2:4), 2) gc1(:, 5)], [], 2);
[maxval2, maxlab2] = max([gc2(:, 1) sum(gc2(:, 2:4), 2) gc2(:, 5)], [], 2);
tx(f+3*ng+2) = (maxlab1-1)*3+ maxlab2;

f= f + 17;

author: Cloud-CV

repo: object-proposals


text:
function play(db,index)
%PLAY - plays a MatlabADT sentence.
%play(ADTobj,index)
%See also query, filterdb, read, play.
[data,smpr] = read(db,index);
for i=1:length(data)
    sound(data{i},smpr);
end
end

author: yueyuzhao

repo: gyrophone

save_path: github-repos/MATLAB/yueyuzhao-gyrophone


text:
function Y = selb(M1, M2, mp)
%Y = selb(M1, M2, mp) coefficient selection for base image
%
%    M1  - coefficients A
%    M2  - coefficients B
%    mp  - switch for selection type
%          mp == 1: select A
%          mp == 2: select B
%          mp == 3: average A and B
%
%    Y   - combined coefficients

%    (Oliver Rockinger 16.08.99)

switch (mp)
  case 1, Y = M1;
  case 2, Y = M2;
  case 3, Y = (M1 + M2) / 2;
  otherwise, error('unknown option');
end;

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function F= GetCharacteristicBasedFilter(level,dBand,dataClass)

F=cell(3,level);
for l=1:level
  cubeSize=8;
  [P ,PF]=GeneratePyramidSection(cubeSize);

  numDir=dBand{level}{l}(1,1);
  shift=cubeSize/numDir;
  F{1,l}=cell(numDir,numDir);
  F{2,l}=cell(numDir,numDir);
  F{3,l}=cell(numDir,numDir);
  SF=PF{1}+PF{2}+PF{3};
  mRadial =ones(shift,cubeSize,shift);

  A=zeros(cubeSize,cubeSize,cubeSize);

   for l2=1:numDir
    for l1=1:numDir
      for c=1:3
      mRadialIdx=[(l1-1)*shift+1  l1*shift ;1 cubeSize ; (l2-1)*shift+1  l2*shift  ];

      F{c,l}{l2,l1}=PolarToRec(mRadialIdx,mRadial,cubeSize,P{c},SF);
       A=A+F{c,l}{l2,l1};
%       F{c,l}{l2,l1}=squeeze(real(fftshift(ifftn(fftshift(F{c,l}{l2,l1})))));

%       F{2,l}{l2,l1}=PolarToRec(mRadialIdx,mRadial,cubeSize,P{2},SF);
%       A=A+F{2,l}{l2,l1};
%       F{2,l}{l2,l1}=real(fftshift(ifftn(fftshift(F{2,l}{l2,l1}))));
%       F{3,l}{l2,l1}=PolarToRec(mRadialIdx,mRadial,cubeSize,P{3},SF);
%       A=A+F{3,l}{l2,l1};
%       F{3,l}{l2,l1}=real(fftshift(ifftn(fftshift(F{3,l}{l2,l1}))));
      end
    end
   end
for c=1:3
  for l2=1:numDir
    for l1=1:numDir
       F{c,l}{l2,l1}=F{c,l}{l2,l1}./A;
        F{c,l}{l2,l1}=squeeze(real(fftshift(ifftn(fftshift(F{c,l}{l2,l1})))));
     end
   end
 end
end

Temp=F(1,:);
F(1,:)=F(2,:);
F(2,:)=Temp;

author: thomaskuestner


text:
function T=transitivity_bd(A)
%TRANSITIVITY_BD    Transitivity
%
%   T = transitivity_bd(A);
%
%   Transitivity is the ratio of 'triangles to triplets' in the network.
%   (A classical version of the clustering coefficient).
%
%   Input:      A       binary directed connection matrix
%
%   Output:     T       transitivity scalar
%
%   Reference:  Rubinov M, Sporns O (2010) NeuroImage 52:1059-69
%               based on Fagiolo (2007) Phys Rev E 76:026107.
%
%
%   Contributors:
%   Mika Rubinov, UNSW/University of Cambridge
%   Christoph Schmidt, Friedrich Schiller University Jena
%   Andrew Zalesky, University of Melbourne
%   2007-2015

%   Modification history:
%   2007: original (MR)
%   2013, 2015: removed tests for absence of nodewise 3-cycles (CS,AZ)

%   Methodological note: In directed graphs, 3 nodes generate up to 8
%   triangles (2*2*2 edges). The number of existing triangles is the main
%   diagonal of S^3/2. The number of all (in or out) neighbour pairs is
%   K(K-1)/2. Each neighbour pair may generate two triangles. "False pairs"
%   are i<->j edge pairs (these do not generate triangles). The number of
%   false pairs is the main diagonal of A^2. Thus the maximum possible
%   number of triangles = (2 edges)*([ALL PAIRS] - [FALSE PAIRS])
%                       = 2 * (K(K-1)/2 - diag(A^2))
%                       = K(K-1) - 2(diag(A^2))

S    = A+A.';                           % symmetrized input graph
K    = sum(S,2);                        % total degree (in + out)
cyc3 = diag(S^3)/2;                     % number of 3-cycles (ie. directed triangles)
CYC3 = K.*(K-1)-2*diag(A^2);            % number of all possible 3-cycles
T    = sum(cyc3)./sum(CYC3);            % transitivity

author: fieldtrip

repo: fieldtrip