text:
function data = funInterp2(x_pred, y_pred, x_obs, y_obs, data_obs, fun)
%
% SINTAX:
%   data = funInterp2(x_pred, y_pred, x_obs, y_obs, data_obs, fun)
%
% INPUT:
%   x_pred , y_pred     coordinates of the prediction point (arrays [n x 1])
%   x_obs, y_obs        coordinates of the observation (array [n x 1], sparse points)
%   data_obs            data observation array [n x n_epochs]
%   fun                 virtual function f(dist) as function of the distance
%
% OUTPUT:
%   data                interpolated data in the point of interest [n x n_epochs]
%
% DEFAULT VALUES:
%   generic interpolator using as correlation function fun
%   <default: fun = @(dist) exp(-dist)>
%
% EXAMPLE:
%   fun = @(dist) 0.2 * exp(-(dist/1e4)) + exp(-(dist/6e3).^2);
%   temp = funInterp2(ep(:), np(:), e_obs(:), n_obs(:), td_obs(:,:), fun);
%

%--- * --. --- --. .--. ... * ---------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0RC1
%
%--------------------------------------------------------------------------
%  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
%  Written by:       Andrea Gatti
%  Contributors:     Andrea Gatti ...
%  A list of all the historical goGPS contributors is in CREDITS.nfo
%--------------------------------------------------------------------------
%
%   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/>.
%
%--------------------------------------------------------------------------
% 01100111 01101111 01000111 01010000 01010011
%--------------------------------------------------------------------------

    narginchk(5, 6);

    if nargin < 6
        % Correlation function
        fun = @(dist) exp(-dist);
    end

    % Init out data
    data = nan(size(data_obs, 2), numel(x_pred));

    [x_mesh, y_mesh] = meshgrid(x_obs, y_obs);
    d_obs = sqrt(abs(x_mesh - x_mesh').^2 + abs(y_mesh - y_mesh').^2);
    q_fun_obs = fun(d_obs);
    %q_fun_obs = exp(-d_obs/0.4e4);
    [xv, ~, xi] = unique(x_pred);
    [yv, ~, yi] = unique(y_pred);
    x2 = (repmat(x_obs, 1, numel(xv)) - repmat(xv', numel(x_obs), 1)).^2;
    y2 = (repmat(y_obs, 1, numel(yv)) - repmat(yv', numel(y_obs), 1)).^2;
    for i = 1 : numel(x_pred)
        d_pred = sqrt(x2(:,xi(i)) + y2(:,yi(i)));
        c_mat = q_fun_obs .* repmat(fun(d_pred)', size(q_fun_obs,1),1);
        %c_mat = bsxfun(@times, q_fun_obs, fun(d_pred)');
        c_mat = triu(c_mat) + triu(c_mat, 1)';

        trans = sum(c_mat);
        w = trans / sum(trans);
        data(:, i) = (w * data_obs)';
    end
end

author: goGPS-Project


text:
classdef PTKMaximumFissurePoints < PTKPlugin
    % PTKMaximumFissurePoints. Plugin which is part of the lobar segmentation.
    %
    %     This is a plugin for the Pulmonary Toolkit. Plugins can be run using
    %     the gui, or through the interfaces provided by the Pulmonary Toolkit.
    %     See PTKPlugin.m for more information on how to run plugins.
    %
    %     Plugins should not be run directly from your code.
    %
    %     PTKMaximumFissurePoints is an intermediate stage in segmenting the
    %     lobes.
    %
    %     For more information, see
    %     [Doel et al., Pulmonary lobe segmentation from CT images using
    %     fissureness, airways, vessels and multilevel B-splines, 2012]
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    properties
        ButtonText = 'Fissureness <br>Maxima'
        ToolTip = ''
        Category = 'Fissures'

        AllowResultsToBeCached = true
        AlwaysRunPlugin = false
        PluginType = 'ReplaceOverlay'
        HidePluginInDisplay = false
        FlattenPreviewImage = true
        PTKVersion = '1'
        ButtonWidth = 6
        ButtonHeight = 2
        GeneratePreview = true
        Visibility = 'Developer'
    end

    methods (Static)
        function results = RunPlugin(application, reporting)
            fissure_approximation = application.GetResult('PTKFissureApproximation');

            left_and_right_lungs = application.GetResult('PTKLeftAndRightLungs');

            results_left = PTKMaximumFissurePoints.GetLeftLungResults(application, fissure_approximation);
            results_right = PTKMaximumFissurePoints.GetRightLungResults(application, fissure_approximation);

            results = PTKCombineLeftAndRightImages(application.GetTemplateImage(PTKContext.LungROI), results_left, results_right, left_and_right_lungs);
            results.ImageType = PTKImageType.Colormap;
        end

        function results = GenerateImageFromResults(results, ~, ~)
        end
    end

    methods (Static, Access = private)
        function left_results = GetLeftLungResults(application, fissure_approximation)
            fissureness_roi = application.GetResult('PTKFissurenessROI');

            lung_mask = application.GetResult('PTKLeftAndRightLungs');
            lung_mask.ChangeRawImage(uint8(lung_mask.RawImage == 2));

            left_lung_roi = application.GetResult('PTKGetLeftLungROI');
            fissure_approximation = fissure_approximation.Copy;
            fissure_approximation.ResizeToMatch(left_lung_roi);
            lung_mask.ResizeToMatch(left_lung_roi);

            [high_fissure_indices, ref_image] = PTKGetMaxFissurePoints(fissure_approximation.RawImage == 6, lung_mask, fissureness_roi.LeftMainFissure, left_lung_roi, left_lung_roi.ImageSize);

            left_results = left_lung_roi.BlankCopy;
            left_results.ChangeRawImage(ref_image);
        end

        function right_results = GetRightLungResults(application, fissure_approximation)
            fissureness_roi = application.GetResult('PTKFissurenessROI');

            lung_mask = application.GetResult('PTKLeftAndRightLungs');
            lung_mask.ChangeRawImage(uint8(lung_mask.RawImage == 1));

            right_lung_roi = application.GetResult('PTKGetRightLungROI');
            fissure_approximation = fissure_approximation.Copy;
            fissure_approximation.ResizeToMatch(right_lung_roi);
            lung_mask.ResizeToMatch(right_lung_roi);

            [high_fissure_indices, ref_image_o] = PTKGetMaxFissurePoints(fissure_approximation.RawImage == 2, lung_mask, fissureness_roi.RightMainFissure, right_lung_roi, right_lung_roi.ImageSize);

            [high_fissure_indices_m, ref_image_m] = PTKGetMaxFissurePoints(fissure_approximation.RawImage == 3, lung_mask, fissureness_roi.RightMidFissure, right_lung_roi, right_lung_roi.ImageSize);

            ref_image = ref_image_m;
            ref_image(ref_image_o == 1) = 1;
            ref_image(ref_image_m == 1) = 8;

            right_results = right_lung_roi.BlankCopy;
            right_results.ChangeRawImage(ref_image);
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
% The COBRAToolbox: testPartition.m
%
% Purpose:
%     - tests the basic functionality of swiftcc++ algorithm
%
% Authors:
%     - Original file: Mojtaba Tefagh, March 2019
%

global CBTDIR

% require the specified toolboxes and solvers
solvers = prepareTest('needsLP', true, 'useSolversIfAvailable', {'gurobi'}, ...
    'requiredToolboxes', {'bioinformatics_toolbox'});

% save the current path
currentDir = pwd;

% initialize the test
testPass = fileparts(which('testPartition.m'));
cd(testPass);

% load the model
model = getDistributedModel('ecoli_core_model.mat');
model.rev = double(model.lb < 0);
A = fastcc(model, 1e-4, 0);

% set the cobra solver
solverLPOK = changeCobraSolver(solvers.LP{1}, 'LP', 0);

fprintf('\n -- Running swiftcc++ using the %s solver...\n\n', solvers.LP{1});
component = partition(model, solvers.LP{1}, 'swift');
assert(all(component(A)));
assert(length(A) == sum(component));

% output a success message
fprintf('\nDone.\n');

% change the directory
cd(currentDir)

author: opencobra

repo: cobratoolbox


text:
function out=bsxfun(f,x,y)
%BSXFUN

dim1 = size(x);
dim2 = size(y);

X = x;
Y = y;

for i = 1:length(dim1)
    if dim1(i) < dim2(i)
        d = ones(1,length(dim1));
        d(i) = dim2(i);
        X = repmat(X,d);
    elseif dim1(i) > dim2(i)
        d = ones(1,length(dim1));
        d(i) = dim1(i);
        Y = repmat(Y,d);
    end
end

% And compute
out = f(X,Y);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
%-------------------------------------------------------------------------------------------------------------------%
%
% IB2d is an Immersed Boundary Code (IB) for solving fully coupled non-linear
% 	fluid-structure interaction models. This version of the code is based off of
%	Peskin's Immersed Boundary Method Paper in Acta Numerica, 2002.
%
% Author: Nicholas A. Battista
% Email:  nick.battista@unc.edu
% Date Created: May 27th, 2015
% Institution: UNC-CH
%
% This code is capable of creating Lagrangian Structures using:
% 	1. Springs
% 	2. Beams (*torsional springs)
% 	3. Target Points
%	4. Muscle-Model (combined Force-Length-Velocity model, "HIll+(Length-Tension)")
%
% One is able to update those Lagrangian Structure parameters, e.g., spring constants, resting %%	lengths, etc
%
% There are a number of built in Examples, mostly used for teaching purposes.
%
% If you would like us %to add a specific muscle model, please let Nick (nick.battista@unc.edu) know.
%
%--------------------------------------------------------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: finds the indices on the Eulerian grid where the 1D Dirac-delta
% kernel is possibly non-zero in BOTH (x,y) directions
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [xInds,yInds] = give_NonZero_Delta_Indices_XY(xLag, yLag, Nx, Ny, dx, dy, supp)

%xLag: gives x-coordinate of Lagrangian position
%yLag: gives y-coordinate of Lagrangian position
%Nx:   # of Eulerian grid pts. in x-dimension
%Ny:   # of Eulerian grid pts. in y-dimension
%dx:   spatial-step along x-dimension of Eulerian grid
%dy:   spatial-step along y-dimension of Eulerian grid
%supp: size of support of the Dirac-delta kernel (should be even)

%Give x-dimension Non-Zero Delta Indices
xIndsAux = give_1D_NonZero_Delta_Indices(xLag, Nx, dx, supp);

%Repeat x-Indices for Non-Zero y-Indices!
xInds = [];
for i=1:supp
   xInds = [xInds xIndsAux]; %Sets up x-INDEX matrix bc we consider BOTH dimensions
end

%Give y-dimension Non-Zero Delta Indices
yIndsAux = give_1D_NonZero_Delta_Indices(yLag, Ny, dy, supp);

%Repeat y-Indices for Non-Zero x-Indices!
yInds = [];
for i=1:supp
    for j=1:supp
        yInds = [yInds yIndsAux(:,i)]; %Sets up y-INDEX matrix bc we consider BOTH dimensions
    end
end

author: nickabattista


text:
function ypost = fmri_spatfil(ypre, SpatialFil)
%
% ypost = fmri_spatfil(ypre, SpatialFilter)
%
% Spatially filters a functional slice.
%
% -------- Input Arguments ------
% 1. ypre - functional slice to be processed. Dimension is
%       [nRows nCols nTP nRuns].
% 2. SpatialFilter - coefficients of spatial filter (2D).
%
% -------- Output Arguments ------
% 1. ypost - spatially filtererd slice of the same dimension
%       as the preprocessed slice, ie, [nRows nCols nTP nRuns].
%
% Douglas Greve (greve@nmr.mgh.harvard.edu)
% January 25, 1999
% February 12, 1999
%
% $Id: fmri_spatfilter.m,v 1.1 2004/03/11 01:30:19 sayres Exp $
%

% Check that the number of input arguments is correct %
if(nargin ~= 2)
  msg = 'ypost = fmri_spatfil(ypre, SpatialFilter)';
  qoe(msg);error(msg);
end

% Get function slice dimensions %
nRows = size(ypre,1);
nCols = size(ypre,2);
nTP   = size(ypre,3);
nRuns = size(ypre,4);
nVoxels = nRows*nCols;

y = ones(nRows,nCols);
EdgeCorrection = conv2(y,SpatialFil,'same');

%%% Initialize matricies %%
ypost = zeros(nRows,nCols,nTP,nRuns);

%% Pass through each run %%
for r = 1:nRuns,

  %% Put slice in a temporary variable %%
  y = ypre(:,:,:,r);

  %% Pass through each time point %%
  for n = 1:nTP,
     z = conv2(y(:,:,n),SpatialFil,'same') ./ EdgeCorrection;
     y(:,:,n) = z;
   end

  ypost(:,:,:,r) = reshape(y, [nRows nCols nTP]);

end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function result = file_exist(file_name)

    result = (exist(file_name, 'file') == 2);

end

author: layumi

repo: Image-Text-Embedding


text:
% Set system limits
lims = mr.opts('MaxGrad',32,'GradUnit','mT/m',...
    'MaxSlew',130,'SlewUnit','T/m/s',...
    'rfRingdownTime', 30e-6, 'rfDeadtime', 100e-6);

seq=mr.Sequence(lims);          % Create a new sequence object
fov=220e-3; Nx=256; Ny=256;     % Define FOV and resolution

foe=200e-3;             % Field of excitation
targetWidth=22.5e-3;    % Diameter of target excitation pattern
n=8;                    % Number of spiral turns
T=8e-3;                 % Pulse duration

% Define spiral k-space trajectory
kMax=(2*n)/foe/2;       % Units of 1/m (not rad/m)
tk=0:seq.gradRasterTime:T-seq.gradRasterTime;
kx=kMax*(1-tk/T).*cos(2*pi*n*tk/T);
ky=kMax*(1-tk/T).*sin(2*pi*n*tk/T);

% Define RF pulse
tr=0:seq.rfRasterTime:T-seq.rfRasterTime;
kxRf=interp1(tk,kx,tr,'linear','extrap');
kyRf=interp1(tk,ky,tr,'linear','extrap');
beta=2*pi*kMax*targetWidth/2/sqrt(2);  % Gaussian width in k-space
signal0 = exp(-beta.^2.*(1-tr/T).^2).*sqrt((2*pi*n*(1-tr/T)).^2+1);
signal = signal0.*(1 + exp(-1j.*2*pi*5e-2*(kxRf + kyRf)));

% Add gradient ramps
[kx,ky,signal]=mr.addRamps({kx,ky},'rf',signal);

rf = mr.makeArbitraryRf(signal,20*pi/180,'system',lims);
gxRf = mr.makeArbitraryGrad('x',mr.traj2grad(kx));
gyRf = mr.makeArbitraryGrad('y',mr.traj2grad(ky));

% Define other gradients and ADC events
deltak=1/fov;
gx = mr.makeTrapezoid('x','FlatArea',Nx*deltak,'FlatTime',6.4e-3);
adc = mr.makeAdc(Nx,'Duration',gx.flatTime,'Delay',gx.riseTime);
gxPre = mr.makeTrapezoid('x','Area',-gx.area/2,'Duration',2e-3);
phaseAreas = ((0:Ny-1)-Ny/2)*deltak;

% Refocusing pulse and spoiling gradients
%[rf180, gz] = mr.makeBlockPulse(pi,'Duration',1e-3,'SliceThickness',5e-3);
[rf180, gz] = mr.makeSincPulse(pi,'system',lims,'Duration',3e-3,...
    'SliceThickness',5e-3,'apodization',0.5,'timeBwProduct',4);

gzSpoil = mr.makeTrapezoid('z','Area',gx.area,'Duration',2e-3);

%%% Calculate timing (TE=20ms, TR=500ms)
delayTE1=ceil((20e-3/2 - mr.calcDuration(gzSpoil) - mr.calcDuration(rf180)/2)/seq.gradRasterTime)*seq.gradRasterTime;
delayTE2=delayTE1 - mr.calcDuration(gxPre) - mr.calcDuration(gx)/2;
delayTR=500e-3 - 20e-3 - mr.calcDuration(rf) - mr.calcDuration(gx)/2;

%% Loop over phase encodes and define sequence blocks
for i=1:Ny
    seq.addBlock(rf,gxRf,gyRf);
    seq.addBlock(mr.makeDelay(delayTE1));
    seq.addBlock(gzSpoil);
    seq.addBlock(rf180,gz);
    seq.addBlock(gzSpoil);
    seq.addBlock(mr.makeDelay(delayTE2));
    gyPre = mr.makeTrapezoid('y','Area',phaseAreas(i),'Duration',2e-3);
    seq.addBlock(gxPre,gyPre);
    seq.addBlock(gx,adc);
    seq.addBlock(mr.makeDelay(delayTR));
end

%% check whether the timing of the sequence is correct
[ok, error_report]=seq.checkTiming;

if (ok)
    fprintf('Timing check passed successfully\n');
else
    fprintf('Timing check failed! Error listing follows:\n');
    fprintf([error_report{:}]);
    fprintf('\n');
end

%%
seq.setDefinition('Name', 'se_selRF');

seq.write('selectiveRf.seq');   % Write to pulseq file
seq.plot();

return

%% Write to file
% The sequence is written to file in compressed form according to the file
% format specification using the |write| method.
%seq.setDefinition('Scan_ID',2068);
%seq.setDefinition('Recon_Mode',1);
%seq.write('external.seq')

%seq.plot

author: pulseq

repo: pulseq

save_path: github-repos/MATLAB/pulseq-pulseq


text:
% hh = plot_onsets(d1, [color], [miny], [scalef], [duration])
%
% plots non-zero elements of an indicator vector or a list of onset times as vertical lines
%
% d1 is indicator vector; uses entries as scaling for height
%
% optional inputs:
% color, e.g., 'r'
% miny, shifts bars up/down to miny
% scalef, scales height of bars
%
% e.g.,
% plot_onsets(d1, 'k', -.5, .4)
%
%

function hh = plot_onsets(d1, varargin)
    if ~any(d1(2:end) == 0)    % we have an onset list, get indicator
        n = round(max(d1));
        z = zeros(n, 1);
        z(round(d1)+1) = 1;   % 0 is start of 1st onset, 1 of 1st element
        d1 = z;             % convert to indicator
        %wh = d1;
    end

    % make a row vector
    if size(d1, 1) > size(d1, 2)
        d1 = d1';
    end

    % indicator
    wh = find(d1);

    color = 'k';
    miny = 0;
    scalef = 1;
    dur = 0;

    if length(varargin) > 0, color = varargin{1}; end
    if length(varargin) > 1, miny = varargin{2}; end
    if length(varargin) > 2, scalef = varargin{3}; end
    if length(varargin) > 3, dur = varargin{4}; end

    % adjust by -1 so first element is time 0
    if dur == 0
        % events
        hh = plot([wh; wh]-1, [zeros(size(wh))+miny; miny+scalef.*d1(wh)], color);
    else
        % epochs
        hh = [];
        for i = 1:length(wh)
            hh(end+1) = drawbox(wh(i)-1, dur, miny, scalef.*d1(wh(i)), color);
        end
    end
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
% Test file for @chebfun/isequal.m.

function pass = test_isequal(pref)

% Get preferences.
if ( nargin < 1 )
    pref = chebfunpref();
end

pref.splitting = 1;

% Check empty case.
f = chebfun();
pass(1) = isequal(f, f);

% Check self equality.
f = chebfun(@(x) sin(x).*abs(x - 0.1), [-1 1], pref);
pass(2) = isequal(f, f);

% Check inequality for chebfuns with different row/column orientation.
pass(3) = ~isequal(f, f.');

% Check inequality for chebfuns of different dimensions.
g = chebfun(@(x) sin(x), [-1 1], pref);
pass(4) = ~isequal(f, g);

g = chebfun(@(x) [sin(x).*abs(x - 0.1) cos(x)], [-1 1], pref);
pass(5) = ~isequal(f, g);

% Check inequality for chebfuns with different domains.
g = chebfun(@(x) sin(x).*abs(x - 0.1), [-1+eps, 1+eps], pref);
pass(6) = ~isequal(f, g);

% Check inequality for chebfuns built from different functions.
g = chebfun(@(x) sin(x).*abs(x - 0.1) + eps, [-1 1], pref);
pass(7) = ~isequal(f, g);

%% Test on singular function:
dom = [-2 7];
pow = -1.64;
f = chebfun(@(x) sin(100*x).*(x-dom(1)).^pow, dom, 'exps', [pow 0], ...
    'splitting', 'on');
g = chebfun(@(x) sin(100*x).*(x-dom(1)).^pow, dom, 'exps', [pow 0], ...
    'splitting', 'off');
pass(8) = ~isequal(f, g);

%% Tests for functions defined on unbounded domain:

% Functions on [-inf inf]:

% Set the domain:
dom = [-Inf Inf];

op = @(x) (1-exp(-x.^2))./x;
f = chebfun(op, dom);
pass(7) = isequal(f, f);

% Blow-up function:
op = @(x) x.^2.*(1-exp(-x.^2));
g = chebfun(op, dom, 'exps', [2 2]);
pass(8) = ~isequal(f, g);
pass(9) = ~isequal(g, f);

%% Functions on [-inf b]:

% Set the domain:
dom = [-Inf -3*pi];

% Array-valued function:
op = @(x) [exp(x) x.*exp(x) (1-exp(x))./x];
f = chebfun(op, dom);
pass(10) = isequal(f, f);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function poolBconv()
    global config mem;
    curr_layer_idx = config.misc.current_layer;
    mem.grads{curr_layer_idx+1} = mem.deltas{curr_layer_idx+1} * mem.activations{curr_layer_idx}';
    mem.grads{curr_layer_idx+1+config.layer_num} = sum(mem.deltas{curr_layer_idx+1}, 2);
	mem.delta_act = config.weights{curr_layer_idx+1}' * mem.deltas{curr_layer_idx+1};
end

author: jimmy-ren


text:
function square_hex_grid_test05 ( )

%*****************************************************************************80
%
%% SQUARE_HEX_GRID_TEST05 tests HEX_GRID_01_APPROXIMATE_N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 January 2009
%
%  Author:
%
%    John Burkardt
%
  test_num = 4

  n_goal_test = [ 100, 200, 500, 10000 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'SQUARE_HEX_GRID_TEST05\n' );
  fprintf ( 1, '  For a hexagonal grid of points in the unit box,\n' );
  fprintf ( 1, '  HEX_GRID_01_APPROXIMATE_N seeks the value of\n' );
  fprintf ( 1, '  NODES_PER_LAYER that produces a mesh of about N points.\n' );
  fprintf ( 1, '\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  N_GOAL  NODES_PER_LAYER       N\n' );
  fprintf ( 1, '\n' );

  for test = 1 : test_num

    n_goal = n_goal_test ( test );

    [ nodes_per_layer, n ] = hex_grid_01_approximate_n ( n_goal );

    fprintf ( 1, '  %6d           %6d  %6d\n', n_goal, nodes_per_layer, n );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function itree = pruefer_to_tree_2 ( nnode, iarray )

%*****************************************************************************80
%
%% PRUEFER_TO_TREE_2 produces the edge list of a tree from its Pruefer code.
%
%  Discussion:
%
%    One can thus exhibit all trees on N nodes, produce
%    one at random, find the M-th one on the list, etc, by
%    manipulating the Pruefer codes.
%
%    For every labeled tree on N nodes, there is a unique N-2 tuple
%    of integers A1 through AN-2, with each A between 1 and N.  There
%    are N^(N-2) such sequences, and each one is associated with exactly
%    one tree.
%
%    This routine apparently assumes that the Pruefer code is
%    generated by taking the LOWEST labeled terminal node each time.
%    This is not consistent with PRUEFER_TO_TREE and TREE_TO_PRUEFER.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 June 2013
%
%  Author:
%
%    Original FORTRAN77 version by Albert Nijenhuis, Herbert Wilf.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Albert Nijenhuis. Herbert Wilf,
%    Combinatorial Algorithms,
%    Academic Press, 1978, second edition,
%    ISBN 0-12-519260-6.
%
%  Parameters:
%
%    Input, integer NNODE, number of nodes in desired tree.
%
%    Input, integer IARRAY(NNODE).  IARRAY(I), I = 1, NNODE-2
%    is the Pruefer code for the tree.
%
%    Output, integer ITREE(NNODE); the I-th edge of the tree
%    joins nodes I and ITREE(I).
%
  itree(1:nnode) = 0;

  for i = nnode-2 : -1 : 1

    l = iarray(i);

    if ( itree(l) == 0 )
      iarray(i) = - l;
      itree(l) = - 1;
    end

  end

  iarray(nnode-1) = nnode;
%
%  Find next index K so that ITREE(K) is 0.
%
  k = 1;
  j = 0;

  while ( itree(k) ~= 0 )
    k = k + 1;
  end

  kp = k;

  while ( 1 )

    j = j + 1;
    ir = abs ( iarray(j) );
    itree(kp) = ir;

    if ( j == nnode - 1 )
      break
    end

    if ( 0 < iarray(j) )
      while ( itree(k) ~= 0 )
        k = k + 1;
      end
      kp = k;
      continue
    end

    if ( k < ir )
      itree(ir) = 0;
      while ( itree(k) ~= 0 )
        k = k + 1;
      end
      kp = k;
      continue
    end

    kp = ir;

  end
%
%  Restore the signs of IARRAY.
%
  iarray(1:nnode-2) = abs ( iarray(1:nnode-2) );

  return
end


text:
%--- help for mdd/swz ---
%
%mdd/swz is a function.
%    [log_mdd, extras] = swz(obj, swz_pvalue, opts)
%

author: jmaih


text:
function ErrorCode = parseDebugSignals(mloc_statemachine_debug, ...
    mvdc_path_matching_debug, mvdc_trajectory_driver_debug, flag_s_request)
%_________________________________________________________________
%% Documentation
%
% Authors:      Alexander Wischnewski (alexander.wischnewski@tum.de)
%
% Start Date:   31.07.2019
%
% Description:  generates error codes based on debug signals and buffers them until they are send
%               to the outside. Error codes are set bitwise, therefore they have to be decoded
%               accordingly.
%
% Inputs:
%   mloc_statemachine_debug         Debug data structure for state machine
%   mvdc_path_matching_debug        Debug data structure for path matching
%   mvdc_trajectory_driver_debug    Debug data structure for trajectory driver
%   flag_s_request                  Debug flag from TMPC controller
%
% Outputs:
%   ErrorCode                       error code integer

ErrorCode = uint32(0);

%% Fault Detection Logic
% Error Code - 001 - Performance Trajectory not ok
if(~mvdc_path_matching_debug.PerformanceTrajOK_b)
    ErrorCode = bitset(ErrorCode, 1);
end

% Error Code - 002 - Emergency Trajectory not ok
if(~mvdc_path_matching_debug.EmergencyTrajOK_b)
    ErrorCode = bitset(ErrorCode, 2);
end

% Error Code - 004 - Path Deviation too high (below critical)
if(mvdc_trajectory_driver_debug.PathMatchingStatus == TUMPathMatchingState.PM_VALID_HIGHDEVIATION)
    ErrorCode = bitset(ErrorCode, 3);
end

% Error Code - 008 - Path Deviation too high (above critical)
if(mvdc_trajectory_driver_debug.PathMatchingStatus == TUMPathMatchingState.PM_VALID_VERYHIGHDEVIATION)
    ErrorCode = bitset(ErrorCode, 4);
end

% Error Code - 016 - Trajectory Comms not ok
if(~(mvdc_trajectory_driver_debug.Trajectories_CommsStatus == TUMHealthStatus.OK))
    ErrorCode = bitset(ErrorCode, 5);
end

% Error Code - 032 - Strategy Comms not ok
if(~(mvdc_trajectory_driver_debug.Strategy_CommsStatus == TUMHealthStatus.OK))
    ErrorCode = bitset(ErrorCode, 6);
end

% Error Code - 064 - Controller not ok
if (~(mvdc_trajectory_driver_debug.ControllerStatus == TUMHealthStatus.OK))
    ErrorCode = bitset(ErrorCode, 7);
end

% Error Code - 128 - Trajectory too short
if (flag_s_request)
    ErrorCode = bitset(ErrorCode, 8);
end

% Error Code - 256 - Vehicle unstable
if (mloc_statemachine_debug.VehicleUnsafe)
    ErrorCode = bitset(ErrorCode, 9);
end

% Error Code - 512 - Software stack not OK
if (~((mvdc_trajectory_driver_debug.Strategy_Status < 50) && (mvdc_trajectory_driver_debug.Strategy_Status >= 30)))
    ErrorCode = bitset(ErrorCode, 10);
end

end

author: TUMFTM


text:
function stroud_test052 ( )

%*****************************************************************************80
%
%% TEST052 tests BALL_VOLUME_3D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 April 2009
%
%  Author:
%
%    John Burkardt
%
  n = 3;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST052\n' );
  fprintf ( 1, '  In 3 dimensions:\n' );
  fprintf ( 1, '  BALL_VOLUME_3D computes the volume of a unit ball.\n' );
  fprintf ( 1, '  BALL_VOLUME_ND will be called for comparison.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  N    R      Volume    Method\n' );
  fprintf ( 1, '\n' );

  r = 1.0;

  for i = 1 : 3

    fprintf ( 1, '  %1d  %12f  %12f  %s\n', ...
      n, r, ball_volume_3d ( r ), 'BALL_VOLUME_3D' );

    fprintf ( 1, '  %1d  %12f  %12f  %s\n', ...
      n, r, ball_volume_nd ( n, r ), 'BALL_VOLUME_ND' );

    r = r * 2.0;

  end

  return
end


text:
function out = get(f, prop)
%GET   GET method for the TRIGTECH class
%   P = GET(F, PROP) returns the property P specified in the string PROP from
%   the TRIGTECH F.  The string PROP may be the name of a TRIGTECH property (see
%   the TRIGTECH classdef files for lists) or any of the following strings:
%       'POINTS'          - Equally spaced points where F is sampled.
%       'VALUES'          - Values of F at trig points.
%       'LVAL'            - Value of F at -1.
%       'RVAL'            - Value of F at +1.
%       'TECH'            - Handle to the TRIGTECH constructor. *
%
% See also CHEBTECH, CHEBTECH2.

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DEVELOPER NOTE:
%  * Currently get(f, 'tech') returns a function handle to the tech constructor.
%    This may change in future to return instead an empty instance of the tech.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

switch prop
    case fieldnames(f)
        % Allow access to any of the properties of F via GET:
        out = f.(prop);
    case 'vscale'
        % Get the vscale:
        out = vscale(f);
    case 'points'
        % Get the trigier grid corresponding to the VALUES:
        out = f.points();
    case 'lval'
        % The value at -1:
        out = feval(f, -1);
    case 'rval'
        % The value at 1:
        out = feval(f, 1);
    case 'tech'
        % TODO: Return function handle, or empty instance of the tech?
        out = @trigtech;
    otherwise
        error('CHEBFUN:TRIGTECH:GET:proname', ...
            'Unknown property name ''%s'' for object of type trigtech.', prop);
end

end


text:
function [TrainX, TestX, TrainY, TestY, numSource, numTrain, numTarget, numTest] = CDPLSA_loadData(inputPath)
    load(inputPath);
    numSource = size(TrainData,2);                              % number of train file
    numTrain = zeros(1,numSource);
    TrainX = [];
    TrainY = [];
    for i = 1 : length(numTrain)
        numTrain(1,i) = size(TrainData{1,i},2);
        TrainX = [TrainX TrainData{1,i}];
        TrainY = [TrainY TrainLabel{1,i}];
    end

    numTarget = size(TestData,2);                              % number of train file
    numTest = zeros(1,numTarget);
    TestX = [];
    TestY = [];
    for i = 1 : length(numTest)
        numTest(1,i) = size(TestData{1,i},2);
        TestX = [TestX TestData{1,i}];
        TestY = [TestY TestLabel{1,i}];
    end
end

author: FuzhenZhuang

repo: Transfer-Learning-Toolkit


text:
function [Eft, Varft, lpyt, Eyt, Varyt] = gpla_loopred(gp, x, y, varargin)
%GPLA_LOOPRED  Leave-one-out predictions with Laplace approximation
%
%  Description
%    [EFT, VARFT, LPYT, EYT, VARYT] = GPLA_LOOPRED(GP, X, Y, OPTIONS)
%    takes a Gaussian process structure GP together with a matrix X
%    of training inputs and vector Y of training targets, and
%    evaluates the leave-one-out predictive distribution at inputs
%    X and returns means EFT and variances VARFT of latent
%    variables, the logarithm of the predictive densities PYT, and
%    the predictive means EYT and variances VARYT of observations
%    at input locations X.
%
%    OPTIONS is optional parameter-value pair
%      z      - optional observed quantity in triplet (x_i,y_i,z_i)
%               Some likelihoods may use this. For example, in case of
%               Poisson likelihood we have z_i=E_i, that is, expected value
%               for ith case.
%
%    Laplace leave-one-out is approximated in linear response style
%    by expressing the solutions for LOO problem in terms of
%    solution for the full problem. The computationally cheap
%    solution can be obtained by making the assumption that the
%    difference between these two solution is small such that their
%    difference may be treated as an Taylor expansion truncated at
%    first order (Winther et al 2012, in progress).
%
%  See also
%    GP_LOOPRED, GP_PRED
%
% Copyright (c) 2011-2012  Aki Vehtari, Ville Tolvanen

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

  ip=inputParser;
  ip.FunctionName = 'GPLA_LOOPRED';
  ip.addRequired('gp', @(x) isstruct(x));
  ip.addRequired('x', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
  ip.addRequired('y', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
  ip.addParamValue('z', [], @(x) isreal(x) && all(isfinite(x(:))))
  ip.addParamValue('method', 'cavity', @(x) ismember(x, {'lrs' 'cavity' 'inla'}))
  ip.addParamValue('fcorr', 'off', @(x) ismember(x, {'off', 'fact', 'cm2', 'on'}));
  ip.parse(gp, x, y, varargin{:});
  z=ip.Results.z;
  method = ip.Results.method;
  fcorr=ip.Results.fcorr;
  [tn,nin] = size(x);

  switch method

    case 'lrs'
      % Manfred Opper and Ole Winther (2000). Gaussian Processes for
      % Classification: Mean-Field Algorithms. In Neural
      % Computation, 12(11):2655-2684.
      %
      % Ole Winther et al (2012). Work in progress.

      % latent posterior
      [f, sigm2ii] = gpla_pred(gp, x, y, 'z', z, 'tstind', [], 'fcorr', fcorr);
      if isnan(f)
        Eft=NaN; Varft=NaN; lpyt=NaN; Eyt=NaN; Varyt=NaN;
        return
      end

      deriv = gp.lik.fh.llg(gp.lik, y, f, 'latent', z);
      La = 1./-gp.lik.fh.llg2(gp.lik, y, f, 'latent', z);
      % really large values don't contribute, but make variance
      % computation unstable. 2e15 approx 1/(2*eps)
      La = min(La,2e15);

      switch gp.type
        case 'FULL'
          % FULL GP (and compact support GP)
          K = gp_trcov(gp,x);
          Varft=1./diag(inv(K+diag(La)))-La;

        case 'FIC'
          % FIC
          % Use inverse lemma for FIC low rank covariance matrix approximation
          % Code adapated from gp_pred

          u = gp.X_u;
          m = size(u,1);
          % Turn the inducing vector on right direction
          if size(u,2) ~= size(x,2)
            u=u';
          end
          [Kv_ff, Cv_ff] = gp_trvar(gp, x);  % 1 x f  vector
          K_fu = gp_cov(gp, x, u);   % f x u
          K_uu = gp_trcov(gp, u);     % u x u, noiseles covariance K_uu
          Luu = chol(K_uu,'lower');
          B=Luu\(K_fu');
          Qv_ff=sum(B.^2)';
          % Add also La to the vector of diagonal elements
          Lav = Cv_ff-Qv_ff + La;   % 1 x f, Vector of diagonal elements

          % iLaKfu = diag(inv(Lav))*K_fu = inv(La)*K_fu
          iLaKfu = zeros(size(K_fu));  % f x u,
          n=size(x,1);
          for i=1:n
            iLaKfu(i,:) = K_fu(i,:)./Lav(i);  % f x u
          end
          A = K_uu+K_fu'*iLaKfu;
          A = (A+A')./2;
          L = iLaKfu/chol(A);

          %Varft=1./diag(inv(K+diag(La)))-La;
          Varft=1./(1./Lav - sum(L.^2,2))-La;

        case {'PIC' 'PIC_BLOCK'}
          % PIC
          % Use inverse lemma for PIC low rank covariance matrix approximation
          % Code adapated from gp_pred (here Lab is same La in gp_pred)

          u = gp.X_u;
          ind = gp.tr_index;
          if size(u,2) ~= size(x,2)
            % Turn the inducing vector on right direction
            u=u';
          end

          % Calculate some help matrices
          [Kv_ff, Cv_ff] = gp_trvar(gp, x);  % 1 x f  vector
          K_fu = gp_cov(gp, x, u);         % f x u
          K_uu = gp_trcov(gp, u);    % u x u, noiseles covariance K_uu
          Luu = chol(K_uu)';

          % Evaluate the Lambda (La) for specific model
          % Q_ff = K_fu*inv(K_uu)*K_fu'
          % Here we need only the diag(Q_ff), which is evaluated below
          B=Luu\K_fu';
          iLaKfu = zeros(size(K_fu));  % f x u
          for i=1:length(ind)
            Qbl_ff = B(:,ind{i})'*B(:,ind{i});
            [Kbl_ff, Cbl_ff] = gp_trcov(gp, x(ind{i},:));
            % Add also La to the diagonal
            Lab{i} = Cbl_ff - Qbl_ff + diag(La(ind{i}));
            iLaKfu(ind{i},:) = Lab{i}\K_fu(ind{i},:);
          end
          A = K_uu+K_fu'*iLaKfu;
          A = (A+A')./2;            % Ensure symmetry
          L = iLaKfu/chol(A);

          % From this on evaluate the prediction
          % See Snelson and Ghahramani (2007) for details
          n=size(y,1);
          iCv=zeros(n,1);
          for i=1:length(ind)
            iCv(ind{i},:) = diag(inv(Lab{i}));
          end

          %Varft=1./diag(inv(K+diag(La)))-La;
          Varft=1./(iCv - sum(L.^2,2))-La;

        case 'CS+FIC'
          % CS+FIC
          % Use inverse lemma for CS+FIC
          % Code adapated from gp_pred (Here Las is same as La in gp_pred)

          u = gp.X_u;
          if size(u,2) ~= size(x,2)
            % Turn the inducing vector on right direction
            u=u';
          end

          n = size(x,1);
          m = size(u,1);
          ncf = length(gp.cf);

          % Indexes to all non-compact support and compact support covariances.
          cf1 = [];
          cf2 = [];

          % Loop through all covariance functions
          for i1 = 1:ncf
            if ~isfield(gp.cf{i1},'cs')
              % Non-CS covariances
              cf1 = [cf1 i1];
            else
              % CS-covariances
              cf2 = [cf2 i1];
            end
          end

          % First evaluate needed covariance matrices
          % v defines that parameter is a vector
          [Kv_ff, Cv_ff] = gp_trvar(gp, x, cf1); % f x 1  vector
          K_fu = gp_cov(gp, x, u, cf1);          % f x u
          K_uu = gp_trcov(gp, u, cf1);    % u x u, noiseles covariance K_uu
          K_uu = (K_uu+K_uu')./2;         % ensure the symmetry of K_uu

          Luu  = chol(K_uu)';
          % Evaluate the Lambda (La)
          % Q_ff = K_fu*inv(K_uu)*K_fu'
          B=Luu\(K_fu');       % u x f
          Qv_ff=sum(B.^2)';
          % Add also La to the vector of diagonal elements
          Lav = Cv_ff-Qv_ff + La;   % f x 1, Vector of diagonal elements

          K_cs = gp_trcov(gp,x,cf2);
          Las = sparse(1:n,1:n,Lav,n,n) + K_cs;

          iLaKfu = Las\K_fu;
          A = K_uu+K_fu'*iLaKfu;
          A = (A+A')./2;     % Ensure symmetry
          L = iLaKfu/chol(A);

          %Varft=1./diag(inv(K+diag(La)))-La;
          Varft=1./(diag(inv(Las)) - sum(L.^2,2))-La;

        otherwise
          error('Unknown type of Gaussian process')
      end

      % check if cavity variances are negative
      ii=find(Varft<0);
      if ~isempty(ii)
        warning('gpla_loopred: some LOO latent variances are negative');
        Varft(ii) = gp.jitterSigma2;
      end
      Eft=f-Varft.*deriv;

      if nargout==3
        lpyt = gp.lik.fh.predy(gp.lik, Eft, Varft, y, z);
      elseif nargout>3
        [lpyt,Eyt,Varyt] = gp.lik.fh.predy(gp.lik, Eft, Varft, y, z);
      end

    case 'cavity'
      % using EP equations

      % latent posterior
      [f, sigm2ii] = gpla_pred(gp, x, y, 'z', z, 'tstind', [], 'fcorr', fcorr);
      % remove the jitter variance (see sigm2_t below)
      sigm2ii  = sigm2ii-gp.jitterSigma2;

      % "site parameters"
      W        = -gp.lik.fh.llg2(gp.lik, y, f, 'latent', z);
      deriv    = gp.lik.fh.llg(gp.lik, y, f, 'latent', z);
      % add the jitter variance (see sigm2ii above)
      sigm2_t  = 1./W+gp.jitterSigma2;
      mu_t     = f + sigm2_t.*deriv;

      % "cavity parameters"
      sigma2_i = 1./(1./sigm2ii-1./sigm2_t);
      myy_i    = sigma2_i.*(f./sigm2ii-mu_t./sigm2_t);
      % check if cavity varianes are negative
      ii=find(sigma2_i<0);
      if ~isempty(ii)
        warning('gpla_loopred: some cavity variances are negative');
        sigma2_i(ii) = sigm2ii(ii);
        myy_i(ii) = f(ii);
      end

      % leave-one-out predictions
      Eft=myy_i;
      Varft=sigma2_i;

      if nargout==3
        lpyt = gp.lik.fh.predy(gp.lik, Eft, Varft, y, z);
      elseif nargout>3
        [lpyt,Eyt,Varyt] = gp.lik.fh.predy(gp.lik, Eft, Varft, y, z);
      end

    case 'inla'
      % Leonhard Held and Birgit Schr�dle and H�vard Rue (2010)
      % Posterior and Cross-validatory Predictive Checks: A
      % Comparison of MCMC and INLA. In (eds) Thomas Kneib and
      % Gerhard Tutz, Statistical Modelling and Regression
      % Structures, pp. 91-110. Springer.

      % latent posterior
      [f, sigm2ii, lp] = gpla_pred(gp, x, y, 'z', z, 'tstind', [], 'fcorr', fcorr);

      Eft = zeros(tn,1);
      Varft = zeros(tn,1);
      lpyt = zeros(tn,1);
      minf = f-6.*sqrt(sigm2ii);
      maxf = f+6.*sqrt(sigm2ii);
      for i=1:tn
        if isempty(z)
          z1 = [];
        else
          z1 = z(i);
        end
        [m0, m1, m2] = quad_moments(@(x) norm_pdf(x, f(i), sqrt(sigm2ii(i)))./llvec(gp.lik,y(i),x,z1), minf(i), maxf(i));
        Eft(i) = m1;
        Varft(i) = m2-Eft(i)^2;
        lpyt(i) = -log(m0);
      end

      if nargout>3
        [tmp,Eyt,Varyt] = gp.lik.fh.predy(gp.lik, Eft, Varft, y, z);
      end
      if sum((abs(lpyt)./abs(lp) > 5) == 1) > 0.1*tn;
        warning('Very bad predictive densities, gpla_loopred might not be reliable, check results!');
      end

  end

end

function expll = llvec(gplik, y, f, z)
  for i=1:size(f,2)
    expll(i) = exp(gplik.fh.ll(gplik, y, f(i), z));
  end
end

author: gpstuff-dev

repo: gpstuff

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


text:
function combo_test09 ( )

%*****************************************************************************80
%
%% COMBO_TEST09 tests I4VEC_SEARCH_BINARY_A and I4VEC_SORT_INSERT_A.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 January 2011
%
%  Author:
%
%    John Burkardt
%
  n = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'COMBO_TEST09\n' );
  fprintf ( 1, '  Integer vectors.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  I4VEC_SORT_INSERT_A ascending sorts;\n' );
  fprintf ( 1, '  I4VEC_SEARCH_BINARY_A searches a ascending sorted vector.\n' );

  a(1:n) = [ 6, 7, 1, 0, 4, 3, 2, 1, 5, 8 ]';

  i4vec_print ( n, a, '  Before ascending sort:' );

  a = i4vec_sort_insert_a ( n, a );

  i4vec_print ( n, a, '  After ascending sort:' );

  b = 5;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Now search for an instance of the value %d\n', b );

  index = i4vec_search_binary_a ( n, a, b );

  fprintf ( 1, '\n' );
  if ( index == 0 )
    fprintf ( 1, '  The value does not occur.\n' );
  else
    fprintf ( 1, '  The value occurs at index = %d\n', index );
  end

  return
end


text:
function [tmptypeeeg, tmprt, tmpresponseeeg, n] = loaddat( FILENAME)

% eeg_load_scan_dat - Read Neuroscan .dat format
%
% usage  [typeeeg, rt, response, n] = loaddat( FILENAME)
%
% Inputs:
%   FILENAME     input Neuroscan .dat file
% Outputs:
% 	typeeeg		type of the sweeps
%	rt		reaction time of the subject
%	response	response of the subject
%	n		number of sweeps
%

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

% Licence:  GNU GPL, no implied or express warranty
% History:  05/2001, arno_delorme@salk.edu
%           - Version 1.0
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin<1
	fprintf('Not enought arguments\n');
	help loaddat
	return;
end;

BOOL='int16';
ULONG='int32';
FLOAT='float32';

fid=fopen(FILENAME,'r','ieee-le');
if fid<0
	fprintf(2,['Error LOADEEG: File ' FILENAME ' not found\n']);
	return;
end;

% skip the first 20 lines
% -----------------------
for index=1:20	fgetl(fid); end;

% read the matrix
% ---------------
tmpMat 		    = fscanf(fid, '%f', [5, inf]);
tmptypeeeg  	= tmpMat(3,:);
tmpresponseeeg  = tmpMat(4,:);
tmprt       	= tmpMat(5,:) * 1000;
n        	    = size( tmpMat, 2);

fclose(fid);
return;

author: PatternRecognition

repo: OpenBMI


text:
% EEG_EPOCH2CONTINUOUS - convert epoched dataset to continuous dataset
%                          with data epochs separated by boundary events.
% Usage:
%           >> EEGOUT = eeg_epoch2continuous(EEGIN);
%
% Inputs:
%   EEGIN  - a loaded epoched EEG dataset structure.
%
% Inputs:
%   EEGOUT - a continuous EEG dataset structure.
%
% Authors: Arnaud Delorme, SCCN, INC, UCSD, January, 2012

% Copyright (C) Arnaud Delorme, SCCN, INC, UCSD, October 11, 2012, arno@sccn.ucsd.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function EEG = eeg_epoch2continuous(EEG)

if nargin < 1
    help eeg_epoch2continuous;
    return;
end

EEG.data = reshape(EEG.data, size(EEG.data,1), size(EEG.data,2)*size(EEG.data,3));

eeglab_options;
for index = 1:EEG.trials-1
    EEG.event(end+1).type = eeg_boundarytype(EEG);
    EEG.event(end  ).latency  = index*EEG.pnts-0.5;
    EEG.event(end  ).duration = NaN;
end

EEG.pnts   = size(EEG.data,2);
EEG.trials = 1;
if ~isempty(EEG.event) && isfield(EEG.event, 'epoch')
    EEG.event = rmfield(EEG.event, 'epoch');
end

EEG = eeg_checkset(EEG);

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [GX,dGdX,dGdP] = g_embedAR(X,P,ut,in)
% AR(1) evolution function

n = size(X,1);
x = X(1:n/2);
z = X(n/2+1:n);

[gx,J,dgdP] = VBA_evalFun('g',x,P,ut,in.opt,in.dim,0);

GX = gx;
dGdX = [J;zeros(n/2)];
dGdP = dgdP;

author: MBB-team

repo: VBA-toolbox


text:
% fixedtrx2mtrax(inmatname,outmatname)
% inmatname: name of a matfile created by fixerrorsgui
% outmatname: name of a matfile to save results to
% results will be saved in the format returned by mtrax
function fixedtrx2mtrax(inmatname,outmatname)

load(inmatname);
nflies = length(trx);
nframes = 0;
for fly = 1:nflies,
  nframes = max(nframes,trx(fly).endframe);
end
x_pos = nan(nflies,nframes);
y_pos = nan(nflies,nframes);
angle = nan(nflies,nframes);
maj_ax = nan(nflies,nframes);
min_ax = nan(nflies,nframes);
identity = repmat((1:nflies)',[1,nframes]);
for fly = 1:nflies,
  x_pos(fly,trx(fly).firstframe:trx(fly).endframe) = trx(fly).x;
  y_pos(fly,trx(fly).firstframe:trx(fly).endframe) = trx(fly).y;
  angle(fly,trx(fly).firstframe:trx(fly).endframe) = trx(fly).theta;
  maj_ax(fly,trx(fly).firstframe:trx(fly).endframe) = trx(fly).a;
  min_ax(fly,trx(fly).firstframe:trx(fly).endframe) = trx(fly).b;
end
keep = ~isnan(x_pos);
x_pos = x_pos(keep)'-1;
y_pos = y_pos(keep)'-1;
angle = angle(keep)';
maj_ax = maj_ax(keep)';
min_ax = min_ax(keep)';
identity = identity(keep)' - 1;
ntargets = sum(double(keep),1);

save(outmatname,'x_pos','y_pos','angle','maj_ax','min_ax','identity','ntargets');

author: kristinbranson

repo: JAABA


text:
function get_voc12train_truncation_stats(cls, visualize, debug)

if nargin < 2
    visualize = 0;
    debug = 0;
end
if nargin < 3
    debug = 0;
end

setup_path;

% LOAD CAD MODELS
CADPathName = fullfile(PASCAL3D_DIR, 'CAD', [cls '.mat']);
cad = load(CADPathName);
cad = cad.(cls);

%image_dir = '/orions3-zfs/projects/rqi/Dataset/VOC12/VOCdevkit/VOC2012/JPEGImages/';
image_dir = fullfile(PASCAL3D_DIR, 'Images', [cls '_pascal']);
ids = textread(sprintf(VOCopts.imgsetpath, 'train'), '%s');

% LOGGINGS
azimuths = [];
elevations = [];
thetas = [];
distances = [];
gt_bboxes = [];
cad_bboxes = [];
gt_bboxes_draw  = [];
cad_bboxes_draw = [];
area_ratios = []; % gt bbox / cad bbox

for i = 1:numel(ids)
    if mod(i,100) == 0
        fprintf('%s: %d/%d\n', cls, i, numel(ids));
    end

    % Load annotation for the i'th image
    filename = fullfile(PASCAL3D_DIR, 'Annotations', sprintf('%s_pascal/%s.mat', cls, ids{i}));
    if ~exist(filename,'file')
        continue;
    end
    anno = load(filename);
    objects = anno.record.objects;
    for k = 1:length(objects)
        obj = objects(k);
        % Proceed only when the object is of cls
        if strcmp(obj.class,cls)
            bbox = obj.bbox;
            bbox_draw = [bbox(1) bbox(2) bbox(3)-bbox(1) bbox(4)-bbox(2)];
            im = imread(fullfile(image_dir, [ids{i} '.jpg']));
            % Show original image
            if debug
                imshow(im);
                hold on;
                rectangle('Position', bbox_draw, 'EdgeColor', 'g');
            end
            % Show annotated anchor points
            if debug && isfield(obj, 'anchors') == 1 && isempty(obj.anchors) == 0
                names = fieldnames(obj.anchors);
                for j = 1:numel(names)
                    if obj.anchors.(names{j}).status == 1
                        if isempty(obj.anchors.(names{j}).location) == 0
                            x = obj.anchors.(names{j}).location(1);
                            y = obj.anchors.(names{j}).location(2);
                            plot(x, y, 'ro');
                        else
                            fprintf('anchor point %s is missing!\n', names{j});
                        end
                    end
                end
            end
            % show overlay of the CAD model
            if isfield(obj, 'cad_index') == 1 && isempty(obj.cad_index) == 0 && ...
                    isfield(obj, 'viewpoint') == 1 && isfield(obj.viewpoint, 'azimuth') == 1
                vertices = cad(obj.cad_index).vertices;
                faces = cad(obj.cad_index).faces;
                x2d = project_3d(vertices, obj);
                if isempty(x2d) == 0

                    if debug
                        patch('vertices', x2d ,'faces', faces, ...
                        'FaceColor', 'blue', 'FaceAlpha', 0.2, 'EdgeColor', 'none');
                    end

                    % Find bbox for the CAD model
                    bbox_cad = [min(x2d(:,1)),min(x2d(:,2)),max(x2d(:,1)),max(x2d(:,2))];
                    bbox_cad_draw = [bbox_cad(1), bbox_cad(2), bbox_cad(3)-bbox_cad(1), bbox_cad(4)-bbox_cad(2)];

                    area_ratio = bbox_draw(3)*bbox_draw(4)/(bbox_cad_draw(3)*bbox_cad_draw(4));
                    if debug
                        rectangle('Position', bbox_cad_draw, 'EdgeColor','r');
                        fprintf('Approx. visible ratio: %f\n', area_ratio);
                        hold off;
                    end

                    azimuths = [azimuths obj.viewpoint.azimuth];
                    area_ratios = [area_ratios area_ratio];
                    gt_bboxes = [gt_bboxes; bbox];
                    cad_bboxes = [cad_bboxes; bbox_cad];
                    gt_bboxes_draw = [gt_bboxes_draw; bbox_draw];
                    cad_bboxes_draw = [cad_bboxes_draw; bbox_cad_draw];
                end

            end
            if debug
                fprintf('press button to continue...\n');
                pause();
                close all;
            end
        end
    end
end

widths = cad_bboxes(:,3) - cad_bboxes(:,1);
heights = cad_bboxes(:,4) - cad_bboxes(:,2);
r = zeros(length(azimuths), 4);
gt_cad_bboxes_diff = gt_bboxes - cad_bboxes;
r(:,1) = gt_cad_bboxes_diff(:,1) ./ widths;
r(:,3) = gt_cad_bboxes_diff(:,3) ./ widths;
r(:,2) = gt_cad_bboxes_diff(:,2) ./ heights;
r(:,4) = gt_cad_bboxes_diff(:,4) ./ heights;
left = r(:,1); top = r(:,2);
right = r(:,3); bottom = r(:,4);
mkdir(g_truncation_statistics_folder);
save(fullfile(g_truncation_statistics_folder, [cls '_truncation_stats']), 'left', 'top', 'right', 'bottom');

if visualize
    % Show relation between azimuth and area_ratios
    vnum = 8;
    azimuth_interval = [0 (360/(vnum*2)):(360/vnum):360-(360/(vnum*2))];
    ratios_by_view = zeros(vnum,1);
    for i = 1:length(azimuths)
        azimuths(i) = find_interval(azimuths(i), azimuth_interval);
    end
    for v = 1:vnum
        [tmp, index] = find(azimuths == v);
        ratios_by_view(v) = median(area_ratios(index));
    end
    figure, subplot(1,2,1), bar(ratios_by_view);
    subplot(1,2,2), scatter(azimuths, area_ratios);

    % Show hist of ratios of gt_cad_bboxes_diff/[width,height]
    figure, subplot(2,2,1), hist(left), title('left');
    subplot(2,2,2), hist(top), title('right');
    subplot(2,2,3), hist(right), title('top');
    subplot(2,2,4), hist(bottom), title('bottom');
end

function ind = find_interval(azimuth, a)
for i = 1:numel(a)
    if azimuth < a(i)
        break;
    end
end
ind = i;
if azimuth > a(end)
    ind = 1;
end

author: ShapeNet

repo: RenderForCNN


text:
function [A,U,V] = regutm(m,n,s)
%REGUTM Test matrix for regularization methods.
%
% [A,U,V] = regutm(m,n,s)
%
% Generates a random m-times-n matrix A such that A*A' and A'*A
% are oscillating.  Hence, in the SVD of A,
%    A = U*diag(s)*V',
% the number of sign changes in U(:,i) and V(:,i) is exactly i-1.
%
% The third argument s specifies the singular values of A.  If not
% present, then s = logspace(0,round(log10(eps)),n).

% Reference: P. C. Hansen, "Test matrices for regularization methods",
% SIAM J. Sci. Comput. 16 (1995), 506--512.

% Per Christian Hansen, IMM, 07/30/97.

% Initialization.
if (nargin==1), n = m; end
if (nargin<3), s = logspace(0,round(log10(eps)),min(m,n)); end

% Special treatment of the case m < n.
if (m < n), [A,V,U] = regutm(n,m,s); A = A'; return, end

% Generate random bidiagonal matrix with nonnegative elements.
if (n < 100), mu = .222*n + .0278*n^2; else mu = 3*n; end
B = abs(diag(randn(n,1)+mu) + diag(randn(n-1,1)+mu,1));

% Compute the SVD of B.
[U,dummy,V] = svd(B); clear dummy

% Repeat if m > n.
if (m > n)
  clear U
  B = abs(diag(randn(m,1)+mu) + diag(randn(m-1,1)+mu,1));
  [U,dummy] = svd(B); clear dummy, U = U(:,1:n);
end

% Compute A.
A = U*diag(s)*V';

repo: REST


text:
function test_suite = test_createRhombododecahedron
%TESTCREATERHOMBODODECAHEDRON  One-line description here, please.
%
%   output = testCreateRhombododecahedron(input)
%
%   Example
%   testCreateRhombododecahedron
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2010-12-07,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

test_suite = functiontests(localfunctions);

function testCreation(testCase) %#ok<*DEFNU>

createRhombododecahedron();

function testVEFCreation(testCase)

[v, e, f] = createRhombododecahedron();
testCase.assertTrue(~isempty(v));
testCase.assertTrue(~isempty(e));
testCase.assertTrue(~isempty(f));

[nv, ne, nf] = getMeshElementsNumber;
testCase.assertEqual([nv 3], size(v));
testCase.assertEqual([ne 2], size(e));
testCase.assertEqual(nf, length(f));

function testVFCreation(testCase)

[v, f] = createRhombododecahedron();
testCase.assertTrue(~isempty(v));
testCase.assertTrue(~isempty(f));

[nv, ne, nf] = getMeshElementsNumber; %#ok<ASGLU>
testCase.assertEqual([nv 3], size(v));
testCase.assertEqual(nf, length(f));

function testMeshCreation(testCase)

mesh = createRhombododecahedron();
testCase.assertTrue(isstruct(mesh));
testCase.assertTrue(isfield(mesh, 'vertices'));
testCase.assertTrue(isfield(mesh, 'edges'));
testCase.assertTrue(isfield(mesh, 'faces'));

[nv, ne, nf] = getMeshElementsNumber;
testCase.assertEqual([nv 3], size(mesh.vertices));
testCase.assertEqual([ne 2], size(mesh.edges));
testCase.assertEqual(nf, length(mesh.faces));

function testFacesOutwards(testCase)

[v, e, f] = createRhombododecahedron(); %#ok<ASGLU>

centro = centroid(v);
fc  = meshFaceCentroids(v, f);
fc2 = createVector(centro, fc);
n   = meshFaceNormals(v, f);

testCase.assertEqual(size(n), size(fc2));

dp = dot(fc2, n, 2);

testCase.assertTrue(sum(dp <= 0) == 0);

function [nv, ne, nf] = getMeshElementsNumber

nv = 14;
ne = 24;
nf = 12;

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [ x, seed ] = c_normal_01_sample ( seed )

%*****************************************************************************80
%
%% C_NORMAL_01_SAMPLE samples the complex Normal 01 PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, complex X, a sample of the PDF.
%
  [ v1, seed ] = r8_uniform_01 ( seed );
  [ v2, seed ] = r8_uniform_01 ( seed );

  x_r = sqrt ( - 2.0 * log ( v1 ) ) * cos ( 2.0 * pi * v2 );
  x_c = sqrt ( - 2.0 * log ( v1 ) ) * sin ( 2.0 * pi * v2 );

  x = x_r + i * x_c;

  return
end