text:
function signal = set_infer_markers(varargin)
% Infers markers from an event channel, if possible.
% Signal = set_infer_markers(Signal)
%
% If the data set has one (or more) channels with relatively few unique values, and relatively few
% off-median values (or plateaus thereof), markers are being generated to encode whenever such a
% channel changes its value. By default, the detected event channels are being removed from the
% data.
%
% This funciton is automatically called by io_loadset if no markers are present in the given
% set. io_loadset has a parameter ('markerchannel') which can be utilized to customize the settings
% below (as a cell array of arguments to set_infer_markers).
%
% In:
%   Signal                    : continuous EEGLAB data set from which an interval should be selected
%
%   MaxEvents                 : a channel is only considered for event generation if it would produce
%                               at most this many distinct events (default: 25000)
%
%   MaxTypes                  : a channel is only considered for event generation if it would produce
%                               at most this many distinct event types (default 300)
%
%   IncludeLabel              : when to include the channel label into the generated event types;
%                               can be 'always', 'multiplechans', 'never' (default: 'multiplechans')
%
%   AllPositive               : whether all events must have positive offsets from the baseline
%                               (median) value (default: true)
%
%   EncodePlateaus            : Whether the trigger channel may contain (non-baseline) plateaus, to
%                               be encoded in event duration. (default: true)
%
%   StrictInteger             : only consider channels that would produce integer event types
%                               (otherwise: round to integer) (default: true)
%
%   RelativeBaselineThreshold : Relative baseline threshold. If the event types would be larger than
%                               this threshold, they will be re-coded to be relative to the smallest
%                               event type in the channel. (default: 10000)
%
%   OmitBaselineThreshold     : Baseline omission threshold. If the baseline event takes up more than
%                               this fraction of the data, it will not be encoded as events.
%                               (default: 0.8)
%
%   RemoveEventChannels       : Whether to remove channels that have been identified as event
%                               channels (default: true)
%
% Out:
%   Signal    : data set restricted to the selected range
%
% Examples:
%   % for a data set with one or more trigger channels, fill in contents of the .event field
%   % and remove the trigger channel(s); do nothing if the data set does not contain trigger channels
%   eeg = set_infer_markers(eeg)
%
%   % fill in contents of the .event field from trigger channels, and do not remove these channels
%   eeg = set_infer_markers('Signal',eeg,'RemoveEventChannels',false)
%
%   % if an event channel is not being detected because it produces more than the default MaxEvents
%   % number of events, use a larger cutoff value
%   eeg = set_infer_markers('Signal',eeg,'MaxEvents',30000)
%
%   % if an event channel is not being detected because it produces more than the default MaxTypes
%   % number of distinct event types, use a larger cutoff value
%   eeg = set_infer_markers('Signal',eeg,'MaxTypes',1000)
%
%   % for event channels that have values that far from zero (say, 30000+), the values are by default
%   % made relative to to the baseline value (most frequent or smallest value, depending on the
%   % EncodePlateaus setting) -- if this is not intended, override the cutoff with a larger value
%   eeg = set_infer_markers('Signal',eeg,'RelativeBaselineThreshold',Inf)
%
%   % if a trigger channel contains many subsequent events of the same type (i.e. plateaus), by
%   % default only one event is generated with the appropriate duration; if instead an event should
%   % be generated for each sample, override that default (and possibly increase MaxEvents)
%   eeg = set_infer_markers('Signal',eeg,'EncodePlateaus',false,'MaxEvents',Inf)
%
% See also:
%   io_loadset
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-01-26

declare_properties('independent_channels',false,'independent_trials',false);

arg_define(varargin, ...
    arg_norep({'signal','Signal'}), ...
    arg({'force_processing','ForceProcessing'},false,[], 'Force marker processing. Scan for marker channels even if a data set already has events.'), ...
    arg({'max_events','MaxEvents'},30000,uint32([0 1000000000]), 'Number of events allowed. If a channel would produce more than this many events, it is not considered an event channel.'), ...
    arg({'max_event_fraction','MaxEventFraction'},0.3,[0 1], 'Maximum fraction of events. If a channel would produce more events than this fraction of its total number of samples, the channel is not considered an event channel.'), ...
    arg({'max_types','MaxTypes'},300,uint32([0 1000000000]), 'Number of eventtypes allowed. If a channel would produce more than this many event types, it is not considered an event channel.'), ...
    arg({'include_label','IncludeLabel'},'multiplechans',{'always','multiplechans','never'},'Integrate channel labels. When to integrate channel labels into the event types - by default. only when there are multiple event channels.'), ...
    arg({'all_positive','AllPositive'},true,[],'Events must be positive. Only if all event types produced by a channel would correspond to positive numbers, the channel is considered an event channel.'), ...
    arg({'encode_plateaus','EncodePlateaus'},true,[],'Plateaus encoded. The trigger channel may contain (non-baseline) plateaus, to be encoded in event duration.'), ...
    arg({'strict_integer','StrictInteger'},true,[],'Strictly integer. Only allow channels that have integer event types (otherwise: round to integers).'), ...
    arg({'relative_baseline_thresh','RelativeBaselineThreshold'},10000,[],'Relative baseline threshold. If the event types would be larger than this threshold, they will be re-coded to be relative to the smallest event type in the channel.'), ...
    arg({'omit_baseline_thresh','OmitBaselineThreshold'},0.8,[0 1],'Baseline omission threshold. If the baseline event takes up more than this fraction of the data, it will not be encoded as events.'), ...
    arg({'remove_eventchns','RemoveEventChannels'},true,[],'Remove event channels. Whether to remove those channels that have been identified as event channels.'));

% obtain the chanlocs
if isfield(signal,{'head','parts'})
    signal = exp_eval(signal); end
utl_check_fields(signal,{'data','chanlocs'},'signal','signal');

% figure out which are the event channels
fprintf('Scanning potential marker channel ');
for k=1:size(signal.data,1) %#ok<*NODEF>
    fprintf('%i ',k);
    X = signal.data(k,:);
    X(isnan(X)) = 0;
    % do a few vectorized computations for speed...
    numtypes(k) = length(unique(X));
    if numtypes(k) <= max_types
        if encode_plateaus
            % in the plateaus case, the baseline is 0 or the smallest value
            baseline(k) = min(X);
            if baseline(k) < relative_baseline_thresh && baseline(k) > 0
                baseline(k) = 0; end
            % whether baseline periods are encoded as events depends on the fraction of time in the
            % data set
            omit_baseline(k) = mean(X == baseline(k)) > omit_baseline_thresh;
            eventmask = ([1; diff(X(:))] ~= 0);
            if omit_baseline(k)
                eventmask = eventmask & (X(:) ~= baseline(k)); end
            numevents(k) = nnz(eventmask);
        else
            % without plateaus, the baseline type is the mode
            baseline(k) = mode(X);
            numevents(k) = nnz(X - baseline(k));
            omit_baseline(k) = true;
        end
    else
        % channel will not be considered
        baseline(k) = 0;
        omit_baseline(k) = true;
        numevents(k) = 0;
    end
    allpositive(k) = all(X >= baseline(k));
    allinteger(k) = all(abs(X-round(X)) <= eps(X));
end
fprintf('\n');

% find the mask of event channels
eventchans = (numtypes <= max_types) & (numevents <= max_events) & (numevents <= size(signal.data,2)*max_event_fraction) & (~all_positive | allpositive) & (~strict_integer | allinteger);

% determine whether to include labels
include_label = hlp_rewrite(include_label,'always',true,'never',false,'multiplechans',nnz(eventchans) > 1);

% initialize events
if ~isfield(signal,'event') || isempty(signal.event)
    signal.event = struct('latency',{},'duration',{},'type',{}); end

% generate events
for k = find(eventchans)
    try
        % find event codes per sample
        X = round(signal.data(k,:) - baseline(k));
        X(isnan(X)) = 0;
        % do a run-length encode into events
        lat = find([1; diff(X(:))] ~= 0);
        % optionally mask out baseline events
        if omit_baseline(k)
            mask = X(lat) ~= 0;
        else
            mask = true(1,length(lat));
        end
        % get values, duration, and latency as cell arrays
        dur = num2cell(diff([lat; numel(X)+1]));
        val = cellfun(@num2str,num2cell(X(lat(mask)),1),'UniformOutput',false);
        if include_label
            val = cellfun(@(n) [signal.chanlocs(k).labels '_' n],val,'UniformOutput',false); end
        lat = num2cell(lat(mask));
        dur = dur(mask);
        if ~isempty(lat)
            % make space
            signal.event(end+length(lat)).latency = [];
            % insert new content
            [signal.event(end-length(lat)+1:end).latency] = lat{:};
            [signal.event(end-length(lat)+1:end).duration] = dur{:};
            [signal.event(end-length(lat)+1:end).type] = val{:};
        else
            % the channel is not actually an event channel (but a flatline)
            eventchans(k) = false;
        end
    catch e
        disp(['Could not process channel ' num2str(k)]);
        env_handleerror(e);
    end
end

% sort events
signal.event = signal.event(hlp_getresult(2,@sort,[signal.event.latency]));

% remove event channels
signal.data = signal.data(~eventchans,:,:,:,:,:,:,:);
signal.chanlocs = signal.chanlocs(~eventchans);
signal.nbchan = size(signal.data,1);

repo: REST


text:
function exportTsvFile(data, filename, withQuotes)
% exportTsvFile
%
%   Export structure or table to a tsv (or txt) file. String (char) entries
%   will be enclosed in double quotes ("") by default, whereas numeric
%   (double) entries will not.
%
% Input:
%
%   data       A structure or table containing the annotation information
%              associated with the GEM reactions, metabolites, or genes.
%
%   filename   Name of the file to which the data will be written,
%              delimited by tabs.
%              If filename does not end in ".tsv" or ".txt", the extension
%              will automatically be appended as ".tsv".
%
%   withQuotes Enclose string elements with double quotes (opt, default is
%              TRUE)
%
%
% Usage:
%
%   exportTsvFile(data, filename, withQuotes);
%

if nargin < 3
    withQuotes = true;
end

if isstruct(data)
    data = struct2table(data);
end
if ~istable(data)
    error('Invalid DATA format. Must be a STRUCT or a TABLE.');
end
if ~endsWith(lower(filename), '.tsv') && ~endsWith(lower(filename), '.txt')
    filename = strcat(filename, '.tsv');
end

writetable(data, filename, 'Delimiter', '\t', 'QuoteStrings', withQuotes, 'FileType', 'text');

author: SysBioChalmers

repo: Human-GEM


text:
function [d2f, err] = tapas_riddersdiff2(f, x, varargin)
% Differentiates the function f *twice* at point x according to Ridders' method:
%
% Ridders, CJF. (1982). Accurate computation of F'(x) and F'(x) F''(x). Advances in Engineering
%     Software, 4(2), 75-6.
%
% INPUT:
%    f             Function handle of a scalar real function of one real variable
%    x             Point at which to differentiate f
%
% OUTPUT:
%    d2f           Second derivative of f at x
%    err           Error estimate
%
% OPTIONS:
%    Optionally, the third argument of the function can be a structure containing further
%    settings for Ridder's method.
%
%    varargin{1}.init_h      Initial finite difference (default: 1)
%    varargin{1}.div         Divisor used to reduce h on each step (default: 1.2)
%    varargin{1}.min_steps   Minimum number of steps in h (default: 3)
%    varargin{1}.max_steps   Maximum number of steps in h (default: 100)
%    varargin{1}.tf          Terminate if last step worse than preceding by a factor of tf
%                            (default: 2)
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2012-2013 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is released under the terms of the GNU General Public Licence (GPL), version 3. You can
% redistribute it and/or modify it under the terms of the GPL (either version 3 or, at your option,
% any later version). For further details, see the file COPYING or <http://www.gnu.org/licenses/>.

    % Defaults
    init_h     = 1;
    div        = 1.2;
    min_steps  = 3;
    max_steps  = 100;
    tf         = 2;
    d2f        = NaN;
    err        = realmax;

    % Overrides
    if nargin > 2
        options = varargin{1};

        if isfield(options,'init_h')
            init_h = options.init_h;
        end

        if isfield(options,'div')
            div = options.div;
        end

        if isfield(options,'min_steps')
            min_steps = options.min_steps;
        end

        if isfield(options,'max_steps')
            max_steps = options.max_steps;
        end

        if isfield(options,'tf')
            tf = options.tf;
        end
    end

    % Initialize matrix of polynomial interpolation values
    P = NaN(max_steps);

    % Initialize finite difference step
    h = init_h;

    % Approximate 2nd derivative at initial step
    P(1,1) = (f(x+h)-2*f(x)+f(x-h))/h^2;

    % Loop through rows of P (i.e., steps of h)
    for i = 2:max_steps

        % New step size
        h = h/div;

        % Approximate 2nd derivative at this step
        P(i,1) = (f(x+h)-2*f(x)+f(x-h))/h^2;

        % Use square of div for extrapolation because errors increase
        % quadratically with h (here, of course, they decrease quadratically
        % because we're reducing h...)
        divsq = div^2;
        t = divsq;

        % Fill the current row using Richardson extrapolation
        for j = 2:i

            % Richardson
            P(i,j) = (t*P(i,j-1)-P(i-1,j-1))/(t-1);

            % Increment extrapolation factor
            t = t*divsq;

            % Error on this trial is defined as the maximum absolute difference
            % to the extrapolation parents
            currerr = max(abs(P(i,j)-P(i,j-1)),abs(P(i,j)-P(i-1,j-1)));

            if currerr < err
                err = currerr;
                d2f = P(i,j);
            end
        end

        % Stop if errors start increasing (to be expected for very small
        % values of h)
        if i > min_steps && abs(P(i,i)-P(i-1,i-1)) > tf*err
            return
        end
    end
end

author: translationalneuromodeling

repo: tapas


text:
function H = histcImWin( I, edges, wtMask, shape )
% Calculates local histograms at every point in an image I.
%
% H(i,j,...,k,:) will contain the histogram at location (i,j,...,k), as
% calculated by weighing values in I by placing wtMask at that
% location.  For example, if wtMask is ones(windowSize) then the
% histogram at every location will simply be a histogram of the pixels
% within that window.  See histc2 for more information about histgorams.
% See convnFast for information on shape flags.
%
% USAGE
%  H = histcImWin( I, edges, wtMask, [shape] )
%
% INPUTS
%  I           - image (possibly multidimensional) [see above]
%  edges       - quantization bounds, see histc2
%  wtMask      - numeric array of weights, or cell array of sep kernels
%  shape       - ['full'], 'valid', 'same', or 'smooth'
%
% OUTPUTS
%  H           - [size(I)xnBins] array of size(I) histograms
%
% EXAMPLE
%  load trees; L=conv2(X,filterDog2d(10,4,1,0),'valid'); figure(1); im(L);
%  f1=filterGauss(25,[],25);  f2=ones(1,15);
%  H1 = histcImWin(L, 15, {f1,f1'}, 'same');  figure(2); montage2(H1);
%  H2 = histcImWin(L, 15, {f2,f2'}, 'same');  figure(3); montage2(H2);
%
% See also ASSIGNTOBINS, HISTC2, CONVNFAST, HISTCIMLOC
%
% Piotr's Computer Vision Matlab Toolbox      Version 2.0
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

if( nargin<4 || isempty(shape) ); shape = 'full';  end;
if( ~iscell(wtMask) ); wtMask={wtMask}; end;

% split I into channels
I = assignToBins( I, edges );
nBins=length(edges)-1; if(nBins==0); nBins=edges; end;
nd = ndims(I); siz=size(I);  maxI = max(I(:));
if( nd==2 && siz(2)==1); nd=1; siz=siz(1); end;
QI = false( [siz maxI] );
inds = {':'}; inds = inds(:,ones(1,nd));
for i=1:nBins;  QI(inds{:},i)=(I==i); end;
H = double( QI );

% convolve with wtMask to get histograms, scale appropriately
for i=1:length(wtMask)
  wtMaski = wtMask{i};
  for d=1:ndims(wtMaski); wtMaski = flipdim(wtMaski,d); end;
  wtMaski = wtMaski / sum(wtMaski(:));
  H = convnFast( H, wtMaski, shape );
end;

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
classdef net_model_acc < mp.net_model_ac & mp.form_acc

%   MATPOWER
%   Copyright (c) 2019-2022, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

    properties
        vr = [];
        vi = [];
    end

    methods
        function obj = net_model_acc()
            obj@mp.net_model_ac();
            obj.element_classes = ...
                { @mp.nme_bus_acc, @mp.nme_gen_acc, @mp.nme_load_acc, ...
                    @mp.nme_branch_acc, @mp.nme_shunt_acc };

            %% Due to a bug related to inheritance in constructors in
            %% Octave 5.2 and earlier (https://savannah.gnu.org/bugs/?52614),
            %% INIT_SET_TYPES() cannot be called directly in the
            %% MP_IDX_MANAGER constructor, as desired.
            %%
            %% WORKAROUND:  INIT_SET_TYPES() is called explicitly as needed
            %%              (if obj.node is empty) in BUILD() and DISPLAY(),
            %%              after object construction, but before object use.
        end

        function obj = def_set_types(obj)
            def_set_types@mp.net_model_ac(obj);     % call parent first
            obj.set_types.vr = 'REAL VOLTAGE VARS (vr)';
            obj.set_types.vi = 'IMAG VOLTAGE VARS (vi)';
        end

        function va = initial_voltage_angle(obj, idx)
            vr = obj.params_var('vr');  %% inital value
            vi = obj.params_var('vi');  %% inital value
            if nargin < 2 || isempty(idx)
                va = angle(vr + 1j * vi);
            else
                va = angle(vr(idx) + 1j * vi(idx));
            end
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
% run sequence similarity

% whether to run the network to extract feature vectors of the objects
EXTRACT_FEATURES = true;

% whether to calculate the similarity matrix of the objects across frames
CALC_SIMILARITY = false;

CALC_DEEP_MATCH = false;

im_name = ['MOT-09/000001.jpg'];
img = imread(im_name);
[im_H,im_W,im_C] = size(img);
GT = csvread('MOT-09-gt/gt.txt');
W = load('net_weights.mat');
cls_ind = 16; % 1 background, 3 bike, 16 human
%frames = [110,115,120];
%ids = [7,22,44];
frames = [21];
ids = [19];

%{
GET FT VECTORS OF EACH BBOX:

  for each image:
    1. get rois, process it to be in original image pixel unit
    2. keep track of (object index, roi index) pairs
    3. do forward and backward pass of network
    4. extract feature vector (pool_5, contrib_img), save to dictionary
%}

% preprocess gt
gt = GT(:,1:6);
gt(:,5) = gt(:,5) + gt(:,3);
gt(:,6) = gt(:,6) + gt(:,4);
gt_ind = find_gt_idx_frame_obj(gt, frames, ids);
gt = gt(gt_ind,:);
for r=1:size(gt,1)
    if gt(r,3)>im_W
        gt(r,3) = im_W;
    end
    if gt(r,4)>im_H
        gt(r,4) = im_H;
    end
    if gt(r,5)>im_W
        gt(r,5) = im_W;
    end
    if gt(r,6)>im_H
        gt(r,6) = im_H;
    end
end

% resize img
for frame = frames
    im_name = ['MOT-09/' sprintf('%06d',frame) '.jpg'];
    img = imread(im_name);
    resized(frame).img = resize_img(img);
end

% get feature vector per frame
% FT_VEC is a matrix of struct
if EXTRACT_FEATURES

    for frame = frames
        gt_ind = find(gt(:,1) == frame);
        rois = gt(gt_ind, 2:6);
        rois(1,2:5) = [1200,300,1920,780];
        obj_ids = gt(gt_ind, 2);

        % setup
        im_name = ['MOT-09/' sprintf('%06d',frame) '.jpg'];
        img = imread(im_name);

        %forward & backward pass the network
        tic;
        res = net_forward_pass(img, W, rois, true);toc
        tic;
        contrib = net_reverse_contrib(res, W, cls_ind, 'img');toc
        contrib2 = net_reverse_contrib(res, W, 1, 'img');
        for roi_ind = 1:size(rois, 1)
            obj_id = obj_ids(roi_ind);
            FT_VEC(frame, obj_id).img = img;
            FT_VEC(frame, obj_id).roipool = squeeze(res.pool_5(:,:,:,roi_ind));
            FT_VEC(frame, obj_id).contrib_roipool = squeeze(contrib.roipool(roi_ind,:,:,:));
            FT_VEC(frame, obj_id).contrib_conv1_1 = squeeze(contrib.conv1_1(roi_ind,:,:,:));
            FT_VEC(frame, obj_id).contrib_img = squeeze(contrib.img(roi_ind,:,:,:));
        end
    end

end

% get cropped img of rois
for frame=frames
    for id=ids
        % get cropped matrix of rois
        im_name = ['MOT-09/' sprintf('%06d',frame) '.jpg'];
        img = imread(im_name);
        gt_ind = find_gt_idx_frame_obj(gt, [frame], [id]);
        roi = gt(gt_ind, 3:6);
        roi = ceil(roi);
        mat = img(roi(2):roi(4), roi(1):roi(3), :);
        im_name = ['output/' sprintf('%03d',id) '_' sprintf('%03d',frame) '.jpg'];
        imwrite(mat,im_name);
        cropped(frame,id).mat = mat;
    end
end

% generate similarity matrix:
frs = frames;
num_obj = length(frs)*length(ids);
SIM_roipool = zeros(num_obj,num_obj);
SIM_DM = SIM_roipool;
SIM_contrib_roipool = SIM_DM;

if CALC_SIMILARITY
    count = 0;
    for A_idx = 1:num_obj % A_idx and B_idx are index in the similarity matrix
        for B_idx = A_idx:num_obj
            count = count + 1;
            disp(count);

            A_id = ids( floor((A_idx-1)/length(ids)) + 1 ); % the actual frame number
            A_fr = frs( mod((A_idx-1), length(ids)) + 1 ); % the actual id number of the object
            B_id = ids( floor((B_idx-1)/length(ids)) + 1 );
            B_fr = frs( mod((B_idx-1), length(ids)) + 1 );

            A_mat = cropped(A_fr,A_id).mat;
            B_mat = cropped(B_fr,B_id).mat;

            A = FT_VEC(A_fr,A_id);
            B = FT_VEC(B_fr,B_id);
            SIM_roipool(A_idx,B_idx) = mat_sim(A, B,'roipool');
            SIM_contrib_roipool(A_idx,B_idx) = mat_sim(A, B,'contrib_roipool');

            A_chist = my_color_hist(A_mat);
            B_chist = my_color_hist(B_mat);
            SIM_color_hist(A_idx,B_idx) = mat_sim(A_chist,B_chist,'color_hist');

            A_chist = my_color_hist_weight(resized(A_fr).img,sum(A.contrib_conv1_1,3));
            B_chist = my_color_hist_weight(resized(B_fr).img,sum(B.contrib_conv1_1,3));
            SIM_color_hist_conv1_contrib(A_idx,B_idx) = mat_sim(A_chist,B_chist,'color_hist');
            if CALC_DEEP_MATCH
                 matches = my_deepmatch(A_mat, B_mat);
                 SIM_DM(A_idx,B_idx) = size(matches, 1)/sqrt(length(A_mat(:))*length(B_mat(:)));
            end

        end
    end

end

% post processing
SIM_roipool = matrix_max_percent(SIM_roipool);
SIM_DM = matrix_max_percent(SIM_DM);
SIM_color_hist = matrix_max_percent(SIM_color_hist);
SIM_color_hist_conv1_contrib = matrix_max_percent(SIM_color_hist_conv1_contrib);
SIM_contrib_roipool = matrix_max_percent(SIM_contrib_roipool);

%for roi_ind = 1:size(rois,1)
%    box = squeeze(rois(roi_ind,:));
%    rectangle('Position',[box(2),box(3),box(4)-box(2),box(5)-box(3)]);
%end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [ xp ] = DoAction( force, x )
%MountainCarDoAction: executes the action (a) into the mountain car
%environment
% a: is the force to be applied to the car
% x: is the vector containning the position and speed of the car
% xp: is the vector containing the new position and velocity of the car

position = x(1);
speed    = x(2);

% bounds for position
bpleft=-1.5;
bpright=0.5;

% bounds for speed
bsleft=-0.07;
bsright=0.07;

speedt1= speed + (0.001*force) + (-0.0025 * cos( 3.0*position) );
%speedt1= speedt1 * 0.999; % thermodynamic law, for a more real system with friction.

if(speedt1<bsleft)
    speedt1=bsleft;
end
if(speedt1>bsright)
    speedt1=bsright;
end

post1 = position + speedt1;

if(post1<=bpleft)
    post1=bpleft;
    speedt1=0.0;
end

if(post1>=bpright)
    post1=bpright;
    speedt1=0.0;
end
xp=[];
xp(1) = post1;
xp(2) = speedt1;

author: epfl-lasa


text:
%+========================================================================+
%|                                                                        |
%|            This script uses the GYPSILAB toolbox for Matlab            |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal, Marc Bakry (c) 2017-2019.                 |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%|             marc.bakry@polytechnique.edu                               |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : nrtVibroSlab2d.m                              |
%|    #    |   VERSION    : 0.55                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal & Marc Bakry                  |
%|  ( # )  |   CREATION   : 14.03.2019                                    |
%|  / 0 \  |   LAST MODIF :                                               |
%| ( === ) |   SYNOPSIS   :                                               |
%|  `---'  |                                                              |
%+========================================================================+

% Cleaning
clear all
close all
clc

% Gypsilab path
run('../../addpathGypsilab.m')

% Width of the slab
e = 0.5

% Frequency
f = 900:500:4000

% Incident direction (from bottom)
X0  = [0 1 0];

% Exterior domain (water)
rho0 = 1000;         % density      (kg.m3)
c0   = 1500;         % celerity     (m.s-1)
k0   = 2*pi/c0.*f;   % wave-number  (m-1)
lam0 = c0./f;        % wave-length  (m)

% Interior domain (different from water)
rhoS = 2*rho0;                                 % density                             (kg.m3)
cL   = 2*c0;                                   % celerity of longitudinal waves      (m.s-1)
kL   = 2*pi.*f./cL;                            % wave-number of longitudinal waves   (m-1)
lamL = real(cL)./f;                            % wavelength of longitudinal waves    (m)
cT   = 0;                                      % celerity of transverse waves        (m.s-1)
kT   = 2*pi.*f./cT;                            % wave-number of transverse waves     (m-1)
lamT = real(cT)./f;                            % wavelength of transverse waves      (m)

% Solution (pressure)
sol = zeros(2,length(f));

% Loop for each frequency
for i = 1:length(f)
    % Minimum wavelength
    tmp  = [lam0(i),lamL(i),lamT(i)];
    lmin = min(tmp(tmp>0));

    % Slab mesh
    L    = 20 * lmin;            % 20 wavelength to simulate infinite slab
    nx   = ceil(L/lmin * 6)+1;   % 6 node per wavelength for L
    ny   = ceil(e/lmin * 12)+1;  % 12 node per wavelength for e
    N    = nx * ny;              % Total number of nodes
    mesh = mshSquare(N,[L e])

    % Radiative mesh (fixed number of nodes)
    radiat = mshSquare(1e3,[L L]);

    % Boundary
    bound = swap(mesh.bnd)

    % Measurement points for trans and refl coeff (1 wavelenth from bound)
    Xmes = [0 -e/2-lmin 0 ; 0 e/2+lmin 0];

    % Cut-off function (50% full, 10% decrease)
    cutoff = vibsCutoff(1,L/5,L/10);

    % Green kernel function
    Gxy         = @(X,Y) femGreenKernel(X,Y,'[H0(kr)]',k0(i));
    gradyGxy{1} = @(X,Y) femGreenKernel(X,Y,'grady[H0(kr)]1',k0(i));
    gradyGxy{2} = @(X,Y) femGreenKernel(X,Y,'grady[H0(kr)]2',k0(i));
    gradyGxy{3} = @(X,Y) femGreenKernel(X,Y,'grady[H0(kr)]3',k0(i));

    % Plane wave function
    PW         = @(X) exp(1i*k0(i)*X*X0') .* cutoff(X);
    gradxPW{1} = @(X) 1i*k0(i)*X0(1) .* PW(X);
    gradxPW{2} = @(X) 1i*k0(i)*X0(2) .* PW(X);
    gradxPW{3} = @(X) 1i*k0(i)*X0(3) .* PW(X);

    % Coupling coeff for Brackage-Werner simulation
    beta = 1i*k0(i);

    % Graphical representation
    figure(1); clf;
    plot(radiat,real(PW(radiat.vtx)))
    hold on
    plot(mesh)
    plot(bound,'r')
    %     plotNrm(bound,'r')
    plot(msh(Xmes))
    axis equal
    colorbar
    title('Incident Wave')
    hold off

    % Quadrature and finite elements (volumn)
    omega = dom(mesh,3);
    U     = fem(mesh,'P1');

    % Quadrature and finite elements (boundary)
    sigma = dom(bound,3);
    u     = fem(bound,'P1');

    % Left-hand side
    tic
    [A,B,C,D] = vibsBlockOperator(omega,U,sigma,u,cL,cT,rhoS,c0,rho0,f(i));
    toc

    % Add dirichlet condition to x unknows (penalization)
    A(sub2ind(size(A),1:length(U),1:length(U))) = 1e15;

    % Right-hand side
    V    = cell(3,1);
    V{1} = - integral(sigma,ntimes(U,1),PW);
    V{2} = - integral(sigma,ntimes(U,2),PW);
    V{3} = integral(sigma,ntimes(u),gradxPW);

    % Resolution with Schur complement
    tic
    LHS    = [A B ; C D];
    RHS    = cell2mat(V);
    mu     = LHS \ RHS;
    mu     = mu(2*length(U)+1:end,:);
    lambda = beta*mu;
    toc

    % Measure of refexive and transmitted coeff
    tic
    Pmes = 1i/4 .* integral(Xmes,sigma,Gxy,u)*lambda - ...
        1i/4 .* integral(Xmes,sigma,gradyGxy,ntimes(u))*mu;
    Pmes(2) = Pmes(2) + PW(Xmes(2,:));
    toc

    % Save solution
    sol(:,i) = Pmes;

    % Finite element radiative operator --> \int_Sy G(x,y) psi(y) dy
    tic
    Srad = 1i/4 .* integral(radiat.vtx,sigma,Gxy,u);
    Srad = Srad - 1/(2*pi) .* regularize(radiat.vtx,sigma,'[log(r)]',u);

    % Finite element radiative operator --> \int_Sy dnyG(x,y) psi(y) dy
    Drad = 1i/4 .* integral(radiat.vtx,sigma,gradyGxy,ntimes(u));
    Drad = Drad - 1/(2*pi) .* regularize(radiat.vtx,sigma,'grady[log(r)]',ntimes(u));
    toc

    % Domain solution
    Psca = Srad*lambda - Drad*mu;
    Pinc = PW(radiat.vtx);
    Ptot = Psca + Pinc;

    % Graphical repesentation
    figure(10); clf;
    subplot(1,2,1)
    plot(bound)
    hold on
    plot(radiat,abs(Psca))
    plot(msh(Xmes),abs(Pmes))
    title('Scattered')
    axis equal
    colorbar
    hold off

    subplot(1,2,2)
    plot(bound)
    hold on
    plot(radiat,abs(Ptot))
    plot(msh(Xmes),abs(Pmes))
    title('Total')
    axis equal
    colorbar
    hold off
end

% Analytical solution
tic
c     = ones(length(f),1) * [c0 cL c0];
rho   = [rho0 rhoS rho0];
[R,T] = slabVibro(f,rho,c,e);
toc

% Comparison (db)
ref = 20*log10(abs([R ; T]));
sol = 20*log10(abs(sol));

norm(ref-sol)/norm(ref)

% Graphical representation
figure(100)
subplot(1,2,1)
plot(f,ref(1,:),'r',f,sol(1,:),'b+')
grid on
title('Reflexion coeffiscient')
legend({'Analytical','Numerical'})
xlabel('Frequency (Hz)')
ylabel('Amplitude (dB)')

subplot(1,2,2)
plot(f,ref(2,:),'r',f,sol(2,:),'b+')
grid on
title('Transmission coeffiscient')
legend({'Analytical','Numerical'})
xlabel('Frequency (Hz)')
ylabel('Amplitude (dB)')

disp('~~> Michto gypsilab !')

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
%SerialLink.plot3d Graphical display and animation of solid model robot
%
% R.plot3d(Q, options) displays and animates a solid model of the robot.
% The robot is displayed at the joint angle Q (1xN), or
% if a matrix (MxN) it is animated as the robot moves along the M-point trajectory.
%
% Options::
%
% 'color',C         A cell array of color names, one per link.  These are
%                   mapped to RGB using colorname().  If not given, colors
%                   come from the axis ColorOrder property.
% 'alpha',A         Set alpha for all links, 0 is transparant, 1 is opaque
%                   (default 1)
% 'path',P          Overide path to folder containing STL model files
% 'workspace', W    Size of robot 3D workspace, W = [xmn, xmx ymn ymx zmn zmx]
% 'floorlevel',L    Z-coordinate of floor (default -1)
%-
% 'delay',D         Delay betwen frames for animation (s)
% 'fps',fps         Number of frames per second for display, inverse of 'delay' option
% '[no]loop'        Loop over the trajectory forever
% '[no]raise'       Autoraise the figure
% 'movie',M         Save frames as files in the folder M
% 'trail',L         Draw a line recording the tip path, with line style L.
%                   L can be a cell array, eg. {'r', 'LineWidth', 2}
%-
% 'scale',S         Annotation scale factor
% 'ortho'           Orthographic view (default)
% 'perspective'     Perspective view
% 'view',V          Specify view V='x', 'y', 'top' or [az el] for side elevations,
%                   plan view, or general view by azimuth and elevation
%                   angle.
%-
% '[no]wrist'       Enable display of wrist coordinate frame
% 'xyz'             Wrist axis label is XYZ
% 'noa'             Wrist axis label is NOA
% '[no]arrow'       Display wrist frame with 3D arrows
%-
% '[no]tiles'       Enable tiled floor (default true)
% 'tilesize',S      Side length of square tiles on the floor (default 0.2)
% 'tile1color',C   Color of even tiles [r g b] (default [0.5 1 0.5]  light green)
% 'tile2color',C   Color of odd tiles [r g b] (default [1 1 1] white)
%-
% '[no]jaxes'       Enable display of joint axes (default true)
% '[no]joints'      Enable display of joints
%-
% '[no]base'        Enable display of base shape
%
% Notes::
% - Solid models of the robot links are required as STL files (ascii or
%   binary) with extension .stl.
% - The solid models live in RVCTOOLS/robot/data/meshes.
% - Each STL model is called 'linkN'.stl where N is the link number 0 to N
% - The specific folder to use comes from the SerialLink.model3d property
% - The path of the folder containing the STL files can be overridden using
%   the 'path' option
% - The height of the floor is set in decreasing priority order by:
%   - 'workspace' option, the fifth element of the passed vector
%   - 'floorlevel' option
%   - the lowest z-coordinate in the link1.stl object
%
% Making an animation::
%
% The 'movie' options saves the animation as a movie file or separate frames in a folder
% - 'movie','file.mp4' saves as an MP4 movie called file.mp4
% - 'movie','folder' saves as files NNNN.png into the specified folder
%   - The specified folder will be created
%   - NNNN are consecutive numbers: 0000, 0001, 0002 etc.
%   - To convert frames to a movie use a command like:
%          ffmpeg -r 10 -i %04d.png out.avi
%
% Authors::
% - Peter Corke, based on existing code for plot().
% - Bryan Moutrie, demo code on the Google Group for connecting ARTE and
%   RTB.
%
% Acknowledgments::
% - STL files are from ARTE: A ROBOTICS TOOLBOX FOR EDUCATION by Arturo Gil
%   (https://arvc.umh.es/arte) are included, with permission.
% - The various authors of STL reading code on file exchange, see stlRead.m
%
% See also SerialLink.plot, plotbotopt3d, SerialLink.animate, SerialLink.teach, stlRead.

% Copyright (C) 1993-2015, 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

function plot3d(robot, q, varargin)

    assert( ~robot.mdh, 'RTB:plot3d:badmodel', '3D models are defined for standard, not modified, DH parameters');

    opt = plot_options(robot, varargin);

    robot.movie = Animate(opt.movie);

    %-- load the shape if need be

    nshapes = robot.n+1;

    if isempty(robot.faces)
        % no 3d model defined, let's try to load one

        if isempty(opt.path)
            % first find the path to the models
            pth = mfilename('fullpath');
            pth = fileparts(pth);
            % peel off the last folder
            s = regexp(pth, filesep, 'split');
            pth = join(s(1:end-1), filesep);

            % find the path to this specific model
            pth = fullfile(pth{1}, 'data/meshes', robot.model3d);
            assert(exist(pth, 'dir') > 0, 'RTB:plot3d:nomodel', 'no 3D model found, install the RTB contrib zip file');
        else
            pth = opt.path;
        end

        % now load the STL files
        robot.points = cell(1, robot.n+1);
        robot.faces = cell(1, robot.n+1);
        fprintf('Loading STL models from ARTE Robotics Toolbox for Education  by Arturo Gil (http://arvc.umh.es/arte)');
        for i=1:nshapes
            [P,F] = stlRead( fullfile(pth, sprintf('link%d.stl', i-1)) );
            robot.points{i} = P;
            robot.faces{i} = F;
            fprintf('.');

        end
        fprintf('\n');
    end

    % if a base is specified set the floor height to this
    if isempty(opt.workspace)
        % workspace not provided, fall through the options for setting floor level
        if ~isempty(opt.floorlevel)
            opt.ws(5) = opt.floorlevel;
        elseif opt.base
            mn = min( robot.points{1} );
            opt.ws(5) = mn(3);
        end
    end
    opt.floorlevel = opt.ws(5);

% TODO
% should test if the plot exists, pinch the logic from plot()

    %-- set up to plot
    % create an axis
    ish = ishold();
    if ~ishold
        % if hold is off, set the axis dimensions
        axis(opt.ws);
        set(gca, 'ZLimMode', 'manual');
        axis(opt.ws);
        hold on
    end

    if opt.raise
        % note this is a very time consuming operation
        figure(gcf);
    end

    if strcmp(opt.projection, 'perspective')
        set(gca, 'Projection', 'perspective');
    end

    grid on
    xlabel('X'); ylabel('Y'); zlabel('Z');

    %--- create floor if required
    if opt.tiles
        create_tiled_floor(opt);
    end

    %--- configure view and lighting
    if isstr(opt.view)
        switch opt.view
            case 'top'
                view(0, 90);
            case 'x'
                view(0, 0);
            case 'y'
                view(90, 0)
            otherwise
                error('rtb:plot3d:badarg', 'view must be: x, y, top')
        end
    elseif isnumeric(opt.view) && length(opt.view) == 2
        view(opt.view)
    else
        campos([2 2 1]);
    end

    daspect([1 1 1]);
    light('Position', [0 0 opt.reach*2]);
    light('Position', [1 0.5 1]);

    %-- figure the colors for each shape
    if isempty(opt.color)
        % if not given, use the axis color order
        C = get(gca,'ColorOrder');
    else
        C = [];
        for c=opt.color
            C = [C; colorname(c{1})];
        end
    end

    %--- create the robot
    %  one patch per shape, use hgtransform to animate them later
    group = hggroup('Tag', robot.name);
    ncolors = numrows(C);

    h = [];
    for link=0:robot.n
        if link == 0
            if ~opt.base
                continue;
            end

            patch('Faces', robot.faces{link+1}, 'Vertices', robot.points{link+1}, ...
                'FaceColor', C(mod(link,ncolors)+1,:), 'EdgeAlpha', 0, 'FaceAlpha', opt.alpha);
        else
            h.link(link) = hgtransform('Tag', sprintf('link%d', link), 'Parent', group);

            patch('Faces', robot.faces{link+1}, 'Vertices', robot.points{link+1}, ...
                'FaceColor', C(mod(link,ncolors)+1,:), 'EdgeAlpha', 0, 'FaceAlpha', opt.alpha, ...
            'Parent', h.link(link));
        end
    end

        % display the wrist coordinate frame
    if opt.wrist
        if opt.arrow
            h.wrist = trplot(eye(4,4), 'labels', upper(opt.wristlabel), ...
                'arrow', 'rgb', 'length', 0.4);
        else
            h.wrist = trplot(eye(4,4), 'labels', upper(opt.wristlabel), ...
                'rgb', 'length', 0.4);
        end
    else
        h.wrist = [];
    end

    if ~isempty(opt.trail)
        h.trail = plot(0, 0, opt.trail{:});
        robot.trail = [];
    end

    % enable mouse-based 3D rotation
    rotate3d on

    h.robot = robot;
    h.link = [0 h.link];
    set(group, 'UserData', h);

    robot.trail = []; % clear the previous trail

    robot.animate(q);

    if opt.movie
        robot.movie.close();
    end

    if ~ish
        hold off
    end
end

function opt = plot_options(robot, optin)
    opt.color = [];
    opt.path = [];  % override path
    opt.alpha = 1;

        % timing/looping
    opt.delay = 0.1;
    opt.fps = [];
    opt.loop = false;

    opt.raise = false;

    % general appearance
    opt.scale = 1;
    opt.trail = [];

    opt.workspace = [];
    opt.floorlevel = [];

    opt.name = true;
    opt.projection = {'ortho', 'perspective'};
    opt.view = [];

    % tiled floor
    opt.tiles = true;
    opt.tile1color = [0.5 1 0.5];  % light green
    opt.tile2color = [1 1 1];  % white
    opt.tilesize = 0.2;

    % the base or pedestal
    opt.base = true;

    % wrist
    opt.wrist = true;
    opt.wristlabel = {'xyz', 'noa'};
    opt.arrow = true;

    % joint rotation axes
    opt.jaxes = false;

    % misc
    opt.movie = [];

    % build a list of options from all sources
    %   1. the M-file plotbotopt if it exists
    %   2. robot.plotopt
    %   3. command line arguments
    if exist('plotbotopt3d', 'file') == 2
        options = [plotbotopt3d robot.plotopt3d optin];
    else
        options = [robot.plotopt3d optin];
    end

    % parse the options
    [opt,args] = tb_optparse(opt, options);
    if ~isempty(args)
        error(['unknown option: ' args{1}]);
    end

    if ~isempty(opt.fps)
        opt.delay = 1/opt.fps;
    end
    % figure the size of the figure

    if isempty(opt.workspace)
        %
        % simple heuristic to figure the maximum reach of the robot
        %
        L = robot.links;
        if any(L.isprismatic)
            error('Prismatic joint(s) present: requires the ''workspace'' option');
        end
        reach = 0;
        for i=1:robot.n
            reach = reach + abs(L(i).a) + abs(L(i).d);
        end

%         if opt.wrist
%             reach = reach + 1;
%         end

        % if we have a floor, quantize the reach to a tile size
        if opt.tiles
            reach = opt.tilesize * ceil(reach/opt.tilesize);
        end

        % now create a 3D volume based on this reach
        opt.ws = [-reach reach -reach reach -reach reach];

        % if a floorlevel has been given, ammend the 3D volume
        if ~isempty(opt.floorlevel)
            opt.ws(5) = opt.floorlevel;
        end
    else
        % workspace is provided
        reach = min(abs(opt.workspace));
        if opt.tiles
            % set xy limits to be integer multiple of tilesize
            opt.ws(1:4) = opt.tilesize * round(opt.workspace(1:4)/opt.tilesize);
            opt.ws(5:6) = opt.workspace(5:6);
        else
            opt.ws=opt.workspace;
        end
    end

    opt.reach = reach;

    % update the fundamental scale factor (given by the user as a multiplier) by a length derived from
    % the overall workspace dimension
    %  we need that a lot when creating the robot model
    opt.scale = opt.scale * reach/40;

    % deal with a few options that need to be stashed in the SerialLink object

    robot.delay = opt.delay;
    robot.loop = opt.loop;
end

% draw a tiled floor in the current axes
function create_tiled_floor(opt)

    xmin = opt.ws(1);
    xmax = opt.ws(2);
    ymin = opt.ws(3);
    ymax = opt.ws(4);

    % create a colored tiled floor
    xt = xmin:opt.tilesize:xmax;
    yt = ymin:opt.tilesize:ymax;
    Z = opt.floorlevel*ones( numel(yt), numel(xt));
    C = zeros(size(Z));
    [r,c] = ind2sub(size(C), 1:numel(C));
    C = bitand(r+c,1);
    C = reshape(C, size(Z));
    C = cat(3, opt.tile1color(1)*C+opt.tile2color(1)*(1-C), ...
        opt.tile1color(2)*C+opt.tile2color(2)*(1-C), ...
        opt.tile1color(3)*C+opt.tile2color(3)*(1-C));
    [X,Y] = meshgrid(xt, yt);
    surface(X, Y, Z, C, ...
        'FaceColor','texturemap',...
        'EdgeColor','none',...
        'CDataMapping','direct');
end

author: petercorke

repo: robotics-toolbox-matlab


text:
function H = higuchi(sequence,isplot)
%
% 'higuchi' estimate the hurst parameter of a given sequence with
%     higuchi's method.
%
% Inputs:
%     sequence: the input sequence for estimate
%     isplot: whether display the plot. without a plot if isplot equal to 0
% Outputs:
%     H: the estimated hurst coeffeient of the input sequence

%  Author: Chu Chen
%  Version 1.0,  03/10/2008
%  chen-chu@163.com
%

if nargin == 1
    isplot = 0;
end

sequence = cumsum(sequence);
N = length(sequence);
mlarge = floor(N/5);
M = [floor(logspace(0,log10(mlarge),50))];
M = unique(M(M>1));
n = length(M);
cut_min = ceil(n/10);
cut_max = floor(6*n/10);

curve_length = zeros(1,n);
for h = 1:n
    m = M(h);
    k = floor((N-m)/m);
    temp_length = zeros(m,k);

    for i = 1:m
        for j = 1:k
            temp_length(i,j) = abs(sequence(i+j*m)-sequence(i+(j-1)*m));
        end
    end

    curve_length(h) = sum(mean(temp_length,2)) * ((N-1)/m^3);
end

x = log(M);
y = log(curve_length);
X = x(cut_min:cut_max);
Y = y(cut_min:cut_max);
p1 = polyfit(X,Y,1);
Yfit = polyval(p1,X);
yfit = polyval(p1,x);
H = 2 + (Yfit(end)-Yfit(1))/(X(end)-X(1));

if isplot ~= 0
    figure,hold on;
    plot(x,y,'b*');
    plot(X,Yfit,'r-','LineWidth',2);
    plot(x(1:cut_min),yfit(1:cut_min),'r:','LineWidth',2);
    plot(x(cut_max:end),yfit(cut_max:end),'r:','LineWidth',2);
    xlabel('Log10(Aggregate Level)'),ylabel('Log10(Curve Legnth)'),title('Higuchi Method');
end

author: Sable

repo: mcbench-benchmarks


text:
function nav=readerp(nav,fname)

global glc gls
flag=0; n=0; erp=gls.erp;

idx=strfind(fname,glc.sep); fname0=fname(idx(end)+1:end);
fprintf('Info:reading erp file %s...',fname0);
fid=fopen(fname);

% read data
while ~feof(fid)
    line=fgets(fid,266);

    if ~isempty(strfind(line,'MJD')),flag=1;end
    if flag==0,continue;end
    if flag==1
        line=fgets(fid); %#ok<NASGU>
        flag=2; continue;
    end
    if flag==2
        value=str2double(strsplit(line));
        erp.data(n+1).mjd    =value(1);
        erp.data(n+1).xp     =value(2)*1e-6*glc.AS2R;
        erp.data(n+1).yp     =value(3)*1e-6*glc.AS2R;
        erp.data(n+1).ut1_utc=value(4)*1e-7;
        erp.data(n+1).lod    =value(5)*1e-7;
        erp.data(n+1).xpr    =value(13)*1e-6*glc.AS2R;
        erp.data(n+1).ypr    =value(14)*1e-6*glc.AS2R;
        n=n+1;
    end
end

fclose(fid);
fprintf('over\n');

erp.n=n;
erp.nmax=n;
nav.erp=erp;

return

repo: GINav


text:
function [ph_uw,msd]=uw_3d(ph,xy,day,ifgday_ix,bperp,options)
%UW_3D unwrap phase time series (single or multiple master)
%   PH_UW = UW_3D(PH,XY,DAY,IFGDAY_IX,BPERP,OPTIONS)
%
%   PH  = N x M matrix of wrapped phase values (real phase or complex phasor)
%        where N is number of pixels and M is number of interferograms
%   XY  = N x 2 matrix of coordinates in metres
%        (optional extra column, in which case first column is ignored)
%   DAY = vector of image acquisition dates in days, relative to master
%   IFGDAY_IX = M x 2 matrix giving index to master and slave date in DAY
%        for each interferogram
%   BPERP  = M x 1 vector giving perpendicular baselines
%   OPTIONS structure optionally containing following fields:
%      la_flag    = look angle error estimation flag (def='y')
%      scf_flag   = spatial cost function estimation flag (def='y')
%      master_day = serial date number of master (used for info msg only, def=0)
%      grid_size  = size of grid in m to resample data to (def=5)
%      prefilt_win = size of prefilter window in resampled grid cells (def=16)
%      time_win   = size of time filter window in days (def=365)
%      unwrap_method (def='3D' for single master, '3D_FULL' otherwise)
%      goldfilt_flag = Goldstein filtering, 'y' or 'n' (def='n')
%      gold_alpha    = alpha value for Goldstein filter (def=0.8)
%      lowfilt_flag  = Low pass filtering, 'y' or 'n' (def='n')
%      n_trial_wraps = no. phase cycles poss between neighbouring points due to la error (def=6)
%      temp          = optional M x 1 vector of temperature difference for each ifg (def=[])
%      n_temp_wraps  = no. phase cycles poss between neighbouring points due to temp diff (def=2)
%      variance      = N x 1 matrix of variance values
%   PH_UW = unwrapped phase
%
%   Andy Hooper, Jun 2007

% ============================================================
% 08/2009 AH: Goldstein alpha value added to options
% 01/2012 AH: Changes for easier running stand-alone
% 01/2012 AH: Add bperp for new look angle error estimation
% 02/2012 AH: Add new method 3D_NEW
% 11/2012 AH: Several new options
% 02/2014 AH: Add predefined ph_uw option
% ============================================================
tic;
if nargin<4
    help uw_3d
    error('not enough arguments')
end

if nargin<5
    bperp=[];
end

if nargin<6
    options=struct;
end

%if isempty(ifgday_ix) | length(unique(ifgday_ix(:,1)))==1
if isempty(ifgday_ix)
    single_master_flag=1;
else
    single_master_flag=0;
end

valid_options={'la_flag','scf_flag','master_day','grid_size','prefilt_win','time_win','unwrap_method','goldfilt_flag','lowfilt_flag','gold_alpha','n_trial_wraps','temp','n_temp_wraps','max_bperp_for_temp_est','variance','ph_uw_predef'};
invalid_options=setdiff(fieldnames(options),valid_options);
if length(invalid_options)>0
    error(['"',invalid_options{1}, '" is an invalid option'])
end

if ~isfield(options,'master_day')
    options.master_day=0;
end

if ~isfield(options,'grid_size')
    options.grid_size=5;
end

if ~isfield(options,'prefilt_win')
    options.prefilt_win=16;
end

if ~isfield(options,'time_win')
    options.time_win=365;
end

if ~isfield(options,'unwrap_method')
    if single_master_flag==1
        options.unwrap_method='3D';
    else
        options.unwrap_method='3D_FULL';
    end
end

if ~isfield(options,'goldfilt_flag')
    options.goldfilt_flag='n';
end

if ~isfield(options,'lowfilt_flag')
    options.lowfilt_flag='n';
end

if ~isfield(options,'gold_alpha')
    options.gold_alpha=0.8;
end

if ~isfield(options,'n_trial_wraps')
    options.n_trial_wraps=6;
end

if ~isfield(options,'la_flag')
    options.la_flag='y';
end

if ~isfield(options,'scf_flag')
    options.scf_flag='y';
end

if ~isfield(options,'temp')
    options.temp=[];
else
    if ~isempty(options.temp) & length(options.temp)~=size(ph,2)
        error('options.temp must be M x 1 vector where M is no. of ifgs')
    end
end

if ~isfield(options,'n_temp_wraps')
    options.n_temp_wraps=2;
end

if ~isfield(options,'max_bperp_for_temp_est')
    options.max_bperp_for_temp_est=100;
end

if ~isfield(options,'variance')
    options.variance=[];
end

if ~isfield(options,'ph_uw_predef')
    options.ph_uw_predef=[];
end

if size(xy,2)==2
   xy(:,2:3)=xy(:,1:2);
end

if size(day,1)==1
    day=day';
end

if strcmpi(options.unwrap_method,'3D') | strcmpi(options.unwrap_method,'3D_NEW')
    if length(unique(ifgday_ix(:,1)))==1
        options.unwrap_method='3D_FULL';
    else
        options.lowfilt_flag='y';
    end
end

uw_grid_wrapped(ph,xy,options.grid_size,options.prefilt_win,options.goldfilt_flag,options.lowfilt_flag,options.gold_alpha,options.ph_uw_predef);
clear ph
uw_interp;
%if single_master_flag==1
%    uw_unwrap_space_time(day,options.unwrap_method,options.time_win,options.master_day,options.la_flag,bperp,options.n_trial_wraps,options.prefilt_win,options.scf_flag,options.temp,options.n_temp_wraps);
%else
    uw_sb_unwrap_space_time(day,ifgday_ix,options.unwrap_method,options.time_win,options.la_flag,bperp,options.n_trial_wraps,options.prefilt_win,options.scf_flag,options.temp,options.n_temp_wraps,options.max_bperp_for_temp_est);
%end
uw_stat_costs(options.unwrap_method,options.variance);
[ph_uw,msd]=uw_unwrap_from_grid(xy,options.grid_size);

author: dbekaert

repo: StaMPS

save_path: github-repos/MATLAB/dbekaert-StaMPS


text:
function G = GCVstopfun(alpha, u, s, beta, m, n)
%
%  G = GCVstopfun(alpha, u, s, beta, n, insolv)
%  This function evaluates the GCV function G(i, alpha), that will be used
%     to determine a stopping iteration.
%
% Input:
%   alpha - regularization parameter at the kth iteration of HyBR
%       u - P_k^T e_1 where P_k contains the left singular vectors of B_k
%       s - singular values of bidiagonal matrix B_k
%    beta - norm of rhs b
%     m,n - size of the ORIGINAL problem (matrix A)

% Silvia Gazzola, University of Bath
% Per Christian Hansen, Technical University of Denmark
% James G. Nagy, Emory University
% April, 2018.

% This file is part of the IR Tools package and is distributed under the
% 3-Clause BSD License. A separate license file should be provided as part
% of the package.

k = length(s);
beta2 = beta^2;

s2 = abs(s) .^ 2;
alpha2 = alpha^2;

t1 = 1 ./ (s2 + alpha2);
t2 = abs(alpha2*u(1:k) .* t1) .^2;
t3 = s2 .* t1;

num = beta2*(sum(t2) + abs(u(k+1))^2)/n;
den = ( (m - sum(t3))/n )^2;
G = num / den;

repo: IRtools


text:
function s=openRawGeneralVolumeWindow
%
% Opens a new volume window and initializes the corresponding data
% structure.
%
% VOLUME is a cell array of volume structures.
% s is the index of the new one.
%
% This function should ONLY be called by either:
% openVolumeWindow or openGrayWindow.
%
% Modifications:
%
% djh and baw, 7/98
% djh, 2/99 modified to enable the window to be opened under
% either volume or gray mode
%
% djh, 4/99
% - Eliminate overlayClip sliders
% - Added mapWin sliders to show overlay only for pixels with parameter
%   map values that are in the appropriate range.
% bw, 12/29/00
% - scan slider instead of buttons
% djh 2/13/2001
% - open multiple volume windows simultaneously

% Make sure the global variables exist

mrGlobals
disp('Initializing Volume view')

% s is the index of the new volume structure.
s = getNewViewIndex(VOLUME);

% Set name, viewType, & subdir
VOLUME{s}.name=['VOLUME{',num2str(s),'}'];
VOLUME{s}.viewType='generalGray';
VOLUME{s}.subdir='generalGray';
if(isfield(mrSESSION,'sessionCode'))
  VOLUME{s}.sessionCode=mrSESSION.sessionCode;
else
  VOLUME{s}.sessionCode='';
end

% Refresh function, gets called by refreshScreen
VOLUME{s}.refreshFn = 'refreshView';

%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize data slots %
%%%%%%%%%%%%%%%%%%%%%%%%%

% Should be VOLUME{s} = viewInit('VOLUME');

%Initialize slot for anat
VOLUME{s}.anat = [];

% Initialize slots for co, amp, and ph
VOLUME{s}.co = [];
VOLUME{s}.amp = [];
VOLUME{s}.ph = [];
VOLUME{s}.map = [];
VOLUME{s}.mapName = '';
VOLUME{s}.mapUnits = '';

% Initialize slots for tSeries
VOLUME{s}.tSeries = [];
VOLUME{s}.tSeriesScan = NaN;
VOLUME{s}.tSeriesSlice = NaN;

% Initialize ROIs
VOLUME{s}.ROIs = [];
VOLUME{s}.selectedROI = 0;

% Initialize curDataType
VOLUME{s}.curDataType = 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Make displayModes and color maps %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Initialize displayModes
VOLUME{s}=resetDisplayModes(VOLUME{s});
VOLUME{s}.ui.displayMode='anat';

%%%%%%%%%%%%%%%%%%%
% Open the window %
%%%%%%%%%%%%%%%%%%%
% Figure number for volume view window
figName = sprintf('Volume %i: %s  [%s]',s,...
                  mrSESSION.sessionCode,...
                  mrSESSION.description);

% open the figure
VOLUME{s}.ui.figNum=figure('MenuBar','none','Units','Normalized',...
                    'Color',[.9 .9 .9],'Position',[0.25 0.25 0.5 0.5],...
                    'NumberTitle','off','Name',figName);

% Handle for volume view window
VOLUME{s}.ui.windowHandle = gcf;

% Handle for main axis of volume view
VOLUME{s}.ui.mainAxisHandle = gca;

% Adjust position and size of main axes so there's room for
% buttons, colorbar, and sliders
set(VOLUME{s}.ui.mainAxisHandle,'position',[0.1 0.1 0.725 0.7]);

% Set minColormap property so there's potentially room for 128
% colors
set(VOLUME{s}.ui.windowHandle,'minColormap',128)

% Sharing of colors seems like it might be OK, but I'm turning it
% off just to be sure (djh, 1/26/98).
set(gcf,'sharecolors','off');

% Set closeRequestFcn so we can clean up when the window is closed
set(gcf,'CloseRequestFcn','closeVolumeWindow');

% Set selectedVOLUME when click in this window
set(gcf,'WindowButtonDownFcn',['selectedVOLUME =',num2str(s),';']);

%%%%%%%%%%%%%
% Add Menus %
%%%%%%%%%%%%%
disp('Attaching menus')
VOLUME{s}=filesMenu(VOLUME{s});
VOLUME{s}=editMenu(VOLUME{s});
VOLUME{s}=windowMenu(VOLUME{s});
VOLUME{s}=analysisMenu(VOLUME{s});
VOLUME{s}=viewMenu(VOLUME{s});
VOLUME{s}=roiMenu(VOLUME{s});
VOLUME{s}=plotMenu(VOLUME{s});
VOLUME{s}=colorMenu(VOLUME{s});
VOLUME{s}=xformVolumeMenu(VOLUME{s});
VOLUME{s}=segmentationMenu(VOLUME{s});

%%%%%%%%%%%%%%%%%%%%%%%%%
% Add Annotation String %
%%%%%%%%%%%%%%%%%%%%%%%%%
VOLUME{s} = makeAnnotationString(VOLUME{s});

%%%%%%%%%%%%%%%%%
% Add Color Bar %
%%%%%%%%%%%%%%%%%
% Make color bar and initialize it to 'off'
VOLUME{s}.ui.colorbarHandle=makeColorBar(VOLUME{s});
setColorBar(VOLUME{s},'off');
VOLUME{s}.ui.cbarRange = [];
% move to nicer location
set(VOLUME{s}.ui.colorbarHandle,'Position',[0.2 0.85 0.6 0.03]);

%%%%%%%%%%%%%%%
% Add Buttons %
%%%%%%%%%%%%%%%
disp('Attaching buttons')

% Buttons and editable text fields for choosing slice number and
% slice orientation
VOLUME{s} = makeVolSliceUI(VOLUME{s});
setCurSliceOri(VOLUME{s},1);

% zoom buttons
VOLUME{s} = makeZoomButtons(VOLUME{s});
VOLUME{s} = makeDomainButton(VOLUME{s});

%%%%%%%%%%%%%%%%%%%
% Add popup menus %
%%%%%%%%%%%%%%%%%%%
disp('Attaching popup menus')
VOLUME{s} = makeROIPopup(VOLUME{s});
VOLUME{s} = makeDataTypePopup(VOLUME{s});

%%%%%%%%%%%%%%%
% Add sliders %
%%%%%%%%%%%%%%%
disp('Attaching sliders')

% Scan number slider
w = 0.12; h = 0.04; l = 0; b = 0.95;
VOLUME{s} = makeSlider(VOLUME{s},'scan',[],[l b w h]);
VOLUME{s} = initScanSlider(VOLUME{s},1);
VOLUME{s} = selectDataType(VOLUME{s},VOLUME{s}.curDataType);

% correlation threshold:
VOLUME{s} = makeSlider(VOLUME{s},'cothresh',[0,1],[.85,.85,.15,.05]);
setCothresh(VOLUME{s},0);

% phase window:
VOLUME{s} = makeSlider(VOLUME{s},'phWinMin',[0,2*pi],[.85,.75,.15,.05]);
VOLUME{s} = makeSlider(VOLUME{s},'phWinMax',[0,2*pi],[.85,.65,.15,.05]);
setPhWindow(VOLUME{s},[0 2*pi]);

% parameter map window:
VOLUME{s} = makeSlider(VOLUME{s},'mapWinMin',[0,1],[.85,.55,.15,.05]);
VOLUME{s} = makeSlider(VOLUME{s},'mapWinMax',[0,1],[.85,.45,.15,.05]);
setMapWindow(VOLUME{s},[0 1]);

% anatClip: determines clipping of the anatomy base-image
%           values to fill the range of available grayscales.
VOLUME{s} = makeSlider(VOLUME{s},'anatMin',[0,1],[.85,.2,.15,.05]);
VOLUME{s} = makeSlider(VOLUME{s},'anatMax',[0,1],[.85,.1,.15,.05]);
setAnatClip(VOLUME{s},[0 .5]);

% Time or freq vals
% This is a complicated slider: It chooses values from the cell array in
% dataTYPES. These values indicate the current data to display in the
% window.
% For example, dataTYPES{1}.generalAnalysisParams(1).temporalAnalysis=[1:256];
%              dataTYPES{1}.generalAnalysisParams(1).frequencyAnalysis={[1 3 5],[2 4
%              6],[1:10],[60]}
% Depending on whether FFTCheckBox is set, this slider can assume values
% from 1 to 256 or from 1 to 4. As the slider is moved, the correct data is
% computed on the fly from the timeseries, the inverse and view.mesh.vertexGrayMap.
% It might be smart to store both the timeseries and the FFT separately in
% so that we don't have to compute the FFT each time. On the other hand,
% for reasonable size datasets (128*512 say..) the FFT takes only
% 1/10 sec.
% Note that unlike the regular gray view, this view >must< have the mesh loaded.
% You can do this (without displaying it) using meshLoad
w = 0.12; h = 0.04; l = 0; b = 0.65;
VOLUME{s} = makeSlider(VOLUME{s},'datavals',[0,1],[l b w h]);
VOLUME{s} = initDataValsSlider(VOLUME{s},1);
%VOLUME{s} = selectDataType(VOLUME{s},VOLUME{s}.curDataType);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize display parameters %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize image field
VOLUME{s}.ui.image = [];

% Show all ROIs
VOLUME{s}.ui.showROIs = 2;

%%%%%%%%%%%%%%%%%%
% Load Anatomies %
%%%%%%%%%%%%%%%%%%
VOLUME{s} = loadAnat(VOLUME{s});

%%%%%%%%%%%%%%%%%%
% Refresh screen %
%%%%%%%%%%%%%%%%%%
VOLUME{s}=refreshScreen(VOLUME{s});
disp('Done initializing Volume view')

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%compute orientation of tailcentralvector

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

larvae=trx.exp2flies{n};
numlarvae=numel(larvae);
tailcentralang=cell(1,numlarvae);

for i=1:numlarvae
    larva=larvae(i);
    tailcentralang{1,i}=bsxfun(@atan2,trx(larva).ycentral_mm-trx(larva).ytail_mm,trx(larva).xcentral_mm-trx(larva).xtail_mm);
end
units=parseunits('rad');
data=tailcentralang;

author: kristinbranson

repo: JAABA


text:
% clear all;
% close all;
% %%%pavia University
% input_files.data_file='./data/PaviaU.mat';
% input_files.gt_file='./data/PaviaU_gt.mat';
% output_file_name='./results/Pavia_U_pca.mat';
% pca_dim=10;
% prepare_data_gt(input_files, output_file_name, pca_dim);

addpath('./pca_ica/');

clear all;
close all;
%%%pavia University
input_files.data_file='./data/Indian_pines_corrected.mat';
input_files.gt_file='./data/indian_pines_gt.mat';
output_file_name='./results/Indian_pines_pca.mat';
pca_dim=30;
prepare_data_gt(input_files, output_file_name, pca_dim);

% clear all;
% close all;
%
% input_files.data_file='./data/Salinas_corrected.mat';
% input_files.gt_file='./data/Salinas_gt.mat';
% output_file_name='./results/Salinas_pca.mat';
% pca_dim=10;
% prepare_data_gt(input_files, output_file_name, pca_dim);


text:
function Q = ForwardStepBC2D(xin, yin, nxin, nyin, mapI, mapO, mapW, mapC, Q, time);

%  function [Q] = ForwardStepBC2D(xin, yin, nxin, nyin, mapI, mapO, mapW, mapC, Q, time);
% Purpose: Impose channel boundary conditions on 2D Euler equations on weak form

% Example is Mach ** 0.3 ** flow in wind tunnel
gamma = 1.4;

% extract conserved variables
rho = Q(:,:,1); rhou = Q(:,:,2); rhov = Q(:,:,3); Ener = Q(:,:,4);

% Inflow conditions -- uniform inflow
rhoin = gamma; uin = 3.0; vin = 0.0; pin = 1.0;
Ein = pin/(gamma-1.0) + 0.5*rhoin*(uin^2+vin^2);

rho(mapI) = rhoin; rhou(mapI) = rhoin*uin; rhov(mapI) = rhoin*vin; Ener(mapI) = Ein;

% Outflow conditions -- supersonic outflow ( do nothing )

% Wall conditions -- reflective, isothermal, i.e., n.u=0, T=T(t=0)
rhoW = rho(mapW); rhouW = rhou(mapW); rhovW = rhov(mapW);
nxW = nxin(mapW);   nyW = nyin(mapW);

% reverse flow in normal direction in ghost elements
rhou(mapW) = rhouW - 2*nxW.*(nxW.*rhouW + nyW.*rhovW);
rhov(mapW) = rhovW - 2*nyW.*(nxW.*rhouW + nyW.*rhovW);

% pack modified conserved variables
Q(:,:,1) = rho; Q(:,:,2) = rhou; Q(:,:,3) = rhov; Q(:,:,4) = Ener;
return

repo: Aero-matlab


text:
function S = tt_matrix2cellstr(M)
%TT_MATRIX2CELLSTR Convert a matrix to a cell array of strings.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

fmt = get(0,'FormatSpacing');
format compact
S = evalc('disp(M)');
set(0,'FormatSpacing',fmt)
S = textscan(S,'%s','delimiter','\n','whitespace','');
S = S{1};
end

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
%% UNO Rosenbrock
clc
clear
%Objective
obj = @(x) (1-x(1))^2 + 100 *(x(2)-x(1)^2)^2;
%Setup Options
opts = optiset('solver','m1qn3','display','iter');
%Build & Solve
Opt = opti('obj',obj,'ndec',2,'options',opts)
x0 = [0 0]';
[x,fval,exitflag,info]= solve(Opt,x0)

%%
% grad = Opt.nlprob.funcs.gradient;
grad = @(x) autoJac(obj,x);
opts.display = 2;
opts.nupdates = 5;
opts.maxtime = 0.5;
opts.iterfun = [];

[x,f,e,i,ii] = m1qn3(obj,grad,x0,opts);

%% 4x Rosenbrock [1,1,1,1]
clc
%Objective
fun = @(x) sum(100*(x(1:end-1).^2-x(2:end)).^2+(x(1:end-1)-1).^2);
grad = @(x) mklJac(fun,x);
x0 = [randn(1,100) 1];
opts.display = 2;
[x,fval,ef,iter] = m1qn3(fun,grad,x0,opts);

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
%QUADRC Trainable quadratic classifier
%
%   W = QUADRC(A,R,S)
%   W = A*QUADRC([],R,S)
%   W = A*QUADRC(R,S)
%
% INPUT
% 	A 		Dataset
%   R,S		0 <= R,S <= 1, regularization parameters (default: R = 0, S = 0)
%
% OUTPUT
% 	W     Quadratic Discriminant Classifier mapping
%
% DESCRIPTION
% Computation of the quadratic classifier between the classes of the dataset
% A based on different class covariances. R and S are regularization
% parameters used for finding the covariance matrix as
%
%		G = (1-R-S)*G + R*diag(diag(G)) + S*mean(diag(G))*eye(size(G,1))
%
% NOTE
% This routine differs from QDC; instead of using the densities, it is based
% on the class covariances. The multi-class problem is solved by multiple
% two-class quadratic discriminants. It is, thereby, a quadratic equivalent
% of FISHERC.
%
% SEE ALSO (PRTools Guide)
% MAPPINGS, DATASETS, FISHERC, NMC, NMSC, LDC, UDC, QDC

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

% $Id: quadrc.m,v 1.5 2010/02/08 15:29:48 duin Exp $

function w = quadrc(varargin)

  mapname = 'Quadr';
  argin = shiftargin(varargin,'scalar');
  argin = setdefaults(argin,[],0,0);

  if mapping_task(argin,'definition')

   w = define_mapping(argin,'untrained',mapname);

	elseif mapping_task(argin,'training')			% Train a mapping.

    [a,arg2,s] = deal(argin{:});
    [m,k,c] = getsize(a);

		islabtype(a,'crisp');
		isvaldfile(a,2,2); % at least 2 objects per class, 2 classes
		a = testdatasize(a,'features');
		a = setprior(a,getprior(a));
		r = arg2;

		if (min(classsizes(a)) < 2)
			error('Classes should contain more than one vector.')
		end

		if (c == 2)

			% 2-class case: calculate quadratic discriminant parameters.

			%p = getprior(a); pa = p(1); pb = p(2);

			JA = findnlab(a,1); JB = findnlab(a,2);

			ma = mean(a(JA,:)); mb = mean(a(JB,:));

			GA = covm(a(JA,:)); GB = covm(a(JB,:));
			GA = (1-r-s) * GA + r * diag(diag(GA)) + ...
															s * mean(diag(GA))*eye(size(GA,1));
			GB = (1-r-s) * GB + r * diag(diag(GB)) + ...
															s*mean(diag(GB))*eye(size(GB,1));
			DGA = det(GA); DGB = det(GB);

			GA = prinv(GA);
			GB = prinv(GB);

			par1 = 2*ma*GA-2*mb*GB; par2 = GB - GA;

			% If either covariance matrix is nearly singular, substitute FISHERC.
			% Otherwise construct the mapping.

			if (DGA <= 0) | (DGB <= 0)
				prwarning(1,'Covariance matrix nearly singular, regularization needed; using FISHERC instead')
				w = fisherc(a);
			else
				%par0 = (mb*GB*mb'-ma*GA*ma') + 2*log(pa/pb) + log(DGB) -log(DGA);
				par0 = (mb*GB*mb'-ma*GA*ma') + log(DGB) -log(DGA);
				w = prmapping(mfilename,'trained',{par0,par1',par2},getlablist(a),k,2);
				w = cnormc(w,a);
				w = setname(w,mapname);
				w = setcost(w,a);
			end

		else

			% For C > 2 classes, recursively call this function, using MCLASSC.

			pars = feval(mfilename,[],r,s);
			w = mclassc(a,pars);

		end

	else % Evaluation

		% Second argument is a trained mapping: test. Note that we can only
		% have a 2-class case here. W's output will be [D, -D], as the distance
		% of a sample to a class is the negative distance to the other class.

    [a,v] = deal(argin{1:2});
    [m,k] = size(a);
		pars = getdata(v);
		d = +sum((a*pars{3}).*a,2) + a*pars{2} + ones(m,1)*pars{1};
		w = setdat(a,[d, -d],v);

	end

return

function p = logdet(cova,covb)
%
% Numerically feasible computation of log(det(cova)/det(covb)))
% (seems not to help)

k = size(cova,1);

sa = mean(diag(cova));
sb = mean(diag(covb));

deta = det(cova./sa)
detb = det(covb./sb)

p = k*log(sa) - k*log(sb) + log(deta) - log(detb);

author: marianux

repo: ecg-kit

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


text:
function rr = ACF(rel_data, wins, up)
%ACF calculates the autocorrelation function, and finds the RR.
%
%	            ACF(rel_data, wins, up)
%
%	Inputs:
%       rel_data    .t  vector of times
%                   .v  vector of resp Sig values
%                   .fs sampling freq
%       up              universal parameters structure
%       wins        .t  vector of start times
%
%	Outputs:
%       rr          .t  vector of times of estimated RRs
%                   .v  vector of estimated RRs
%                   .f  vector of freqs of power spectrum
%                   .p  vector of powers of power spectrum
%

%% Setup

downsample_freq = up.paramSet.fft_resample_freq;    % it would be worth changing this - it changes the answer (try 1 Hz e.g.)
true_fs = rel_data.fs;

%% Cycle through windows
rr.t = mean([wins.t_start(:)' ; wins.t_end(:)']); rr.t = rr.t(:);
rr.v = nan(length(rr.t),1);

for win_no = 1 : length(wins.t_start)

    % extract relevant data
    rel_els = find(rel_data.t >= wins.t_start(win_no) & rel_data.t < wins.t_end(win_no));
    data.v = rel_data.v(rel_els);
    data.t = rel_data.t(rel_els);

    good_els = ~isnan(data.v);
    data.v = data.v(good_els);
    data.t = data.t(good_els);

    % Downsample
    data.filt.t = downsample(data.t, true_fs/downsample_freq);
    data.filt.v = decimate(data.v, true_fs/downsample_freq);
    data.filt.v = detrend(data.filt.v);

    data.filt.v = data.filt.v(:);
    data.filt.t = data.filt.t(:);

    % Find the cross correlation

    [r,lags] = xcorr(data.filt.v);

    % find the fft of the cross correlation

    % Find FFT
    WINLENGTH = length(r);      % data.filt is a column vector, non?
    NFFT = 2^nextpow2(WINLENGTH);           % Size of the FFT
    HAMMWIN = hamming(WINLENGTH);          % Gots me my Hamming window.  Mmm... Hamm.
    HAMMWIN = HAMMWIN(:);
    f_nyq = downsample_freq/2;                           % Nyquist frequency.  "It's a Nyquist thing."
    FREQS = f_nyq.*linspace(0, 1, NFFT/2+1);            % Array of correspondent FFT bin frequencies, in BR (RPM)
    % NB: the 60 stays regardless of window length because this is used to calculate the freq in bpm

    WINDATA = detrend(r);                      % Remove the LSE straight line from the data
    % used to be: WINDATA = detrend(data.filt.v, 'constant');                      % Remove the LSE straight line from the data
    WINDATA = WINDATA .* HAMMWIN;                           % Apply the Hamm
    myFFT = fft(WINDATA, NFFT);                             % Congratulations, madame - it's an FFT
    myFFT = myFFT(1 : NFFT/2 + 1);                          % Discrard the upper half, which reflects the lower half
    myFFT = 2.*abs(myFFT/NFFT);                             % Single-sided amplitude spectrum
    psdx = (1/(downsample_freq*NFFT)) * abs(myFFT).^2;
    psdx(2:end-1) = 2*psdx(2:end-1);

    data.power = 10*log10(psdx); data.power = data.power(:);
    data.freqs = FREQS; data.freqs = data.freqs(:);

    %% Find spectral peak
    [rr.v(win_no), rr.f{win_no}, rr.p{win_no}] = find_spectral_peak(data, up);

end

end

author: peterhcharlton

repo: RRest


text:
%SWIVELDATA Create look-up tables of swivel angle median and range
%
% Will save the file SwivelData.mat to the same location as this
% function. The median and range look-up tables are called phi_med and
% phi_range respectively. Angles are in radians, and non-reachable
% points are NaN. They are made with 1deg resolution, and are nxnxn in
% size, where n may be set, but default 101.
%
% Copyright (C) Bryan Moutrie, 2013-2014
% Licensed under the GNU Lesser General Public License
% see full file for full statement
%
% Known issues:
%  - The model for the shoulder range of motion seems to give some
%     outlier results, which causes in a large swivel angle range
%  - When interpolating to find a value, when near the surface of the
%     arm's reacahble workspace, a NaN may be returned though it is
%     reachable

% LICENSE STATEMENT:
%
% This file is part of pHRIWARE.
%
% pHRIWARE 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.
%
% pHRIWARE 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 Lesser General Public
% License along with pHRIWARE.  If not, see <http://www.gnu.org/licenses/>.

w = which('wrapToPi');
if isempty(w)
    error(pHRIWARE('error', ...
        ['The MATLAB wrapToPi function could not be detected.', ...
        'Replace it with your own version']));
end

n = 101;

hal = HAL();
AP = hal.AP;
reach = sum([AP{2} AP{3}]);

X = linspace(-reach,reach,n);
Y = linspace(-reach,reach,n);
Z = linspace(-reach,reach,n);
PHI = d2r*(-180:179);

phi_med = zeros(n,n,n,'single');
phi_range = zeros(n,n,n,'single');

Q1 = zeros(length(PHI), 7);
Q2 = zeros(length(PHI), 7);

for x = 1:n
    for y = 1:n
        for z = 1:n
            [~,~,Tu] = h2fsu(AP,[X(x),Y(y),Z(z)]',PHI);
            [Q1(:,1:3), Q2(:,1:3)] = gikine(AP{1},Tu);
            [~, u] = hal.reachable(Q1, Q2);

            valid = logical(~u');
%             v = find(valid);
%             if ~isempty(v)
%                 d = diff([v, v(end)+v(1)]);
%                 if ~all(d==1)
%                     longest = [0 0 0];
%                     count = 0;
%                     for i=1:length(d)
%                         if d(i) == 1
%                             count = count+1;
%                         else
%                             if count > longest(3)
%                                 longest = [i-count i-1 count];
%                             end
%                             count = 0;
%                         end
%                     end
%                     if ~isequal(longest,[0 0 0])
%                         valid = v(longest(1):longest(2));
%                     else
%                         valid = false(size(PHI));
%                     end
%                 else
%                     valid = v;
%                 end
%             end

                reachablePhi = PHI(valid);

                if ~isempty(reachablePhi)
                    phiM = max(reachablePhi);
                    phim = min(reachablePhi);

                    if phim == PHI(1) && phiM == PHI(end)
                        unreachablePhi = PHI(logical(~valid));
                        phiM = min(unreachablePhi) + 1*d2r + 2*pi;
                        phim = max(unreachablePhi) - 1*d2r;
                    end
                else
                    phiM = NaN;
                    phim = NaN;
                end

            phi_med(x,y,z) = (phiM+phim)/2;
            phi_range(x,y,z) = (phiM-phim)/2;
        end
    end
    fprintf('.');
    if ~mod(x,10), fprintf('\n'); end
end

phi_med = wrapToPi(phi_med);

dir = fileparts(which('swivelData.m'));
save([dir,'\SwivelData.mat'],'phi_med','phi_range');
fprintf('saved!\n');

author: Allopart

repo: rbpf-gmapping


text:
function f = abs(f, varargin)
%ABS   Absolute value of a CHEBTECH object.
%   ABS(F) returns the absolute value of F, where F is a CHEBTECH object with no
%   roots in [-1 1]. If ~isempty(roots(F)), then ABS(F) will return garbage
%   with no warning. F may be complex.

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

if ( isreal(f) || isreal(1i*f) )
    % Convert to values and then compute ABS().
    values = f.coeffs2vals(f.coeffs);
    values = abs(values);
    f.coeffs = f.vals2coeffs(values);
else
    f = compose(f, @abs, [], [], varargin{:});
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
filename='Bridge_quad_fine';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance';'perimeter'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
optimizer = 'IPOPT'; incrementFactor = 1;
filterType = 'P1';

nsteps = 5;
Vfrac_final = 0.2;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-2;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;


text:
%% DEMO_import_FEB_export_INP
% Below is a demonstration of how import a FEB file and subsequently export
% the geometry into an INP file.

%%

clear; close all; clc;

%%
% Plot settings
fontSize=15;
faceAlpha1=0.5;
faceAlpha2=0.5;
edgeColor=0.25*ones(1,3);
edgeWidth=1.5;
markerSize1=50;

%% Importing .feb file

%Set main folders
defaultFolder = fileparts(fileparts(mfilename('fullpath')));
pathName_FEB=fullfile(defaultFolder,'data','FEB'); %Where to load the FEB file
pathName_INP=fullfile(defaultFolder,'data','INP'); %Where to export the INP file

febFileNamePart='example_HEX_QUAD.feb';
febFileName=fullfile(pathName_FEB,febFileNamePart);
[febXML,nodeStruct,elementCell]=import_FEB(febFileName);

%% Plotting model

% Plotting the example model surfaces
hf1=cFigure;
title('Visualizing model','FontSize',fontSize);
xlabel('X','FontSize',fontSize);ylabel('Y','FontSize',fontSize); zlabel('Z','FontSize',fontSize);
hold on;

uniqueMaterialIndices=[];
for q=1:1:numel(elementCell)
    uniqueMaterialIndices=unique([uniqueMaterialIndices(:); elementCell{q}.E_mat(:)]);
    switch elementCell{q}.E_type
        case {'tri3', 'quad4'}
            F=elementCell{q}.E;
            V=nodeStruct.N;
            C=elementCell{q}.E_mat;
       case {'hex8', 'tet4'}
            [F,C]=element2faces(elementCell{q}.E,elementCell{q}.E_mat); %Creates faces and colors (e.g. stress) for patch based plotting
    end
    hp=patch('Faces',F,'Vertices',V,'EdgeColor','k','FaceColor','flat','Cdata',C,'FaceAlpha',0.8);
end

colormap(jet(numel(uniqueMaterialIndices))); hc=colorbar; caxis([min(uniqueMaterialIndices)-0.5 max(uniqueMaterialIndices)+0.5]);
axis equal; view(3); axis tight;  grid on; set(gca,'FontSize',fontSize);
camlight('headlight');
drawnow;

%% EXPORTING INP FILES FOR EACH ELEMENT TYPE

%You can change this example to do this for material type instead. Just use
%the material indices to select the elements from the lists. However the
%export_INP function can only handle 1 element type at a time at the moment

for q=1:1:numel(elementCell)

    inpFileNamepart=[febFileNamePart(1:end-4),'_',num2str(q),'.inp']; %filename for inp file
    inpFileName=fullfile(pathName_INP,inpFileNamepart);

    elementStruct=elementCell{q};
    %Setting appropriate element type line for ABAQUS. CHECK THESE!
    switch elementStruct.E_type
        case 'tri3'
            elementStruct.E_type='*ELEMENT, TYPE=STRI3, ELSET=PART-DEFAULT_1_EB1';
        case 'quad4'
            elementStruct.E_type='*ELEMENT, TYPE=S4R, ELSET=PART-DEFAULT_1_EB1';
        case 'tet4'
            elementStruct.E_type='*ELEMENT, TYPE=C3D4, ELSET=PART-DEFAULT_1_EB1';
        case 'hex8'
            elementStruct.E_type='*ELEMENT, TYPE=C3D8R, ELSET=PART-DEFAULT_1_EB1';
    end
    export_INP(elementStruct,nodeStruct,inpFileName);
end

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [idx, dist]=closestnode(node,p)
%
% [idx, dist]=closestnode(node,p)
%
% Find the closest point in a node list and return its index

% author: Qianqian Fang (q.fang at neu.edu)
%
% input:
%    node: each row is an N-D node coordinate
%    p: a given position in the same space
%
% output:
%    idx: the index of the position in the node list that has the shortest
%         Euclidean distance to the position p
%    dist: the distances between p and each node
%
%
% -- this function is part of brain2mesh toolbox (http://mcx.space/brain2mesh)
%    License: GPL v3 or later, see LICENSE.txt for details
%

dd=node-repmat(p,size(node,1),1);
[dist,idx]=min(sum(dd.*dd,2));

author: fangq


text:
function boolean = ExitConditionEvaluation(T_new,T_reference,x_intervals,y_intervals)

boolean = 0;
Tempcount = 0;

for x_index = 2:1:x_intervals-1
    for y_index = 2:1:y_intervals-1
         Tempcount = Tempcount + T_new(x_index,y_index);
    end
end

if(round(Tempcount/((x_intervals-1)*(y_intervals-1))) <= T_reference)
boolean = 1;
end

end


text:
function [VariableType,F_x,F_w,F_xw,h] = robust_classify_variables_newest(F,h,ops,w);

Dependency = iterateDependance( yalmip('monomtable') | yalmip('getdependence') |  yalmip('getdependenceUser'));
DependsOnw = find(any((Dependency(:,getvariables(w))),2));

h_variables = getvariables(h);
h_w = find(ismember(h_variables,DependsOnw));
if ~isempty(h_w)
    % If we have Q(x) + w*h(x) we can replace w*h(x) with t and have a
    % quadratic objective, and the linearly parameterized w*h(x) <= t
    % However, if h is quadratic in w such as x^2+2*x*w+w^2 we must move
    % everything to the cone epigraph, otherwise it is not convex in x
    if ~is(h,'compound') && max(degree(h,w)) == 1 && min(degree(h,w)) == 1
        base = getbase(h);
        h0 = base(1);
        base = base(2:end);base = base(:);
        sdpvar t

        F = [F,base(h_w(:))'*recover(h_variables(h_w)) <= t];
        base(h_w) = 0;
        h = h0 + base(:)'*recover(h_variables) + t;
    else
        sdpvar t
        F = [F, h <= t];
        h = t;
    end
    Dependency = iterateDependance(yalmip('monomtable') | yalmip('getdependence') | yalmip('getdependenceUser'));
    DependsOnw = find(any((Dependency(:,getvariables(w))),2));
end

DoesNotDependOnw = find(~any((Dependency(:,getvariables(w))),2));
[notused,x_variables] = find(Dependency(DoesNotDependOnw,:));

F_w = [];
F_x = [];
F_xw = [];
for i = 1:length(F)
    F_vars = getvariables(F(i));
    F_vars = find(any((Dependency(F_vars,:)),1));
    if all(ismember(F_vars,DependsOnw))
        F_w = F_w + F(i);
    elseif all(ismember(F_vars,DoesNotDependOnw))
       F_x = F_x + F(i);
    else
       F_xw = F_xw + F(i);
    end
end
ops.removeequalities = 0;
[F_x,failure,cause] = expandmodel(F_x,h,ops);
[F_w,failure,cause] = expandmodel(F_w,[],ops);
ops.expandbilinear = 1;
ops.reusemodel = 1; % Might be case x+norm(w)<1, norm(w)<1
[F_xw,failure,cause] = expandmodel(F_xw,h,ops,w);

w_variables = depends(F_w);
x_variables = unique([depends(F_x) depends(F_xw) depends(h)]);
x_variables = setdiff(x_variables,w_variables);

% After exanding the conic represntable, we have introduced new variables
Dependency = iterateDependance(yalmip('monomtable') | yalmip('getdependence') | yalmip('getdependenceUser'));

auxiliary = unique([yalmip('extvariables') yalmip('auxvariables')]);
if 1%~isempty(auxiliary)
    DependsOnw = find(any((Dependency(:,getvariables(w))),2));
    DoesNotDependOnw = find(~any((Dependency(:,getvariables(w))),2));

    temp = intersect(DependsOnw,x_variables);
    x_variables = setdiff(x_variables,DependsOnw);
    aux_with_w_dependence = temp;
else
    aux_with_w_dependence = [];
end

% aux_w_or_w = union(aux_with_w_dependence,w_variables);
% old_w_variables = [];
% while ~isequal(w_variables,old_w_variables);
%     old_w_variables = w_variables;
%     for i = 1:length(F_xw)
%         if all(ismember(depends(F_xw(i)),aux_w_or_w))
%             if ~any(ismember(depends(F_xw(i)),x_variables))
%                 new_w = intersect(depends(F_xw(i)),aux_w_or_w);
%                 w_variables = union(w_variables,new_w);
%                 aux_with_w_dependence = setdiff(aux_with_w_dependence,new_w);
%                 goon = 1;
%             end
%         end
%     end
%     aux_w_or_w = union(aux_with_w_dependence,w_variables);
% end

x = recover(x_variables);
w = recover(w_variables);

if ~isempty(F_xw)
    F_xw_scalar =  F_xw(find(is(F_xw,'elementwise') | is(F_xw,'equality')));
    F_xw_multi = F_xw - F_xw_scalar;
else
    F_xw_scalar = [];
    F_xw_multi = F_xw;
end

[MonomTable,Nonlinear] = yalmip('monomtable');
Dependency = yalmip('getdependenceUser');
evar = yalmip('extvariables');
if length(F_xw_scalar)>0
    % Optimize dependency graph
    X = sdpvar(F_xw_scalar);
    Xvar = getvariables(X);
    Xbase = getbase(X);Xbase = Xbase(:,2:end);
    for i = 1:size(Xbase,1)
        used = Xvar(find(Xbase(i,:)));
        if any(Nonlinear(used))
        used = find(any(MonomTable(used,:),1));
        end
        auxUsed = intersect(used,aux_with_w_dependence);

        if ~isempty(auxUsed)
            wUsed = intersect(used,w_variables);
            if ~isempty(wUsed)
                Dependency(auxUsed,wUsed) = 1;
            end
            eUsed = intersect(used,evar);
            if ~isempty(eUsed)
                Dependency(eUsed,auxUsed) = 1;
            end
        end
    end
end
if length(F_xw_multi) > 0
    for i = 1:length(F_xw_multi)
        used = getvariables(F_xw_multi(i));
        used = find(any(MonomTable(used,:),1));
        auxUsed = intersect((used),aux_with_w_dependence);
        wUsed = intersect((used),w_variables);
        if ~isempty(auxUsed) & ~isempty(wUsed)
            Dependency(auxUsed,wUsed) = 1;
        end
        eUsed = intersect((used),evar);
        if ~isempty(auxUsed) & ~isempty(eUsed)
            Dependency(eUsed,auxUsed) = 1;
        end
    end
end
UserDependency = yalmip('getdependenceUser');
fixed = find(any(UserDependency,2));
Dependency(fixed,:) = UserDependency(fixed,:);
Dependency = iterateDependance(Dependency);

VariableType.Graph = Dependency;
VariableType.x_variables = x_variables;
VariableType.w_variables = w_variables;
VariableType.aux_with_w_dependence = aux_with_w_dependence;

function Graph = iterateDependance(Graph)

Graph = Graph + speye(length(Graph));
Graph0 = double(Graph*Graph ~=0);
while ~isequal(Graph,Graph0)
    Graph = Graph0;
    Graph0 = double(Graph*Graph~=0);
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function fh = copyfig(fh)
%COPYFIG Create a copy of a figure, without changing the figure
%
% Examples:
%   fh_new = copyfig(fh_old)
%
% This function will create a copy of a figure, but not change the figure,
% as copyobj sometimes does, e.g. by changing legends.
%
% IN:
%    fh_old - The handle of the figure to be copied. Default: gcf.
%
% OUT:
%    fh_new - The handle of the created figure.

% Copyright (C) Oliver Woodford 2012, Yair Altman 2015

% 26/02/15: If temp dir is not writable, use the dest folder for temp
%           destination files (Javier Paredes)
% 15/04/15: Suppress warnings during copyobj (Dun Kirk comment on FEX page 2013-10-02)
% 09/09/18: Fix issue #252: Workaround for cases where copyobj() fails for any reason

    % Set the default
    if nargin == 0
        fh = gcf;
    end
    % Is there a legend?
    useCopyobj = isempty(findall(fh, 'Type', 'axes', 'Tag', 'legend'));
    if useCopyobj
        % Safe to copy using copyobj
        oldWarn = warning('off'); %Suppress warnings during copyobj (Dun Kirk comment on FEX page 2013-10-02)
        try
            fh = copyobj(fh, 0);
        catch
            % Fix issue #252: Workaround for cases where copyobj() fails for any reason
            useCopyobj = false;  % if copyobj() croaks, use file save/load below
        end
        warning(oldWarn);
    end
    if ~useCopyobj
        % copyobj will change the figure, so save and then load it instead
        tmp_nam = [tempname '.fig'];
        try
            % Ensure that the temp dir is writable (Javier Paredes 26/2/15)
            fid = fopen(tmp_nam,'w');
            fwrite(fid,1);
            fclose(fid);
            delete(tmp_nam);  % cleanup
        catch
            % Temp dir is not writable, so use the current folder
            [dummy,fname,fext] = fileparts(tmp_nam); %#ok<ASGLU>
            fpath = pwd;
            tmp_nam = fullfile(fpath,[fname fext]);
        end
        hgsave(fh, tmp_nam);
        fh = hgload(tmp_nam);
        delete(tmp_nam);
    end
end

author: altmany


text:
function yrange = GetAxesYRangeForStimPlot(hAxes)

if ~exist('hAxes','var') || isempty(hAxes)
    hAxes = gca;
end

ylim = get(hAxes, 'ylim');

% Fill up the ylim range with the stim by adding and
% subtracting a tiny amount so that the ylim isn't
% automatically changed
d = (1e-4)*(ylim(2)-ylim(1));
% d = 0;
yrange = [ylim(1)+d, ylim(2)-d];

author: BUNPC