text:
classdef L1_ADMM_YZ < handle
% Solver for various L1 minimization problems
% based on the 2011 paper by
% Junfeng Yang and Yin Zhang

properties
    % weight for the quadratic penalty term
    rho
     % verbosity
    verbose
    % Maximum number of ADMM iterations
    max_iterations
    % relative tolerance
    tolerance
    % additional scaling factor for x updates
    gamma
end % properties

properties(SetAccess=private)
    % The sensing operator
    A
    % The signals being recovered
    B
    % The solution representation vectors
    X
    % The dimension of signal space
    M
    % The dimension of representation space
    N
    % Details of execution (intermediate values)
    details
end % properties

% Following properties are meant for implementation only
properties(Access=private)
    X0
    Y0
    Z0
end

methods

function self = L1_ADMM_YZ(A, options)
    if nargin < 1
        error('A must be specified.');
    end
    if isa(A, 'spx.dict.Operator')
        self.A = A;
    elseif ismatrix(A)
        self.A = spx.dict.MatrixOperator(A);
    else
        error('Unsupported operator.');
    end
    [self.M, self.N] = size(self.A);
    if nargin < 2
        options = struct;
    end
    self.init_options();
    self.process_options(options);
end % function

function [X] = solve_bp(self, B, options)
    % Solves the problem min ||x||_1 s.t. Ax = b
    if nargin < 2
        error('B must be specified.');
    end
    if nargin < 3
        options = struct;
    end
    self.process_options(options);
    num_problems = size(B, 2);
    if self.verbose > 0
        fprintf('Solving l1 minimization problem: BP.\n');
    end
    self.init_details(num_problems);
    max_iterations = self.max_iterations;
    A  = self.A;
    X = zeros(self.N, num_problems);
    % import relevant functions from other modules
    import spx.opt.projections.proj_linf_ball;
    % iterate over problems
    for prob=1:num_problems
        tstart = tic;
        b = B(:, prob);
        b_max = norm(b, 'inf');
        terminated = false;
        % Check for zero solution condition
        if b_max < self.tolerance
            % There is no need to proceed further
            % zero solution is best solution
            self.details.terminated(prob) = true;
            self.details.iterations(prob) = 0;
            self.details.elapsed_times(prob) = toc(tstart);
            X(:, prob) = zeros(self.N, 1);
            continue;
        end
        if self.rho > 0
            % we will use user defined quadratic term weight
            rho = self.rho;
        else
            rho = mean(abs(b));
        end
        gamma = self.gamma;
        % scale the problem
        b = b / b_max;
        % Initialize solution
        x = A.adjoint(b);
        z = zeros(self.N, 1);
        y = zeros(self.M, 1);
        % compute current primary residual
        r_primal = A.apply(x) - b;
        for iter=1:max_iterations
            % check if we need to iterate further
            if terminated; break; end;
            % update z
            z_prev = z;
            z = proj_linf_ball(A.adjoint(y) + (x / rho));
            % update y
            y_prev = y;
            Az = A.apply(z);
            y = (Az - r_primal / rho);
            % update dual residual
            Aty = A.adjoint(y);
            r_dual =  Aty - z;
            % update x
            x_prev = x;
            x = x + (gamma * rho)*r_dual;
            % update primary residual
            r_primal = A.apply(x) - b;
            % primal objective
            primal_obj = sum(abs(x));
            % dual objective
            dual_obj = real(b' * y);
            terminated = self.check_termination(x, x_prev, y, z, ...
                r_primal, r_dual, primal_obj, dual_obj, iter, prob);
        end % iteration loop
        self.details.terminated(prob) = terminated;
        self.details.iterations(prob) = iter-1;
        self.details.elapsed_times(prob) = toc(tstart);
        % put the final solution back in result
        X(:, prob) = x * b_max;
    end % problem loop
end % function

function [X] = solve_bpic(self, B, delta, options)
    % Solves the problem min ||x||_1 s.t. ||Ax - b||_2 <= delta
    if nargin < 2
        error('B must be specified.');
    end
    if nargin < 3
        error('delta must be specified.');
    end
    if nargin < 4
        options = struct;
    end
    self.process_options(options);
    num_problems = size(B, 2);
    if self.verbose > 0
        fprintf('Solving constrained l1 minimization problem: BPIC.\n');
    end
    self.init_details(num_problems);
    max_iterations = self.max_iterations;
    A  = self.A;
    X = zeros(self.N, num_problems);
    % import relevant functions from other modules
    import spx.opt.projections.proj_linf_ball;
    import spx.opt.projections.proj_l2_ball;
    % iterate over problems
    for prob=1:num_problems
        tstart = tic;
        b = B(:, prob);
        b_max = norm(b, 'inf');
        terminated = false;
        % Check for zero solution condition
        if norm(b) < delta
            % There is no need to proceed further
            % zero solution is best solution
            self.details.terminated(prob) = true;
            self.details.iterations(prob) = 0;
            self.details.elapsed_times(prob) = toc(tstart);
            X(:, prob) = zeros(self.N, 1);
            continue;
        end
        if self.rho > 0
            % we will use user defined quadratic term weight
            rho = self.rho;
        else
            rho = mean(abs(b));
        end
        gamma = self.gamma;
        % scale the problem
        b = b / b_max;
        % Initialize solution
        x = A.adjoint(b);
        z = zeros(self.N, 1);
        y = zeros(self.M, 1);
        % compute current primary residual
        r_primal = A.apply(x) - b;
        delta_by_rho = delta / rho;
        for iter=1:max_iterations
            % check if we need to iterate further
            if terminated; break; end;
            % update z
            z_prev = z;
            z = proj_linf_ball(A.adjoint(y) + (x / rho));
            % update y
            y_prev = y;
            Az = A.apply(z);
            y = (Az - r_primal / rho);
            y = y - proj_l2_ball(y, delta_by_rho);
            % update dual residual
            Aty = A.adjoint(y);
            r_dual =  Aty - z;
            % update x
            x_prev = x;
            x = x + (gamma * rho)*r_dual;
            % update primary residual
            r_primal = A.apply(x) - b;
            % primal objective
            primal_obj = sum(abs(x));
            % dual objective
            dual_obj = real(b' * y) - delta * norm(y);
            terminated = self.check_termination(x, x_prev, y, z, ...
                r_primal, r_dual, primal_obj, dual_obj, iter, prob);
        end % iteration loop
        self.details.terminated(prob) = terminated;
        self.details.iterations(prob) = iter-1;
        self.details.elapsed_times(prob) = toc(tstart);
        % put the final solution back in result
        X(:, prob) = x * b_max;
    end % problem loop
end % function

function [X] = solve_bpdn_l2(self, B, mu, options)
    % Solves the problem min ||x||_1  + ||Ax - b||_2 / (2 mu)
    if nargin < 2
        error('B must be specified.');
    end
    if nargin < 3
        error('mu must be specified.');
    end
    if nargin < 4
        options = struct;
    end
    self.process_options(options);
    num_problems = size(B, 2);
    if self.verbose > 0
        fprintf('Solving l1-l2 unconstrained minimization problem: BPDN l2.\n');
    end

    self.init_details(num_problems);
    max_iterations = self.max_iterations;
    A  = self.A;
    X = zeros(self.N, num_problems);
    % import relevant functions from other modules
    import spx.opt.projections.proj_linf_ball;
    % iterate over problems
    for prob=1:num_problems
        tstart = tic;
        b = B(:, prob);
        %  Compute A' b
        Atb = A.adjoint(b);
        terminated = false;
        % Check for zero solution condition
        if norm(Atb, 'inf') <= mu
            % There is no need to proceed further
            % zero solution is best solution
            self.details.terminated(prob) = true;
            self.details.iterations(prob) = 0;
            self.details.elapsed_times(prob) = toc(tstart);
            X(:, prob) = zeros(self.N, 1);
            continue;
        end
        if self.rho > 0
            % we will use user defined quadratic term weight
            rho = self.rho;
        else
            rho = mean(abs(b));
        end
        gamma = self.gamma;
        % scale the problem
        b_max = norm(b, 'inf');
        b = b / b_max;
        mu = mu / b_max;
        % Initialize solution
        x = A.adjoint(b);
        z = zeros(self.N, 1);
        y = zeros(self.M, 1);
        % compute current primary residual
        r_primal = A.apply(x) - b;
        for iter=1:max_iterations
            % check if we need to iterate further
            if terminated; break; end;
            % update z
            z_prev = z;
            z = proj_linf_ball(A.adjoint(y) + (x / rho));
            % update y
            y_prev = y;
            Az = A.apply(z);
            y = (rho/(mu +  rho)) * (Az - r_primal / rho);
            % update dual residual
            Aty = A.adjoint(y);
            r_dual =  Aty - z;
            % update x
            x_prev = x;
            x = x + (gamma * rho)*r_dual;
            % update primary residual
            r_primal = A.apply(x) - b;
            % primal objective
            primal_obj = sum(abs(x)) + (1/(2*mu))*(r_primal'*r_primal);
            % dual objective
            dual_obj = real(b' * y) - (mu / 2) * (y' * y);
            terminated = self.check_termination(x, x_prev, y, z, ...
                r_primal, r_dual, primal_obj, dual_obj, iter, prob);
        end % iteration loop
        self.details.terminated(prob) = terminated;
        self.details.iterations(prob) = iter-1;
        self.details.elapsed_times(prob) = toc(tstart);
        % put the final solution back in result
        X(:, prob) = x * b_max;
    end % problem loop
end % function

function [X] = solve_bpdn_l1(self, B, nu, options)
    % Solves the problem min ||x||_1  + ||Ax - b||_1 / (nu)
end % function

end % methods

methods(Access=private)

function init_options(self)
    self.rho = 0;
    self.verbose = 0;
    self.max_iterations = 100;
    % self.eps_abs = 1e-4;
    self.tolerance = 1e-2;
    self.gamma = 1;
end % function

function process_options(self, options)
    % weight for the quadratic penalty term
    if isfield(options, 'rho')
        self.rho = options.rho;
    end
    % verbosity
    if isfield(options, 'verbose')
        self.verbose = options.verbose;
    end
    % Maximum number of ADMM iterations
    if isfield(options, 'max_iterations')
        self.max_iterations = options.max_iterations;
    end
    % gamma for primal variable update
    if isfield(options, 'gamma')
        self.gamma = options.gamma;
    end
    % relative tolerance
    if isfield(options, 'tolerance')
        self.tolerance = options.tolerance;
    end
end % function

function init_xyz(self, num_problems, options)
    if isfield(options, 'X0')
        self.X0 = X0;
    else
        self.X0 = zeros(self.N, num_problems);
    end
    if isfield(options, 'Y0')
        self.Y0 = Y0;
    else
        self.Y0 = zeros(self.N, num_problems);
    end
    if isfield(options, 'Z0')
        self.Z0 = Z0;
    else
        self.Z0 = zeros(self.N, num_problems);
    end
end % function

function init_details(self, num_problems)
    max_iterations = self.max_iterations;
    self.details.iterations = zeros(1, num_problems);
    self.details.terminated = zeros(1, num_problems);
    self.details.elapsed_times = zeros(1, num_problems);
    self.details.x_norms = zeros(max_iterations, num_problems);
    self.details.z_norms = zeros(max_iterations, num_problems);
    self.details.r_primal_norms = zeros(max_iterations, num_problems);
    self.details.r_dual_norms = zeros(max_iterations, num_problems);
    self.details.primal_objectives = zeros(max_iterations, num_problems);
    self.details.dual_objectives = zeros(max_iterations, num_problems);
end % function

function terminate = check_termination(self, x, x_prev, y, z, ...
    r_primal, r_dual, primal_obj, dual_obj, iter, prob)
    terminate = false;
    % norm of primal variable x
    x_norm = norm(x);
    % norm of change in x
    x_diff_norm = norm(x - x_prev);
    % relative change in x
    x_rel_change = x_diff_norm / x_norm;
    % norm of dual variable z
    z_norm = norm(z);
    % primal residual norm
    r_primal_norm  = norm(r_primal);
    % dual residual norm
    r_dual_norm = norm(r_dual);
    % relative change in dual norm
    r_dual_rel_change = r_dual_norm / z_norm;
    % gap in objectives
    gap = abs(primal_obj  - dual_obj);
    relative_gap = gap / abs(primal_obj);

    self.details.x_norms(iter, prob) = x_norm;
    self.details.z_norms(iter, prob) = z_norm;
    self.details.r_primal_norms(iter, prob) = r_primal_norm;
    self.details.r_dual_norms(iter, prob) = r_dual_norm;
    self.details.primal_objectives(iter, prob) = primal_obj;
    self.details.dual_objectives(iter, prob) = dual_obj;

    if self.verbose > 1
        fprintf('Objective primal: %.4f, dual: %.4f\n', primal_obj, dual_obj);
    end

    tolerance = self.tolerance;
    if x_rel_change < tolerance
        terminate = true;
    end
    if x_rel_change >= tolerance*1.1
        % we will continue for a while
        return;
    end
    if relative_gap < tolerance
        terminate = true;
    end
    if r_dual_rel_change < tolerance
        terminate = true;
    end
end % function

end % methods

end % class

author: indigits

repo: sparse-plex


text:
function [ params ] = setEdgeBoxesParamsFromConfig( config)

%%Loading Default Parameters here
%(1) main parameters
%   .name           - [] target filename (if specified return is 1)
%   .alpha          - [.65] step size of sliding window search
%   .beta           - [.75] nms threshold for object proposals
%   .minScore       - [.01] min score of boxes to detect
%   .maxBoxes       - [1e4] max number of boxes to detect

config.params.name=[];
if isfield(config.params, 'alpha')
	config.params.alpha=[config.params.alpha];
end
if isfield(config.params, 'beta')
	config.params.beta=[config.params.beta];
end
if isfield(config.params, 'minScore')
	config.params.minScore=[config.params.minScore];
end
if isfield(config.params, 'maxBoxes')
	config.params.maxBoxes=[config.params.maxBoxes];
end

%(2)  additional parameters, safe to ignore and leave at default vals

%   .edgeMinMag     - [.1] increase to trade off accuracy for speed
%   .edgeMergeThr   - [.5] increase to trade off accuracy for speed
%   .clusterMinMag  - [.5] increase to trade off accuracy for speed
%   .maxAspectRatio - [3] max aspect ratio of boxes
%   .minBoxArea     - [1000] minimum area of boxes
%   .gamma          - [2] affinity sensitivity, see equation 1 in paper
%   .kappa          - [1.5] scale sensitivity, see equation 3 in paper

if isfield(config.params, 'edgeMinMag')
	config.params.edgeMinMag = [config.params.edgeMinMag];
end
if isfield(config.params, 'edgeMergeThr')
	config.params.edgeMergeThr = [config.params.edgeMergeThr];
end
if isfield(config.params, 'clusterMinMag')
	config.params.clusterMinMag  = [config.params.clusterMinMag];
end
if isfield(config.params, 'maxAspectRatio')
	config.params.maxAspectRatio = [config.params.maxAspectRatio];
end
if isfield(config.params, 'minBoxArea')
	config.params.minBoxArea= [config.params.minBoxArea] ;
end
if isfield(config.params, 'gamma')
	config.params.gamma=[config.params.gamma];
end
if isfield(config.params, 'kappa')
	config.params.kappa=[config.params.kappa];
end

params = config.params;

end

author: Cloud-CV

repo: object-proposals


text:
function [spectrum,ntaper,freqoi] = ft_specest_mtmfft(dat, time, varargin)

% FT_SPECEST_MTMFFT computes a fast Fourier transform using multitapering with
% multiple tapers from the DPSS sequence or using a variety of single tapers.
%
% Use as
%   [spectrum,ntaper,freqoi] = ft_specest_mtmfft(dat,time...)
% where
%   dat        = matrix of chan*sample
%   time       = vector, containing time in seconds for each sample
%   spectrum   = matrix of taper*chan*freqoi of fourier coefficients
%   ntaper     = vector containing number of tapers per element of freqoi
%   freqoi     = vector of frequencies in spectrum
%
% Optional arguments should be specified in key-value pairs and can include
%   taper      = 'dpss', 'hanning' or many others, see WINDOW (default = 'dpss')
%   pad        = number, total length of data after zero padding (in seconds)
%   padtype    = string, indicating type of padding to be used (see ft_preproc_padding, default: zero)
%   freqoi     = vector, containing frequencies of interest
%   tapsmofrq  = the amount of spectral smoothing through multi-tapering. Note: 4 Hz smoothing means plus-minus 4 Hz, i.e. a 8 Hz smoothing box
%   dimord     = 'tap_chan_freq' (default) or 'chan_time_freqtap' for memory efficiency (only when use variable number slepian tapers)
%   polyorder  = number, the order of the polynomial to fitted to and removed from the data prior to the fourier transform (default = 0 -> remove DC-component)
%   taperopt   = additional taper options to be used in the WINDOW function, see WINDOW
%   verbose    = output progress to console (0 or 1, default 1)
%
% See also FT_FREQANALYSIS, FT_SPECEST_MTMCONVOL, FT_SPECEST_TFR, FT_SPECEST_HILBERT, FT_SPECEST_WAVELET

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

% these are for speeding up computation of tapers on subsequent calls
persistent previous_argin previous_tap

% get the optional input arguments
taper     = ft_getopt(varargin, 'taper'); if isempty(taper), ft_error('You must specify a taper'); end
pad       = ft_getopt(varargin, 'pad');
padtype   = ft_getopt(varargin, 'padtype', 'zero');
freqoi    = ft_getopt(varargin, 'freqoi', 'all');
tapsmofrq = ft_getopt(varargin, 'tapsmofrq');
dimord    = ft_getopt(varargin, 'dimord', 'tap_chan_freq');
fbopt     = ft_getopt(varargin, 'feedback');
verbose   = ft_getopt(varargin, 'verbose', true);
polyorder = ft_getopt(varargin, 'polyorder', 0);
tapopt    = ft_getopt(varargin, 'taperopt');

if isempty(fbopt)
  fbopt.i = 1;
  fbopt.n = 1;
end

% throw errors for required input
if isempty(tapsmofrq) && (strcmp(taper, 'dpss') || strcmp(taper, 'sine'))
  ft_error('you need to specify tapsmofrq when using dpss or sine tapers')
end

% this does not work on integer data
dat = cast(dat, 'double');

% Set n's
[nchan,ndatsample] = size(dat);

% This does not work on integer data
if ~isa(dat, 'double') && ~isa(dat, 'single')
  dat = cast(dat, 'double');
end

% Remove polynomial fit from the data -> default is demeaning
if polyorder >= 0
  dat = ft_preproc_polyremoval(dat, polyorder, 1, ndatsample);
end

% Determine fsample and set total time-length of data
fsample = 1./mean(diff(time));
dattime = ndatsample / fsample; % total time in seconds of input data

% Zero padding
if round(pad * fsample) < ndatsample
  ft_error('the padding that you specified is shorter than the data');
end
if isempty(pad) % if no padding is specified padding is equal to current data length
  pad = dattime;
end
postpad    = ceil((pad - dattime) * fsample);
endnsample = round(pad * fsample);  % total number of samples of padded data
endtime    = pad;                   % total time in seconds of padded data

% Set freqboi and freqoi
freqoiinput = freqoi;
if isnumeric(freqoi) % if input is a vector
  freqboi   = round(freqoi ./ (fsample ./ endnsample)) + 1; % is equivalent to: round(freqoi .* endtime) + 1;
  freqboi   = unique(freqboi);
  freqoi    = (freqboi-1) ./ endtime; % boi - 1 because 0 Hz is included in fourier output
elseif strcmp(freqoi,'all') % if input was 'all'
  freqboilim = round([0 fsample/2] ./ (fsample ./ endnsample)) + 1;
  freqboi    = freqboilim(1):1:freqboilim(2);
  freqoi     = (freqboi-1) ./ endtime;
end
nfreqboi = length(freqboi);
nfreqoi  = length(freqoi);
if (strcmp(taper, 'dpss') || strcmp(taper, 'sine')) && numel(tapsmofrq)~=1 && (numel(tapsmofrq)~=nfreqoi)
  ft_error('tapsmofrq needs to contain a smoothing parameter for every frequency when requesting variable number of slepian tapers')
end

% throw a warning if input freqoi is different from output freqoi
if isnumeric(freqoiinput)
  % check whether padding is appropriate for the requested frequency resolution
  rayl = 1/endtime;
  if any(rem(freqoiinput,rayl)) % not always the case when they mismatch
    ft_warning('padding not sufficient for requested frequency resolution, for more information please see the FAQs on www.ru.nl/neuroimaging/fieldtrip');
  end
  if numel(freqoiinput) ~= numel(freqoi) % freqoi will not contain double frequency bins when requested
    ft_warning('output frequencies are different from input frequencies, multiples of the same bin were requested but not given');
  else
    if any(abs(freqoiinput-freqoi) >= eps*1e6)
      ft_warning('output frequencies are different from input frequencies');
    end
  end
end

% determine whether tapers need to be recomputed
current_argin = {time, postpad, taper, tapsmofrq, freqoi, tapopt, dimord}; % reasoning: if time and postpad are equal, it's the same length trial, if the rest is equal then the requested output is equal
if isequal(current_argin, previous_argin)
  % don't recompute tapers
  tap = previous_tap;

else
  % recompute tapers
  switch taper

    case 'dpss'
      if numel(tapsmofrq)==1
        % create a sequence of DPSS tapers, ensure that the input arguments are double precision
        tap = double_dpss(ndatsample,ndatsample*(tapsmofrq./fsample))';
        % remove the last taper because the last slepian taper is always messy
        tap = tap(1:(end-1), :);

        % give error/warning about number of tapers
        if isempty(tap)
          ft_error('datalength to short for specified smoothing\ndatalength: %.3f s, smoothing: %.3f Hz, minimum smoothing: %.3f Hz',ndatsample/fsample,tapsmofrq,fsample/ndatsample);
        elseif size(tap,1) == 1
          ft_warning('using only one taper for specified smoothing');
        end
      elseif numel(tapsmofrq)>1
        tap = cell(1,nfreqoi);
        for ifreqoi = 1:nfreqoi
          % create a sequence of DPSS tapers, ensure that the input arguments are double precision
          currtap = double_dpss(ndatsample, ndatsample .* (tapsmofrq(ifreqoi) ./ fsample))';
          % remove the last taper because the last slepian taper is always messy
          currtap = currtap(1:(end-1), :);

          % give error/warning about number of tapers
          if isempty(currtap)
            ft_error('%.3f Hz: datalength to short for specified smoothing\ndatalength: %.3f s, smoothing: %.3f Hz, minimum smoothing: %.3f Hz',freqoi(ifreqoi), ndatsample/fsample,tapsmofrq(ifreqoi),fsample/ndatsample(ifreqoi));
          elseif size(currtap,1) == 1
            disp([num2str(freqoi(ifreqoi)) ' Hz: WARNING: using only one taper for specified smoothing'])
          end
          tap{ifreqoi} = currtap;
        end
      end

    case 'sine'
      if numel(tapsmofrq)==1
        % create a sequence of sine tapers,
        tap = sine_taper(ndatsample, ndatsample*(tapsmofrq./fsample))';
        % remove the last taper
        tap = tap(1:(end-1), :);

        % give error/warning about number of tapers
        if isempty(tap)
          ft_error('datalength to short for specified smoothing\ndatalength: %.3f s, smoothing: %.3f Hz, minimum smoothing: %.3f Hz',ndatsample/fsample,tapsmofrq,fsample/ndatsample);
        elseif size(tap,1) == 1
          ft_warning('using only one taper for specified smoothing');
        end
      elseif numel(tapsmofrq)>1
        tap = cell(1,nfreqoi);
        for ifreqoi = 1:nfreqoi
          % create a sequence of sine tapers
          currtap = sine_taper(ndatsample, ndatsample .* (tapsmofrq(ifreqoi) ./ fsample))';
          % remove the last taper because the last slepian taper is always messy
          currtap = currtap(1:(end-1), :);

          % give error/warning about number of tapers
          if isempty(currtap)
            ft_error('%.3f Hz: datalength to short for specified smoothing\ndatalength: %.3f s, smoothing: %.3f Hz, minimum smoothing: %.3f Hz',freqoi(ifreqoi), ndatsample/fsample,tapsmofrq(ifreqoi),fsample/ndatsample(ifreqoi));
          elseif size(currtap,1) == 1
            disp([num2str(freqoi(ifreqoi)) ' Hz: WARNING: using only one taper for specified smoothing'])
          end
          tap{ifreqoi} = currtap;
        end
      end

    case 'sine_old'
      % to provide compatibility with the tapers being scaled (which was default
      % behavior prior to 29apr2011) yet this gave different magnitude of power
      % when comparing with slepian multi tapers
      tap = sine_taper_scaled(ndatsample, ndatsample*(tapsmofrq./fsample))';
      tap = tap(1:(end-1), :); % remove the last taper

    case 'alpha'
      ft_error('not yet implemented');

    case 'hanning'
      tap = hanning(ndatsample)';
      tap = tap./norm(tap, 'fro');

    otherwise
      % create the taper and ensure that it is normalized
      if isempty(tapopt) % some windowing functions don't support nargin>1, and window.m doesn't check it
        tap = window(taper, ndatsample)';
      else
        tap = window(taper, ndatsample, tapopt)';
      end
      tap = tap ./ norm(tap,'fro');

  end % switch taper
end % isequal currargin

% set ntaper
if ~((strcmp(taper,'dpss') || strcmp(taper,'sine')) && numel(tapsmofrq)>1) % variable number of slepian tapers not requested
  ntaper = repmat(size(tap,1),nfreqoi,1);
else % variable number of slepian tapers requested
  ntaper = cellfun(@size,tap,repmat({1},[1 nfreqoi]));
end

% determine phase-shift so that for all frequencies angle(t=0) = 0
timedelay = time(1);
if timedelay ~= 0
  angletransform = complex(zeros(1,nfreqoi));
  for ifreqoi = 1:nfreqoi
    missedsamples = round(timedelay * fsample);
    % determine angle of freqoi if oscillation started at 0
    % the angle of wavelet(cos,sin) = 0 at the first point of a cycle, with sine being in upgoing flank, which is the same convention as in mtmconvol
    anglein = (missedsamples) .* ((2.*pi./fsample) .* freqoi(ifreqoi));
    coswav  = cos(anglein);
    sinwav  = sin(anglein);
    angletransform(ifreqoi) = atan2(sinwav, coswav);
  end
  angletransform = repmat(angletransform,[nchan,1]);
end

% compute fft
if ~((strcmp(taper,'dpss') || strcmp(taper,'sine')) && numel(tapsmofrq)>1) % ariable number of slepian tapers not requested
  str = sprintf('nfft: %d samples, datalength: %d samples, %d tapers',endnsample,ndatsample,ntaper(1));
  [st, cws] = dbstack;
  if length(st)>1 && strcmp(st(2).name, 'ft_freqanalysis')
    % specest_mtmfft has been called by ft_freqanalysis, meaning that ft_progress has been initialised
    ft_progress(fbopt.i./fbopt.n, ['processing trial %d/%d ',str,'\n'], fbopt.i, fbopt.n);
  elseif verbose
    fprintf([str, '\n']);
  end
  spectrum = cell(ntaper(1),1);

  for itap = 1:ntaper(1)
    dum = fft(ft_preproc_padding(bsxfun(@times,dat,tap(itap,:)), padtype, 0, postpad),[], 2);
    dum = dum(:,freqboi);
    % phase-shift according to above angles
    if timedelay ~= 0
      dum = dum .* exp(-1i*angletransform);
    end
    dum = dum .* sqrt(2 ./ endnsample);
    spectrum{itap} = dum;
  end

  spectrum = reshape(vertcat(spectrum{:}),[nchan ntaper(1) nfreqboi]); % collecting in a cell-array and later reshaping provides significant speedups
  spectrum = permute(spectrum, [2 1 3]);

else % variable number of slepian tapers requested
  switch dimord

    case 'tap_chan_freq' % default
      % start fft'ing
      spectrum = complex(NaN([max(ntaper) nchan nfreqoi]));
      for ifreqoi = 1:nfreqoi
        str = sprintf('nfft: %d samples, datalength: %d samples, frequency %d (%.2f Hz), %d tapers',endnsample,ndatsample,ifreqoi,freqoi(ifreqoi),ntaper(ifreqoi));
        [st, cws] = dbstack;
        if length(st)>1 && strcmp(st(2).name, 'ft_freqanalysis') && verbose
          % specest_mtmconvol has been called by ft_freqanalysis, meaning that ft_progress has been initialised
          ft_progress(fbopt.i./fbopt.n, ['processing trial %d, ',str,'\n'], fbopt.i);
        elseif verbose
          fprintf([str, '\n']);
        end
        for itap = 1:ntaper(ifreqoi)

          dum = fft(ft_preproc_padding(bsxfun(@times,dat,tap{ifreqoi}(itap,:)), padtype, 0, postpad), [], 2);

          dum = dum(:,freqboi(ifreqoi));
          % phase-shift according to above angles
          if timedelay ~= 0
            dum = dum .* exp(-1i*angletransform(:,ifreqoi));
          end
          dum = dum .* sqrt(2 ./ endnsample);
          spectrum(itap,:,ifreqoi) = dum;
        end
      end % for nfreqoi

    case 'chan_freqtap' % memory efficient representation
      % create tapfreqind
      freqtapind = cell(1,nfreqoi);
      tempntaper = [0; cumsum(ntaper(:))];
      for ifreqoi = 1:nfreqoi
        freqtapind{ifreqoi} = tempntaper(ifreqoi)+1:tempntaper(ifreqoi+1);
      end

      % start fft'ing
      spectrum = complex(zeros([nchan sum(ntaper)]));
      for ifreqoi = 1:nfreqoi
        str = sprintf('nfft: %d samples, datalength: %d samples, frequency %d (%.2f Hz), %d tapers',endnsample,ndatsample,ifreqoi,freqoi(ifreqoi),ntaper(ifreqoi));
        [st, cws] = dbstack;
        if length(st)>1 && strcmp(st(2).name, 'ft_freqanalysis') && verbose
          % specest_mtmconvol has been called by ft_freqanalysis, meaning that ft_progress has been initialised
          ft_progress(fbopt.i./fbopt.n, ['processing trial %d, ',str,'\n'], fbopt.i);
        elseif verbose
          fprintf([str, '\n']);
        end
        for itap = 1:ntaper(ifreqoi)

          dum = fft(ft_preproc_padding(bsxfun(@times,dat,tap{ifreqoi}(itap,:)), padtype, 0, postpad), [], 2);

          dum = dum(:,freqboi(ifreqoi));
          % phase-shift according to above angles
          if timedelay ~= 0
            dum = dum .* exp(-1i*angletransform(:,ifreqoi));
          end
          dum = dum .* sqrt(2 ./ endnsample);
          spectrum(:,freqtapind{ifreqoi}(itap)) = dum;
        end
      end % for nfreqoi
  end % switch dimord
end

% remember the current input arguments, so that they can be
% reused on a subsequent call in case the same input argument is given
previous_argin = current_argin;
previous_tap   = tap;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION ensure that the first two input arguments are of double
% precision this prevents an instability (bug) in the computation of the
% tapers for MATLAB 6.5 and 7.0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [tap] = double_dpss(a, b, varargin)
tap = dpss(double(a), double(b), varargin{:});

author: spm


text:
function pass = test_sum3(pref)
% Test file for @chebfun3t/sum3.

% Obtain preferences.
if ( nargin == 0 )
    pref = chebfunpref();
end
tol = 1e1 * pref.cheb3Prefs.chebfun3eps;

% Check empty case
pass(1) = isempty(sum3(chebfun3t()));

% Check constant
f = chebfun3t(@(x,y,z) 1);
pass(2) = abs(sum3(f) - 8) < tol;

% Runge function
f = chebfun3t(@(x,y,z) 1./(1+x.^2+y.^2+z.^2));
pass(3) = abs(sum3(f) - 4.28685406230184188268) < tol;

% Different domains
f = chebfun3t(@(x,y,z) x, [0, 1, 0, 2, 0, 3]);
pass(4) = abs(sum3(f) - 3) < tol;

f = chebfun3t(@(x,y,z) y, [0, 1, 0, 2, 0, 3]);
pass(5) = abs(sum3(f) - 6) < tol;

f = chebfun3t(@(x,y,z) z, [0, 1, 0, 2, 0, 3]);
pass(6) = abs(sum3(f) - 9) < tol;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function C = spm_dctmtx(N,K,n,f)
% Creates basis functions for Discrete Cosine Transform.
% FORMAT C = spm_dctmtx(N,K,n)
%     OR C = spm_dctmtx(N,K)
%     OR D = spm_dctmtx(N,K,n,'diff')
%     OR D = spm_dctmtx(N,K,'diff')
% N - dimension
% K - order
% n - optional points to sample
%____________________________________________________________________________
% spm_dctmtx creates a matrix for the first few basis functions of a one
% dimensional discrete cosine transform.
% With the 'diff' argument, spm_dctmtx produces the derivatives of the
% DCT.
%
% See:    Fundamentals of Digital Image Processing (p 150-154).
%         Anil K. Jain 1989.
%____________________________________________________________________________
% @(#)spm_dctmtx.m	2.1 John Ashburner MRCCU/FIL 01/08/28

d = 0;

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

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

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

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

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function C=cutcategories(t,j)
%CUTCATEGORIES Categories used for branches in decision tree.
%   C=CUTCATEGORIES(T) returns an N-by-2 cell array C of the categories
%   used at branches in the decision tree T, where N is the number of
%   nodes in the tree.  For each branch node K based on a categorical
%   predictor variable Z, the left child is chosen if Z is among the
%   categories listed in C{K,1}, and the right child is chosen if Z is
%   among those listed in C{K,2}.  Both columns of C are empty for branch
%   nodes based on continuous predictors and for non-branch (leaf) nodes.
%
%   C=CUTCATEGORIES(T,J) takes an array J of node numbers and returns the
%   categories for the specified nodes.
%
%   See also CLASSREGTREE, CLASSREGTREE/CUTVAR, CLASSREGTREE/CUTPOINT, CLASSREGTREE/CUTTYPE.

%   Copyright 2006-2007 The MathWorks, Inc.
%   $Revision: 1.1.6.6 $  $Date: 2009/01/08 18:59:32 $

if nargin>=2 && ~validatenodes(t,j)
    error('stats:classregtree:cutcategories:InvalidNode',...
          'J must be an array of node numbers or a logical array of the proper size.');
end

% Get variable numbers
if nargin<2
    n = t.var;
    cut = t.cut;
else
    n = t.var(j,:);
    cut = t.cut(j);
end

% Set up array for output
C = repmat({[] []},numel(n),1);

% Fetch category list pairs for categorical branches
catvar = find(n<0);
if ~isempty(catvar)
    for i=1:length(catvar)
        C(catvar(i),:) = cut{catvar(i)};
    end
end


text:
function cameras = loadcameradata(dataDir, idx)
%LOADCAMERADATA: Load the dino data
%
%   CAMERAS = LOADCAMERADATA() loads the dinosaur data and returns a
%   structure array containing the camera definitions. Each camera contains
%   the image, internal calibration and external calibration.
%
%   CAMERAS = LOADCAMERADATA(IDX) loads only the specified file indices.
%
%   Example:
%   >> cameras = loadcameradata(1:3);
%   >> showcamera(cameras)
%
%   See also: SHOWCAMERA

%   Copyright 2005-2009 The MathWorks, Inc.
%   $Revision: 1.0 $    $Date: 2006/06/30 00:00:00 $

if nargin<2
    idx = 1:36;
end

cameras = struct( ...
    'Image', {}, ...
    'P', {}, ...
    'K', {}, ...
    'R', {}, ...
    'T', {}, ...
    'Silhouette', {} );

%% First, import the camera Pmatrices
rawP = load( fullfile( dataDir, 'dino_Ps.mat') );

%% Now loop through loading the images
tmwMultiWaitbar('Loading images',0);
for ii=idx(:)'
    % We try both JPG and PPM extensions, trying JPEG first since it is
    % the faster to load
    filename = fullfile( dataDir, sprintf( 'viff.%03d.jpg', ii ) );
    if exist( filename, 'file' )~=2
        % Try PPM
        filename = fullfile( dataDir, sprintf( 'viff.%03d.ppm', ii ) );
        if exist( filename, 'file' )~=2
            % Failed
            error( 'SpaceCarving:ImageNotFound', ...
                'Could not find image %d (''viff.%03d.jpg/.ppm'')', ...
                ii, ii );
        end
    end

    [K,R,t] = spacecarving.decomposeP(rawP.P{ii});
    cameras(ii).rawP = rawP.P{ii};
    cameras(ii).P = rawP.P{ii};
    cameras(ii).K = K/K(3,3);
    cameras(ii).R = R;
    cameras(ii).T = -R'*t;
    cameras(ii).Image = imread( filename );
    cameras(ii).Silhouette = [];
    tmwMultiWaitbar('Loading images',ii/max(idx));
end
tmwMultiWaitbar('Loading images','close');

author: Sable

repo: mcbench-benchmarks


text:
function K = sdlfmaXsdlfmKernCompute(sdlfmaKern1, sdlfmKern2, t1, t2, covIC)

% SDLFMAXSDLFMKERNCOMPUTE Cross kernel between a SDLFMA and a SDLFM kernels.
% FORMAT
% DESC computes cross kernel terms between the accel. and the position of
% two switching dynamical LFM kernels for the multiple output kernel.
% ARG sdlfmaKern1 : the kernel structure associated with the accel. of the
% first SDLFM kernel.
% ARG sdlfmKern2 : the kernel structure associated with the position of the
% second SDLFM kernel.
% ARG t : inputs for which kernel is to be computed.
% ARG covIC : covariance for the initial conditions
% RETURN K : block of values from kernel matrix.
%
% FORMAT
% DESC computes cross kernel terms between the accel. and the position of
% two SDLFM kernels for the multiple output kernel.
% ARG sdlfmaKern1 : the kernel structure associated with the accel. of
% the first SDLFM kernel.
% ARG sdlfmKern2 : the kernel structure associated with the position of
% the second SDLFM kernel.
% ARG t1 : row inputs for which kernel is to be computed.
% ARG t2 : column inputs for which kernel is to be computed.
% ARG covIC : covariance for the initial conditions
% RETURN K : block of values from kernel matrix.
%
% SEEALSO : sdlfmaKernParamInit, sdlfmaKernCompute, sdlfmKernParamInit
%
% COPYRIGHT : Mauricio A. Alvarez, 2010

% KERN

if nargin < 4
    t2 = t1;
end

K = sdlfmXsdlfmKernCompute(sdlfmaKern1, sdlfmKern2, t1, t2, covIC, 'AccelPos');

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [g] = mg2g(mg)
% Convert mass from milligrams to grams.
% Chad Greene 2012
g = mg*0.001 ;


text:
function biosigpathlast

str2doublepath = fileparts( which('str2double') );
sopenpath      = fileparts( which('sopen') );
if strcmp(str2doublepath,sopenpath)
    rmpath(str2doublepath);
    str2doublepat2 = fileparts( which('str2double') );
    addpath(str2doublepath,'-begin');
    addpath(str2doublepat2,'-begin');
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function C = ttimes(A,B)
% TTIMES Tropical multiplication (max-plus algebra)
%
% C = TTIMES(A,B) Computes the tropical multiplication of the matrices A
% and B, AB(i,j) = max(A(i,:) + B(:,j)');
%
% See also tplus

n = size(A,1);
m = size(B,2);
C = reshape(max(kron(ones(m,1),A)+kron(B',ones(n,1)),[],2),n,m);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function res = lt(a,b)
%GE           Implements  a < b  for Taylor, compares only a.x and b.x
%

% written  05/21/09     S.M. Rump
%

  if ~isa(a,'taylor')
    res = ( a<reshape(b.t(1,:),b.size) );
  elseif ~isa(b,'taylor')
    res = ( reshape(a.t(1,:),a.size)<b );
  else
    res = ( reshape(a.t(1,:),a.size)<reshape(b.t(1,:),b.size) );
  end

repo: OpenMAS


text:
function Y = vl_nndwt2( X, dzdy, varargin )
% dwt multi_channel
% only support haart and db2 for 1-level transformation

opts.padding = 0 ;
opts.wavename = 'haart';
opts = vl_argparse(opts, varargin, 'nonrecursive') ;
padding = opts.padding;

if  nargin <= 1 || isempty(dzdy)
    sz = size(X);
    X = X/2;
    if size(X, 3) == 1
        sz(3) = 1;
    end
    if size(X, 4) == 1
        sz(4) = 1;
    end
            im_c1 = X(1:2:end, 1:2:end, :, :);
            im_c2 = X(1:2:end, 2:2:end, :, :);
            im_c3 = X(2:2:end, 1:2:end, :, :);
            im_c4 = X(2:2:end, 2:2:end, :, :);
            Y = zeros([sz(1)/2 sz(2)/2 sz(3)*4 sz(4)], 'like', X);
            Y(:,:,1:sz(3),:)           = im_c1 + im_c2 + im_c3 + im_c4;
            Y(:,:,sz(3)+1:sz(3)*2,:)   = -im_c1 - im_c2 + im_c3 + im_c4;
            Y(:,:,sz(3)*2+1:sz(3)*3,:) = -im_c1 + im_c2 - im_c3 + im_c4;
            Y(:,:,sz(3)*3+1:end,:)     = im_c1 - im_c2 - im_c3 + im_c4;
else
    sz = size(dzdy);
    if size(X, 4) == 1
        sz(4) = 1;
    end

            dzdy = dzdy/2;
            Y = zeros([sz(1)*2 sz(2)*2 sz(3)/4 sz(4)],'like',dzdy);
            Y(1:2:end, 1:2:end, : , :) = dzdy(:,:,1:sz(3)/4,:) - dzdy(:,:,sz(3)/4+1:sz(3)/2,:) - ...
                                         dzdy(:,:,sz(3)/2+1:3*sz(3)/4,:) + dzdy(:,:,sz(3)/4*3+1:end,:);
            Y(1:2:end, 2:2:end, : , :) = dzdy(:,:,1:sz(3)/4,:) - dzdy(:,:,sz(3)/4+1:sz(3)/2,:) + ...
                                         dzdy(:,:,sz(3)/2+1:3*sz(3)/4,:) - dzdy(:,:,sz(3)/4*3+1:end,:);
            Y(2:2:end, 1:2:end, : , :) = dzdy(:,:,1:sz(3)/4,:) + dzdy(:,:,sz(3)/4+1:sz(3)/2,:) - ...
                                         dzdy(:,:,sz(3)/2+1:3*sz(3)/4,:) - dzdy(:,:,sz(3)/4*3+1:end,:);
            Y(2:2:end, 2:2:end, : , :) = dzdy(:,:,1:sz(3)/4,:) + dzdy(:,:,sz(3)/4+1:sz(3)/2,:) + ...
                                         dzdy(:,:,sz(3)/2+1:3*sz(3)/4,:) + dzdy(:,:,sz(3)/4*3+1:end,:);

end

repo: MWCNN


text:
function [fx,fy,ft] = computeDerivatives2(im1,im2)
%
% function [fx,fy,fz] = computeDerivatives2(im1,im2)
%
% im1 and im2 are images
%
% [fx,fy,ft] are images, derivatives of the input images

filter = [0.03504 0.24878 0.43234 0.24878 0.03504];
dfilter = [0.10689 0.28461 0.0  -0.28461  -0.10689];

dx1 = conv2sep(im1,dfilter,filter,'valid');
dy1 = conv2sep(im1,filter,dfilter,'valid');
blur1 = conv2sep(im1,filter,filter,'valid');
dx2 = conv2sep(im2,dfilter,filter,'valid');
dy2 = conv2sep(im2,filter,dfilter,'valid');
blur2 = conv2sep(im2,filter,filter,'valid');

fx=(dx1+dx2)/2;
fy=(dy1+dy2)/2;
ft=(blur2-blur1);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [ Y ] = ktnormalize( Y,nrm,nnproj )
if nargin==1
    nrm=[];nnproj=[];
elseif nargin==2
    nnproj=[];
end
if isempty(nrm)
    nrm=2;
end
if isempty(nnproj)
    nnproj=true;
end

if strcmpi(class(Y),'ktensor')
    I=size(Y);
    N=numel(I);
    J=length(Y.lambda);
    switch nrm
        case 0
            for n=1:N
                ts=0:I(n):(J-1)*I(n);
                [temp idx]=max(abs(Y.U{n}));
                d=Y.U{n}(idx+ts);
                d=sign(d).*max(abs(d),eps);
                Y.U{n}=bsxfun(@rdivide,Y.U{n},d);
                Y.lambda=Y.lambda.*d';
            end
        case 1
            for n=1:N
                d=max(sum(abs(Y.U{n})),eps);
                Y.U{n}=bsxfun(@rdivide,Y.U{n},d);
                Y.lambda=Y.lambda.*d';
            end
        case 2
            for n=1:N
                d=max(sum(Y.U{n}.^2),eps);
                Y.U{n}=bsxfun(@rdivide,Y.U{n},d);
                Y.lambda=Y.lambda.*d';
            end
        otherwise
            disp('Unsupported norm.');
    end

    [temp idx]=sort(abs(Y.lambda),'descend');
    Y.lambda=Y.lambda(idx);
    for n=1:N
        Y.U{n}=Y.U{n}(:,idx);
    end

    %% nonnegative projection
    if nnproj
        ts=0:I(n):(J-1)*I(n);
        [temp idx]=max(abs(Y.U{n}));
        d=sign(Y.U{n}(idx+ts));
        Y.U{n}=bsxfun(@rdivide,Y.U{n},d);
        Y.lambda=Y.lambda.*d';
    end
end

end

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function imdb = setup_imdb_imagenet(dataDir, varargin)
% CNN_IMAGENET_SETUP_DATA  Initialize ImageNet ILSVRC CLS-LOC challenge data
%    This function creates an IMDB structure pointing to a local copy
%    of the ILSVRC CLS-LOC data.
%
%    The ILSVRC data ships in several TAR archives that can be
%    obtained from the ILSVRC challenge website. You will need:
%
%    ILSVRC2012_img_train.tar
%    ILSVRC2012_img_val.tar
%    ILSVRC2012_img_test.tar
%    ILSVRC2012_devkit.tar
%
%    Note that images in the CLS-LOC challenge are the same for the 2012, 2013,
%    and 2014 edition of ILSVRC, but that the development kit
%    is different. However, all devkit versions should work.
%
%    In order to use the ILSVRC data with these scripts, please
%    unpack it as follows. Create a root folder <DATA>, by default
%
%    data/ILSVRC2012
%
%    (note that this can be a simlink). Use the 'dataDir' option to
%    specify a different path.
%
%    Within this folder, create the following hierarchy:
%
%    <DATA>/images/train/ : content of ILSVRC2012_img_train.tar
%    <DATA>/images/val/ : content of ILSVRC2012_img_val.tar
%    <DATA>/images/test/ : content of ILSVRC2012_img_test.tar
%    <DATA>/ILSVRC2012_devkit : content of ILSVRC2012_devkit.tar
%
%    In order to speedup training and testing, it may be a good idea
%    to preprocess the images to have a fixed size (e.g. 256 pixels
%    high) and/or to store the images in RAM disk (provided that
%    sufficient RAM is available). Reading images off disk with a
%    sufficient speed is crucial for fast training.
opts.datasetName = 'ILSVRC2012';
opts.seed = 1 ;
opts = vl_argparse(opts, varargin) ;

if ischar(dataDir) && ~isempty(dataDir), opts.dataDir = dataDir;
else opts.dataDir = fullfile('data', opts.datasetName) ; end
opts.lite = false ;
opts = vl_argparse(opts, varargin) ;

% -------------------------------------------------------------------------
%                                                  Load categories metadata
% -------------------------------------------------------------------------

d = dir(fullfile(opts.dataDir, '*devkit*')) ;
if numel(d) == 0
  error('Make sure that ILSVRC data is correctly installed in %s', ...
    opts.dataDir) ;
end
devkitPath = fullfile(opts.dataDir, d(1).name) ;

% find metadata
mt = dir(fullfile(devkitPath, 'data', 'meta_clsloc.mat')) ;
if numel(mt) == 0
  mt = dir(fullfile(devkitPath, 'data', 'meta.mat')) ;
end
metaPath = fullfile(devkitPath, 'data', mt(1).name) ;

% find validation images labels
tmp = dir(fullfile(devkitPath, 'data', '*_validation_ground_truth*')) ;
valLabelsPath = fullfile(devkitPath, 'data', tmp(1).name) ;

% find validation images blacklist
tmp = dir(fullfile(devkitPath, 'data', '*_validation_blacklist*')) ;
if numel(tmp) > 0
  valBlacklistPath = fullfile(devkitPath, 'data', tmp(1).name) ;
else
  valBlacklistPath = [] ;
  warning('Could not find validation images blacklist file');
end

fprintf('using devkit %s\n', devkitPath) ;
fprintf('using metadata %s\n', metaPath) ;
fprintf('using validation labels %s\n', valLabelsPath) ;
fprintf('using validation blacklist %s\n', valBlacklistPath) ;

meta = load(metaPath) ;
cats = {meta.synsets(1:1000).WNID} ;
descrs = {meta.synsets(1:1000).words} ;

imdb.classes.name = cats ;
imdb.classes.description = descrs ;
imdb.imageDir = fullfile(opts.dataDir, 'images') ;

% -------------------------------------------------------------------------
%                                                           Training images
% -------------------------------------------------------------------------

fprintf('searching training images ...\n') ;
names = {} ;
labels = {} ;
for d = dir(fullfile(opts.dataDir, 'images', 'train', 'n*'))'
  [~,lab] = ismember(d.name, cats) ;
  ims = dir(fullfile(opts.dataDir, 'images', 'train', d.name, '*.JPEG')) ;
  names{end+1} = strcat([d.name, filesep], {ims.name}) ;
  labels{end+1} = ones(1, numel(ims)) * lab ;
  fprintf('.') ;
  if mod(numel(names), 50) == 0, fprintf('\n') ; end
  %fprintf('found %s with %d images\n', d.name, numel(ims)) ;
end
names = horzcat(names{:}) ;
labels = horzcat(labels{:}) ;

if numel(names) ~= 1281167
  warning('Found %d training images instead of 1,281,167. Dropping training set.', numel(names)) ;
  names = {} ;
  labels =[] ;
end

names = strcat(['train' filesep], names) ;

imdb.images.id = 1:numel(names) ;
imdb.images.name = names ;
imdb.images.set = ones(1, numel(names)) ;
imdb.images.label = labels ;

% -------------------------------------------------------------------------
%                                                         Validation images
% -------------------------------------------------------------------------

ims = dir(fullfile(opts.dataDir, 'images', 'val', '*.JPEG')) ;
names = sort({ims.name}) ;
labels = textread(valLabelsPath, '%d') ;

if numel(ims) ~= 50e3
  warning('Found %d instead of 50,000 validation images. Dropping validation set.', numel(ims))
  names = {} ;
  labels =[] ;
else
  if ~isempty(valBlacklistPath)
    black = textread(valBlacklistPath, '%d') ;
    fprintf('blacklisting %d validation images\n', numel(black)) ;
    keep = setdiff(1:numel(names), black) ;
    names = names(keep) ;
    labels = labels(keep) ;
  end
end

names = strcat(['val' filesep], names) ;

imdb.images.id = horzcat(imdb.images.id, (1:numel(names)) + 1e7 - 1) ;
imdb.images.name = horzcat(imdb.images.name, names) ;
imdb.images.set = horzcat(imdb.images.set, 2*ones(1,numel(names))) ;
imdb.images.label = horzcat(imdb.images.label, labels') ;

% -------------------------------------------------------------------------
%                                                               Test images
% -------------------------------------------------------------------------

ims = dir(fullfile(opts.dataDir, 'images', 'test', '*.JPEG')) ;
names = sort({ims.name}) ;
labels = zeros(1, numel(names)) ;

if numel(labels) ~= 100e3
  warning('Found %d instead of 100,000 test images', numel(labels))
end

names = strcat(['test' filesep], names) ;

imdb.images.id = horzcat(imdb.images.id, (1:numel(names)) + 2e7 - 1) ;
imdb.images.name = horzcat(imdb.images.name, names) ;
imdb.images.set = horzcat(imdb.images.set, 3*ones(1,numel(names))) ;
imdb.images.label = horzcat(imdb.images.label, labels) ;

% -------------------------------------------------------------------------
%                                                            Postprocessing
% -------------------------------------------------------------------------

% sort categories by WNID (to be compatible with other implementations)
[imdb.classes.name,perm] = sort(imdb.classes.name) ;
imdb.classes.description = imdb.classes.description(perm) ;
relabel(perm) = 1:numel(imdb.classes.name) ;
ok = imdb.images.label >  0 ;
imdb.images.label(ok) = relabel(imdb.images.label(ok)) ;

if opts.lite
  % pick a small number of images for the first 10 classes
  % this cannot be done for test as we do not have test labels
  clear keep ;
  for i=1:10
    sel = find(imdb.images.label == i) ;
    train = sel(imdb.images.set(sel) == 1) ;
    val = sel(imdb.images.set(sel) == 2) ;
    train = train(1:256) ;
    val = val(1:40) ;
    keep{i} = [train val] ;
  end
  test = find(imdb.images.set == 3) ;
  keep = sort(cat(2, keep{:}, test(1:1000))) ;
  imdb.images.id = imdb.images.id(keep) ;
  imdb.images.name = imdb.images.name(keep) ;
  imdb.images.set = imdb.images.set(keep) ;
  imdb.images.label = imdb.images.label(keep) ;
end

repo: ResNet-Matconvnet


text:
% For FDSA, this code provides values of a, A, and c in the gains a_k=a/(k+1+A)^.602 and
% c_k = c/(k+1)^.101. Code does not have built-in provision for gradient averaging (unlike
% SPSA) due to the relatively large no. of loss evals. required per gradient.
%
% Specify the dimension, loss function, and i.c.
%
global p sigma
p=10;
loss='loss4thnoise';
theta=ones(10,1);
alpha=0.602;
%
% User input on measurement noise level, expected no. of iterations in the SPSA run,
% desired magnitude of change in the theta elements, the number of SPSA gradient approximations
% that will be averaged, and the no. of loss evaluations to be used in the gain calculations
% here (note this no. should be divisible by twice the no. of averaged gradients).
%
step= input('What is the initial desired magnitude of change in the theta elements? ');
A = .10*input('What is the expected number of loss evaluations per run? ')...
   /(2*p);
sigma = input('What is the standard deviation of the measurement noise at i.c.? ');
c = max(sigma, .0001);
NL = input('How many loss function evaluations do you want to use in this gain calculation? ');
%
% Calculate the NL/(2*p) FD gradient estimates
%
rand('seed',31415927)
randn('seed',111113); %used in setting seed for noise in loss measurements
gbar=0;
e=eye(p);
atemp=zeros(p,1);
for i=1:NL/(2*p)
   ghat=0;
   for j=1:p
      thetaplus = theta + c*e(:,j);
      thetaminus = theta - c*e(:,j);
      yplus=feval(loss,thetaplus);
      yminus=feval(loss,thetaminus);
      ghat(j) = (yplus - yminus)/(2*c);
   end
   gbar=gbar+abs(ghat);
end
gbar=gbar/(NL/(2*p));
for i=1:p
   atemp(i)=step*((A+1)^alpha)/gbar(i);
end
atemp
a=min(atemp);
c
A
a


text:
function geometry_test2062 ( )

%*****************************************************************************80
%
%% TEST2062 tests TRIANGLE_AREA_HERON;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  dim_num = 3;

  t = [ ...
    1.0, 0.0, 0.0; ...
    0.0, 1.0, 0.0; ...
    0.0, 0.0, 1.0 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST2062\n' );
  fprintf ( 1, '  For a triangle in any dimension,\n' );
  fprintf ( 1, '  TRIANGLE_AREA_HERON computes the area;\n' );

  r8mat_transpose_print ( dim_num, 3, t, '  Triangle vertices:' );

  for j = 1 : 3

    s(j) = 0.0;

    jp1 = mod ( j, 3 ) + 1;

    for i = 1 : dim_num
      s(j) = s(j) + ( t(i,j) - t(i,jp1) ).^2;
    end

    s(j) = sqrt ( s(j) );

  end

  r8vec_print ( 3, s, '  Side lengths:' );

  area = triangle_area_heron ( s );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The area is %f\n', area );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [Rob, Frm] = rob2frm(Rob, Frm)

% FRM2ROB Robot to Frame information transfer
%   [Rob, Frm] = ROB2FRM(Rob, Frm) Updates the frame Frm from information
%   in Rob.
%
%   See also FRM2ROB.

%   Copyright 2015 Joan Sola @ IRI-CSIC-UPC

% Set 'used'
Frm.used = true;

% State
Frm.state.x = Rob.state.x;

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

%   # START GPL LICENSE

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

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

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function pdf = coupon_complete_pdf ( type_num, box_num )

%*****************************************************************************80
%
%% COUPON_COMPLETE_PDF evaluates the Complete Coupon Collection PDF.
%
%  Discussion:
%
%    PDF(TYPE_NUM;BOX_NUM) is the probability that, given an inexhaustible
%    supply of boxes, inside each of which there is one of TYPE_NUM distinct
%    coupons, which are uniformly distributed among the boxes, that it will
%    require opening exactly BOX_NUM boxes to achieve at least one of each
%    kind of coupon.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 August 2006
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Herbert Wilf,
%    Some New Aspects of the Coupon Collector's Problem,
%    SIAM Review,
%    Volume 48, Number 3, September 2006, pages 549-565.
%
%  Parameters:
%
%    Input, integer BOX_NUM, the number of boxes that had to be opened
%    in order to just get at least one of each coupon.
%    0 <= BOX_NUM.  If BOX_NUM < TYPE_NUM, then PDF is surely 0.
%
%    Input, integer TYPE_NUM, the number of distinct coupons.
%    1 <= TYPE_NUM.
%
%    Output, real PDF, the value of the PDF.
%

%
%  Nonsense cases.
%
  if ( box_num < 0 )

    pdf = 0.0;

  elseif ( type_num < 1 )

    pdf = 0.0;
%
%  Degenerate but meaningful case.
%
  elseif ( type_num == 1 )

    if ( box_num == 1 )
      pdf = 1.0;
    else
      pdf = 0.0;
    end
%
%  Easy cases.
%
  elseif ( box_num < type_num )

    pdf = 0.0;
%
%  General case.
%
  else

    factor = 1.0;
    for i = 1 : type_num
      factor = factor * i / type_num;
    end
    for i = type_num+1 : box_num
      factor = factor / type_num;
    end

    pdf = factor * stirling2_value ( box_num-1, type_num-1 );

  end

  return
end


text:
%  Numerical solution of polynomial systems by hom4ps Homotopy method
%
%  Syntax:  >> [S, var] = psolve( P )
%
%  Input:  P --- (cell array) the polynomial system to be solved
%  For example:
%     >> P = {'-x^5+y^5-3*y-1','5*y^4-3','-20*x+y-z'}
%
%  Output: S --- (matrix) numerical solutions (as columns of S)
%        var --- (cell array) the array of variables
%   For example, output
%      S =
%         -0.8264 + 0.6004i  -0.6092 - 1.0165i
%         -0.8801 + 0.0000i  -0.0000 - 0.8801i
%         15.6482 -12.0086i  12.1831 +19.4496i
%
%       var = {'x','y','z'}
%           means there are two numerical solutions
%       (x,y,z) =
%       (-0.8264+0.6004i, -0.8801+0.0000i, 15.6482-12.0086i)
%       (-0.6092-1.0165i, -0.8801+0.0000i, 12.1831+19.4496i)

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function varargout = interp_ungridded(pos_from, pos_to, varargin)

% INTERP_UNGRIDDED computes an interpolation matrix for two clouds of 3-D points
%
% To get the interpolated data, use as
%   [valto] = interp_ungridded(pos_from, pos_to, 'data', valfrom, ...)
% or to get the interpolation matrix itself, use as
%   [interpmat, distmat] = interp_ungridded(pos_from, pos_to, ...)
% where
%   pos_from  Nx3 matrix with the vertex positions
%   pos_to    Mx3 matrix with the vertex positions onto which the data should be interpolated
%
% Optional arguments are specified in key-value pairs and can be
%    data         = NxK matrix with functional data
%    distmat      = NxM matrix with precomputed distances
%    projmethod   = 'nearest', 'sphere_avg', 'sphere_weighteddistance', 'smudge'
%    triout       = triangulation for the second set of vertices
%    sphereradius = scalar
%    power        = scalar, power parameter as in the Inverse Distance Weighting function proposed by Shepard (default = 1).

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

if nargin<3
  ft_error('Not enough input arguments.');
end

% get the optional arguments
projmethod   = ft_getopt(varargin, 'projmethod');     % required
sphereradius = ft_getopt(varargin, 'sphereradius');   % required for some projection methods
powerparam   = ft_getopt(varargin, 'power', 1);
distmat      = ft_getopt(varargin, 'distmat');        % will be computed if needed and not present
triout       = ft_getopt(varargin, 'triout');
dat          = ft_getopt(varargin, 'data');           % functional data defined at pos_from
inside       = ft_getopt(varargin, 'inside');

hasdat    = ~isempty(dat);
hasinside = ~isempty(inside);
npos_from = size(pos_from, 1);
npos_to   = size(pos_to, 1);
dimres    = sqrt(sum((pos_to(2,:)-pos_to(1,:)).^2,2));

if hasinside
  % convert to boolean vector
  tmp         = false(npos_from,1);
  tmp(inside) = true;
  inside      = tmp;
  clear tmp;
else
  inside      = true(npos_from,1);
end

if isempty(distmat)
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % compute a distance matrix
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  switch projmethod
    case 'nearest'
      if ~isempty(sphereradius)
        ft_warning('sphereradius is not used for projmethod ''nearest''');
      end

      % determine the nearest voxel for each surface point
      ind     = find_nearest(pos_to, pos_from, 5);
      sel     = ind>0;
      indx    = 1:npos_to;
      distmat = sparse(indx(sel), ind(sel), ones(size(ind(sel))), npos_to, npos_from);

    case {'sphere_avg', 'sphere_weighteddistance'}
      if isempty(sphereradius)
        ft_error('sphereradius should be specified');
      end

      % compute the distance between voxels and each surface point
      dpos_fromsq = sum(pos_from.^2,2); % squared distance to origin
      dpos_tosq   = sum(pos_to.^2,2); % squared distance to origin
      maxnpnt = double(npos_to*ceil(4/3*pi*(sphereradius/max(dimres))^3)); % initial estimate of nonzero entries
      maxnpnt = min(maxnpnt, npos_to*npos_from);
      val   = nan(maxnpnt, 1);
      indx1 = nan(maxnpnt, 1);
      indx2 = nan(maxnpnt, 1);
      cnt = 1;
      for j = 1:npos_to

        % d   = dpos_tosq(j) + dpos_fromsq - 2 * pos_from * pos_to(j,:)';
        % sel = find(d<sphereradius.^2);
        % the following lines are equivalent to the previous 2 but uses fewer flops
        d    = sqrt(dpos_fromsq - 2 * pos_from * pos_to(j,:)' + dpos_tosq(j));
        sel  = find(d < sphereradius & inside);

        nsel = numel(sel);
        if nsel>0
          indx1(cnt:(cnt+nsel-1)) = j(ones(nsel,1));
          indx2(cnt:(cnt+nsel-1)) = sel(:);
          val(cnt:(cnt+nsel-1))   = d(sel) + eps('double');
          cnt = cnt + nsel;
        end
      end
      indx1(isnan(indx1)) = [];
      indx2(isnan(indx2)) = [];
      val(isnan(val))     = [];
      distmat = sparse(indx1, indx2, val, npos_to, npos_from);

    case 'smudge'
      if isempty(triout)
        ft_error('the ''smudge'' method needs a triangle definition');
      end

      [datin, loc] = ismember(pos_to, pos_from, 'rows');  % note that small numerical errors can cause them to be different
      [datout, S1] = smudge(datin, triout, 6);            % FIXME 6 is number of iterations, improve here

      sel = find(datin);
      S2  = sparse(sel(:), loc(datin), ones(npos_from,1), npos_to, npos_from);
      distmat = S1 * S2;

    otherwise
      ft_error('unsupported projection method');
  end % case projmethod
end % if isempty distmat

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% do something with the distance matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch projmethod
  case 'nearest'
    projmat         = distmat;

  case 'sphere_avg'
    projmat         = distmat;
    [ind1, ind2]    = find(projmat);
    nnz             = full(sum(spones(projmat),2));
    for k = 1:length(ind1)
      projmat(ind1(k),ind2(k)) = 1./nnz(ind1(k));
    end

  case 'sphere_weighteddistance'
    projmat         = distmat;
    [ind1, ind2, d] = find(projmat);
    projmat         = sparse(ind1, ind2, d.^-powerparam, npos_to, npos_from);
    [ind1, ind2, d] = find(projmat);
    normnz          = full(sum(projmat, 2));
    projmat         = sparse(ind1, ind2, d./normnz(ind1), npos_to, npos_from);

  case 'smudge'
    projmat = distmat;

  otherwise
    ft_error('unsupported projection method');
end  % case projmethod

if hasdat
  % return the interpolated values
  varargout{1} = projmat * dat;
else
  % return the interpolation and the distance matrix
  varargout{1} = projmat;
  varargout{2} = distmat;
end


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:  nickabattista@gmail.com
% 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 (nickabattista@gmail.com) know.
%
%--------------------------------------------------------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: updates the target point positions
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function targets = update_Target_Point_Positions(dt,current_time,targets)

%IDs = targets(:,1);                 % Stores Lag-Pt IDs in col vector
%xPts= targets(:,2);                 % Original x-Values of x-Target Pts.
%yPts= targets(:,3);                 % Original y-Values of y-Target Pts.
%kStiffs = targets(:,4);             % Stores Target Stiffnesses
%N_target = length(targets(:,1));    % Gives total number of target pts!

tR1 = 0.25;                          % Period A->B (and B->A) for RIGHT
tR2 = 0.25;                          % Period B->C (abd C->B) for RIGHT
periodR = (tR1+tR2);                % Time it takes to move to the right

off_Left = 0;%0.10*periodR;             % OFFSET FOR LEFT ARM
tL1 = 0.25;                          % Period A->B (and B->A) for LEFT
tL2 = 0.25;                          % Period B->C (abd C->B) for LEFT
periodL = (tL1+tL2);        % Time it takes to move to the LEFT

tR = rem(current_time,periodR);      % "Time" (moded out) according to right arm

tL = rem(current_time,periodL);      % "Time" (moded out) according to left arm

% Cubic Interpolation Information
p1 = 0.125;
p2 = 0.875;

% a COEFFICIENTS
a0 = 0;
a1 = 0;
a2 = 0;
a3 = 9.1429;

% b COEFFICIENTS
b0 = 0.0238;
b1 = -0.5714;
b2 = 4.5714;
b3 = -3.0476;

% c COEFFICIENTS
c0 = -8.1429;
c1 = 27.4286;
c2 = -27.4286;
c3 = 9.1429;

%
% JELLYFISH STATES:
%   XY(:,1) = 1st X-positions
%   XY(2:150,2) = 2nd X-positions (RIGHT)
%   XY(151:end,2) = 2nd X-positions (LEFT)
%
%   XY(:,3) = 1st Y-positions
%   XY(2:150,4) = 2nd Y-positions (RIGHT)
%   XY(151:end,4) = 2nd Y-positions (LEFT)
%
XY = read_In_Phase_Positions('XY_2Pos.txt');

% Coming from give_Me_Jelly:
NArm = 174;

% 1st and 2nd position of RIGHT arm
xR1 = XY(2:NArm+1,1); yR1 = XY(2:NArm+1,3);
xR2 = XY(2:NArm+1,2); yR2 = XY(2:NArm+1,4);

% 1st and 2nd position of LEFT arm
xL1 = XY(NArm+2:end,1); yL1 = XY(NArm+2:end,3);
xL2 = XY(NArm+2:end,2); yL2 = XY(NArm+2:end,4);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% START THE INTERPOLATING BETWEEN STATES!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%
%
% RIGHT ARM!
%
%%%%%%%%%%%%%

if tR <= tR1+off_Left % STATE A -> STATE B

    % Scaling time for appropriate use in interp. function so tRTilde\in[0,1]
    tRTilde = ( tR / tR1 );

    % Evaluate Piecewise Cubic Interpolation Poly
    if tRTilde<=p1
        gFUNC = a0 + a1*tRTilde + a2*tRTilde^2 + a3*tRTilde^3;
    elseif tRTilde<=p2
        gFUNC = b0 + b1*tRTilde + b2*tRTilde^2 + b3*tRTilde^3;
    else
        gFUNC = c0 + c1*tRTilde + c2*tRTilde^2 + c3*tRTilde^3;
    end

    % MOVING X-POSITIONS
    targets(2:NArm+1,2) = xR1 + gFUNC*( xR2 - xR1 );

    % MOVING Y-POSITIONS
    targets(2:NArm+1,3) = yR1 + gFUNC*( yR2 - yR1 );

elseif tR <= (tR1+tR2) % STATE B -> A

    % Scaling time for appropriate use in interp. function so tRTilde\in[0,1]
    tRTilde = (tR-tR1)/( tR2 );

    % Evaluate Piecewise Cubic Interpolation Poly
    if tRTilde<=p1
        gFUNC = a0 + a1*tRTilde + a2*tRTilde^2 + a3*tRTilde^3;
    elseif tRTilde<=p2
        gFUNC = b0 + b1*tRTilde + b2*tRTilde^2 + b3*tRTilde^3;
    else
        gFUNC = c0 + c1*tRTilde + c2*tRTilde^2 + c3*tRTilde^3;
    end

    % MOVING X-POSITIONS
    targets(2:NArm+1,2) = xR2 + gFUNC*( xR1 - xR2 );

    % MOVING Y-POSITIONS
    targets(2:NArm+1,3) = yR2 + gFUNC*( yR1 - yR2 );

end

%%%%%%%%%%%%%
%
% LEFT ARM!
%
%%%%%%%%%%%%%

if current_time >= off_Left

    if tL <= tL1 % STATE A -> STATE B

        % Scaling time for appropriate use in interp. function so tRTilde\in[0,1]
        tLTilde = (tL/tL1);

        % Evaluate Piecewise Cubic Interpolation Poly
        if tLTilde<=p1
            gFUNC = a0 + a1*tLTilde + a2*tLTilde^2 + a3*tLTilde^3;
        elseif tRTilde<=p2
            gFUNC = b0 + b1*tLTilde + b2*tLTilde^2 + b3*tLTilde^3;
        else
            gFUNC = c0 + c1*tLTilde + c2*tLTilde^2 + c3*tLTilde^3;
        end

        % MOVING X-POSITIONS
        targets(NArm+2:end,2) = xL1 + gFUNC*( xL2 - xL1 );

        % MOVING Y-POSITIONS
        targets(NArm+2:end,3) = yL1 + gFUNC*( yL2 - yL1 );

    elseif tL <= (tL1+tL2) % STATE B -> A

        % Scaling time for appropriate use in interp. function so tRTilde\in[0,1]
        tLTilde = (tL-tL1)/(tL2);

        % Evaluate Piecewise Cubic Interpolation Poly
        if tLTilde<=p1
            gFUNC = a0 + a1*tLTilde + a2*tLTilde^2 + a3*tLTilde^3;
        elseif tLTilde<=p2
            gFUNC = b0 + b1*tLTilde + b2*tLTilde^2 + b3*tLTilde^3;
        else
            gFUNC = c0 + c1*tLTilde + c2*tLTilde^2 + c3*tLTilde^3;
        end

        % MOVING X-POSITIONS
        targets(NArm+2:end,2) = xL2 + gFUNC*( xL1 - xL2 );

        % MOVING Y-POSITIONS
        targets(NArm+2:end,3) = yL2 + gFUNC*( yL1 - yL2 );

    end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: Reads in the # of vertex pts and all the vertex pts from the
%           .vertex file.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function PTS = read_In_Phase_Positions(struct_name)

filename = struct_name;  %Name of file to read in
fileID = fopen(filename);

% Read in the file, use 'CollectOutput' to gather all similar data together
% and 'CommentStyle' to to end and be able to skip lines in file.
C = textscan(fileID,'%f %f %f %f','CollectOutput',1);

fclose(fileID);     %Close the data file.

vertices = C{1};    %Stores all read in data in vertices (N+1,2) array

PTS = vertices(1:end,1:4);

author: nickabattista


text:
% a tutorial adapted from the fipy diffusion 1D example
% see: http://www.ctcms.nist.gov/fipy/examples/diffusion/index.html

% Written by Ali A. Eftekhari
% Last checked: June 2021
clc

%% define the domain
L = 5;  % domain length
Nx = 100; % number of cells
meshstruct = createMesh1D(Nx, L);
BC = createBC(meshstruct); % all Neumann boundary condition structure
BC.left.a(:) = 0; BC.left.b(:)=1; BC.left.c(:)=1; % left boundary
BC.right.a(:) = 0; BC.right.b(:)=1; BC.right.c(:)=0; % right boundary
x = meshstruct.cellcenters.x;
%% define the transfer coeffs
D_val = 1;
alfa = createCellVariable(meshstruct, 1);
Dave = createFaceVariable(meshstruct, D_val);
%% define initial values
c_old = createCellVariable(meshstruct, 0, BC); % initial values
c = c_old;
%% loop
dt = 0.001; % time step
final_t = 0.5;
for t=dt:dt:final_t
    % step 1: calculate divergence term
    RHS = divergenceTerm(Dave.*gradientTerm(c_old));
    % step 2: calculate the new value for internal cells
    c_old_int=internalCells(c_old);
    c_val_int = dt*excludeGhostRHS(meshstruct, RHS+constantSourceTerm(c_old))+c_old_int(:);
    c_analytical = 1-erf(x/(2*sqrt(D_val*t)));
    figure(1);plot(x, c_val_int, x, c_analytical, 'r--');drawnow;
    % Step 3: calculate the ghost cell values
    c.value=cellBoundary(c_val_int, BC);
    c_old = c;
    disp(t)
end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function [rifftype,msg,msgID] = readfourcc(fid)
%READFOURCC read four character code
%   [RIFFTYPE,MSG] = READFOURCC(FID) returns a four character code RIFFTYPE
%   from the file represented by FID. If the desired amount of data was not
%   read, then MSG is a string containing an error message, otherwise MSG is
%   empty.

%   Copyright 1984-2007 The MathWorks, Inc.
%   $Revision: 1.1.6.6 $  $Date: 2011/07/20 00:00:23 $

msg = '';
msgID='';
[rifftype, count] = fread(fid,4,'uchar');
rifftype = [char(rifftype)]';
if (count ~= 4)
  msgID='MATLAB:audiovideo:fourcc:invalidFourCC';
  msg = getString(message(msgID));
end
return;

author: kristinbranson

repo: JAABA


text:
function h = times(f,g)
%.*   Pointwise multiplication for DISKFUN objects.
%   F.*G multiplies DISKFUN objects F and G. Alternatively F or G could be
%   a double.

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

if ( isa(f, 'diskfun') && isa(g, 'diskfun') )
    % Grady's faster times for rank 1 functions:
    if ( length( f ) == 1 )
        [C, D, R] = cdr(f);
        h = g;
        onesForC = sqrt(abs(D)) * ones(1, length(g));
        onesForR = sign(D) * onesForC;
        h.cols = (C * onesForC) .* g.cols;
        h.rows = (R * onesForR) .* g.rows;
        % Switch the parity terms if needed:
        % plus*plus = plus -> do nothing
        % plus*minus = minus -> do nothing
        % minus*minus = plus -> switch minus to plus
        % minus*plus = minus -> switch plus to minus
        if ( isempty(f.idxPlus) )  % f is a minus
            h.idxMinus = g.idxPlus;
            h.idxPlus = g.idxMinus;
        end
        % Both terms have to be non-zero at the poles for the product to be
        % non-zero at the poles.
        h.nonZeroPoles = f.nonZeroPoles & g.nonZeroPoles;
    elseif ( length( g ) == 1 )
         h = times(g, f);
    else
        % Let separableApprox deal with the multiplication
        h = times@separableApprox(f, g);
    end
else
    % Let separableApprox deal with the multiplication
    h = times@separableApprox(f, g);
end

end


text:
function [ar,e,dc]=lpccovar(s,p,t,w)
%LPCCOVAR performs covariance LPC analysis [AR,E,DC]=(S,P,T)
%
%  Inputs:  S(NS)    is the input signal
%           P        is the order (default: 12)
%           T(NF,:)  specifies the frames size details: each row specifies one frame
%                    T can be a cell array if rows have unequal numbers of values
%                       T(:,1) gives the start of the analysis interval: must be >P
%                       T(:,2) gives the end of the anaylsis interval [default: t(:+1,1)-1]
%                       subsequent pairs can be used to specify multiple disjoint segments
%                    If T is omitted, T(1,1)=P+1, T(1,2)=NS;
%                    The elements of t need not be integers.
%           W(NS)    The error at each sample is weighted by W^2 (default: 1)
%
% Outputs:  AR(NF,P+1)  are the AR coefficients with AR(:,1) = 1
%           E(NF,2)     is the energy in the residual and in the original window.
%                       sqrt(E) is often called the 'gain' of the LPC filter.
%           DC          is the DC component of the signal S. If this output is included,
%                       the LPC equations are modified to include a DC offset.

% Notes:
%
% (1) For speech processing P should be at least 2*F*L/C where F is the sampling
%     frequency, L the vocal tract length and C the speed of sound. For a typical
%     male (l=17 cm) this gives f/1000.
%
% (2) Each analysis frame should contain at least 2P samples. If note (1) is followed
%     this implies at least 2 ms of speech signal per frame.
%
% (3) It can be advantageous to restrict the analysis regions to time intervals
%     when the glottis is closed (closed-phase analysis). This can be achieved by
%     setting the T input parameter appropriately. If the closed-phase is shorter than
%     2 ms then two or more successive closed-phases should be used by defining 4 or more
%     elements in the corresponding row of T.
%
% (4) A previous version of this routine allowed T() to have a single row which would
%     be replicated for the entire file length. This has be removed because it gave rise
%     to an ambiguity.

%  Bugs: sould really detect a singular matrix and reduce the order accordingly

%	   Copyright (C) Mike Brookes 1995
%      Version: $Id: lpccovar.m 713 2011-10-16 14:45:43Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 2 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 can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

s = s(:); % make it a column vector
if nargin < 2 p=12; end;
if nargin < 3 t=[p+1 length(s)]; end;
wq = nargin>3;
[nf,ng]=size(t);
if iscell(t)
   t{nf+1}=length(s)+1;
else
   if rem(ng,2)
      t(:,end+1)=[t(2:nf,1)-1; length(s)];
   end
end
ar=zeros(nf,p+1);
ar(:,1)=1;
e=zeros(nf,2);
dc=zeros(nf,1);
d0=nargout >2;

rs=(1:p);

for jf=1:nf
   if iscell(t)
      tj=t{jf};
      if rem(length(tj),2)
         tj(end+1)=t{jf+1}(1)-1;
      end
   else
      tj=t(jf,:);
   end

   ta = ceil(tj(1));
   tb = floor(tj(2));
   cs = (ta:tb).';
   for js=3:2:length(tj)
      ta = ceil(tj(js));
      tb = floor(tj(js+1));
      cs = [cs; (ta:tb).'];
   end
   %disp(cs([logical(1); (cs(2:end-1)~=cs(1:end-2)+1)|(cs(2:end-1)~=cs(3:end)-1); logical(1)])');
   nc = length(cs);
   pp=min(p,nc-d0);
   dm=zeros(nc,pp);	% predefine shape
   dm(:) = s(cs(:,ones(1,pp))-rs(ones(nc,1),1:pp));
   if nargout>2
      if wq
         dm = [ones(nc,1) dm].*w(cs(:,ones(1,1+pp)));
         sc=(s(cs).*w(cs));
         aa = (dm\sc).';
      else
         dm = [ones(nc,1) dm];
         sc=s(cs);
         aa = (dm\sc).';
      end
      ar(jf,2:pp+1) = -aa(2:pp+1);
      e(jf,1)=sc.'*(sc - dm*aa.');
      e(jf,2)=sc.'*sc;
      dc(jf)=aa(1)/sum(ar(jf,:));
   else
      if wq
         dm = dm.*w(cs(:,ones(1,pp)));
         sc=(s(cs).*w(cs));
         aa = (dm\sc).';
      else
         sc=s(cs);
         aa = (dm\sc).';
      end;
      ar(jf,2:pp+1) = -aa;
      if nargout>1
         e(jf,1)=sc.'*(sc - dm*aa.');
         e(jf,2)=sc.'*sc;
      end
   end
end
if ~nargout
   lpcar2db(ar,127);
end

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function varargout = Sim_Single_Voxel_Curve_GUI(varargin)
% SIM_SINGLE_VOXEL_CURVE_GUI MATLAB code for Sim_Single_Voxel_Curve_GUI.fig
%      SIM_SINGLE_VOXEL_CURVE_GUI, by itself, creates a new SIM_SINGLE_VOXEL_CURVE_GUI or raises the existing
%      singleton*.
%
%      H = SIM_SINGLE_VOXEL_CURVE_GUI returns the handle to a new SIM_SINGLE_VOXEL_CURVE_GUI or the handle to
%      the existing singleton*.
%
%      SIM_SINGLE_VOXEL_CURVE_GUI('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in SIM_SINGLE_VOXEL_CURVE_GUI.M with the given input arguments.
%
%      SIM_SINGLE_VOXEL_CURVE_GUI('Property','Value',...) creates a new SIM_SINGLE_VOXEL_CURVE_GUI or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before Sim_Single_Voxel_Curve_GUI_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to Sim_Single_Voxel_Curve_GUI_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help Sim_Single_Voxel_Curve_GUI

% Last Modified by GUIDE v2.5 10-Jul-2017 21:03:41

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @Sim_Single_Voxel_Curve_GUI_OpeningFcn, ...
                   'gui_OutputFcn',  @Sim_Single_Voxel_Curve_GUI_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before Sim_Single_Voxel_Curve_GUI is made visible.
function Sim_Single_Voxel_Curve_GUI_OpeningFcn(hObject, eventdata, handles, varargin)
set(findobj('Name','qMRLab'),'pointer', 'watch'); drawnow;

Model = varargin{1};
setappdata(0,'Model',Model);

if ~isfield(handles,'opened')
    handles.output = hObject;
    handles.Model = varargin{1};
    % OPTIONS
    if isprop(handles.Model,'Sim_Single_Voxel_Curve_buttons')
        opts = handles.Model.Sim_Single_Voxel_Curve_buttons;
    else
        opts = {'SNR',50};
    end

    % Generate Buttons
    handles.options = GenerateButtonsWithPanels(opts,handles.OptionsPanel);

     % Create CALLBACK for buttons
    ff = fieldnames(handles.options);
    for ii=1:length(ff)
        switch get(handles.options.(ff{ii}),'Style')
%             case 'popupmenu'
%                 val =  find(cell2mat(cellfun(@(x) strcmp(x,Model.options.(ff{ii})),get(handles.OptionsPanel_handle.(ff{ii}),'String'),'UniformOutput',0)));
%                 set(handles.OptionsPanel_handle.(ff{ii}),'Value',val);
%             case 'checkbox'
%                 set(handles.OptionsPanel_handle.(ff{ii}),'Value',Model.options.(ff{ii}));
%             case 'edit'
%                 set(handles.OptionsPanel_handle.(ff{ii}),'String',Model.options.(ff{ii}));
            case 'togglebutton'
                set(handles.options.(ff{ii}),'Callback',@(src,event) ModelSimOptions_Callback(handles));
        end
    end

    Nparam = length(handles.Model.xnames);
    FitOptTable(:,1)=handles.Model.xnames(:);
    if isprop(handles.Model,'st') && ~isempty(handles.Model.st)
        FitOptTable(:,2)=mat2cell(handles.Model.st(:),ones(Nparam,1));
    elseif isprop(handles.Model,'lb') && ~isempty(handles.Model.lb) && isprop(handles.Model,'ub') && ~isempty(handles.Model.ub)
        FitOptTable(:,2) = mat2cell(mean([handles.Model.lb(:), handles.Model.ub(:)],2),ones(Nparam,1));
    end
    set(handles.ParamTable,'Data',FitOptTable)

%     % launch plot
%     UpdatePlot_Callback(hObject, eventdata, handles)

    % opened
    handles.opened = true;
end
% Update handles structure
guidata(hObject, handles);
set(findobj('Name','qMRLab'),'pointer', 'arrow'); drawnow;

% --- Executes on button press in UpdatePlot.
function UpdatePlot_Callback(hObject, eventdata, handles)
Model_new = getappdata(0,'Model');
if ~isempty(Model_new) && strcmp(class(Model_new),class(handles.Model))
    handles.Model = Model_new;
end
set(findobj('Name','Single Voxel Curve'),'pointer', 'watch'); drawnow;
if isgraphics(handles.SimCurveAxe)
    axes(handles.SimCurveAxe)
end

xtable = get(handles.ParamTable,'Data');
x=cell2mat(xtable(~cellfun(@isempty,xtable(:,2)),2))';

FitResults = Sim_Single_Voxel_Curve(handles.Model,x,button_handle2opts(handles.options));
hold off;

% put results in table
ff = fieldnames(FitResults);
for ii=1:length(ff)
    index = strcmp(xtable(:,1),ff{ii});
    if find(index)
        xtable{index,3} = FitResults.(ff{ii})(1);
        xtable{index,4} = round((FitResults.(ff{ii})(1) - xtable{index,2})/xtable{index,2}*100);
    else
        xtable{end+1,1} = ff{ii};
        xtable{end,3} = FitResults.(ff{ii})(1);
    end
end

% CRLB
SNR = str2double(get(handles.options.SNR,'String'));
[~,~,~,F] = SimCRLB(handles.Model,handles.Model.Prot.(handles.Model.MRIinputs{1}).Mat,x,1/SNR);

for ii=1:sum(~handles.Model.fx)
    ll=find(~handles.Model.fx);
    xtable{ll(ii),5}=F(ii)*100;
end
set(handles.ParamTable,'Data',xtable);
set(findobj('Name','Single Voxel Curve'),'pointer', 'arrow'); drawnow;

% --- Executes on button press in Options panel.
function ModelSimOptions_Callback(handles)
xtable = get(handles.ParamTable,'Data');
x=cell2mat(xtable(~cellfun(@isempty,xtable(:,2)),2))';
xnew = SimOpt(handles.Model,x,button_handle2opts(handles.options));
if ~isempty(xnew) % update the ParamTable in the GUI
    Nparam = length(handles.Model.xnames);
    xtable(1:Nparam,2) = mat2cell(xnew',ones(Nparam,1));
    set(handles.ParamTable,'Data',xtable);
end

% --- Outputs from this function are returned to the command line.
function varargout = Sim_Single_Voxel_Curve_GUI_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function a = dif1 ( m, n )

%*****************************************************************************80
%
%% DIF1 returns the DIF1 matrix.
%
%  Discussion:
%
%    For a set of N points X(I) with equal spacing H, and a set of data
%    values Y(I) associated with those points, the product
%    1/(2*H) * A * Y returns an approximation to the first derivative
%    of Y(X) at the interior points X(2:N-1).
%
%  Example:
%
%    N = 5
%
%    0 +1  .  .  .
%   -1  0 +1  .  .
%    . -1  0 +1  .
%    .  . -1  0 +1
%    .  .  . -1  0
%
%  Rectangular Properties:
%
%    A is banded, with bandwidth 3.
%
%    A is tridiagonal.
%
%    Because A is tridiagonal, it has property A (bipartite).
%
%    A is integral: int ( A ) = A.
%
%    A is Toeplitz: constant along diagonals.
%
%  Square Properties:
%
%    A is antisymmetric: A' = -A.
%
%    Because A is antisymmetric, it is normal.
%
%    Because A is normal, it is diagonalizable.
%
%    If N is even, then A is nonsingular.
%    If N is odd, then A is singular.
%
%    If N is even, det ( A ) = 1.0.
%    If N is odd, det ( A ) = 0.0.
%
%    If N is odd, a null vector is ( 1, 0, 1, 0, ..., 1, 0, 1 )..
%
%    The family of matrices is nested as a function of N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns of A.
%
%    Output, real A(M,N), the matrix.
%
  a = zeros ( m, n );

  for i = 1 : m

    if ( 0 < i - 1 )
      a(i,i-1) = -1.0;
    end

    if ( i + 1 <= n )
      a(i,i+1) = +1.0;
    end

  end

  return
end


text:
function Z = xkron(X,Y)
%XKRON Kronecker tensor product.
%
%   XKRON(X,Y), where X and Y are matrices, is the Kronecker tensor product
%   of X and Y.  The result is a large matrix formed by taking all possible
%   products between the elements of X and those of Y.
%
%   For example, if X is 2-by-3, and Y is any N-dimensional array, then
%   KRON(X,Y) is
%
%     [ X(1,1)*Y  X(1,2)*Y  X(1,3)*Y
%       X(2,1)*Y  X(2,2)*Y  X(2,3)*Y ]
%
%   If either X or Y is sparse, and neither X nor Y has dimension
%   greater than 2, then only nonzero elements are multiplied in the
%   computation, and the result is sparse.
%
%   XKRON works also when X and/or Y are N-dimensional arrays.  There does
%   not seem to be a unique definition of an N-dimensional Kronecker tensor
%   product, so I simply picked the one I found most logical.
%
%   If X is 3-by-4-by-2, and Y is any N-dimensional array, then
%
%     cat( 3, [ X(1,1,1)*Y  X(1,2,1)*Y  X(1,3,1)*Y  X(1,4,1)*Y ;  ...
%               X(2,1,1)*Y  X(2,2,1)*Y  X(2,3,1)*Y  X(2,4,1)*Y ;  ...
%               X(3,1,1)*Y  X(3,2,1)*Y  X(3,3,1)*Y  X(3,4,1)*Y ], ...
%             [ X(1,1,2)*Y  X(1,2,2)*Y  X(1,3,2)*Y  X(1,4,2)*Y ;  ...
%               X(2,1,2)*Y  X(2,2,2)*Y  X(2,3,2)*Y  X(2,4,2)*Y ;  ...
%               X(3,1,2)*Y  X(3,2,2)*Y  X(3,3,2)*Y  X(3,4,2)*Y ] )

%   Author:      Peter J. Acklam
%   Time-stamp:  2002-03-03 13:50:29 +0100
%   E-mail:      pjacklam@online.no
%   URL:         http://home.online.no/~pjacklam

   xs = size(X);
   ys = size(Y);
   xd = length(xs);
   yd = length(ys);

   if ( issparse(X) | issparse(Y) ) & ( xd <= 2 ) & ( yd <= 2 )

      % When at least one input is sparse, and neither argument has
      % dimension greater than 2, the result is sparse.

      mx = xs(1); nx = xs(2);
      my = ys(1); ny = ys(2);

      [ix, jx, sx] = find(X);
      [iy, jy, sy] = find(Y);
      ix = ix(:); jx = jx(:); sx = sx(:);
      iy = iy(:); jy = jy(:); sy = sy(:);
      kx = ones(size(sx));
      ky = ones(size(sy));
      t = my*(ix-1)';
      ik = t(ky,:)+iy(:,kx);
      t = ny*(jx-1)';
      jk = t(ky,:)+jy(:,kx);
      Z = sparse(ik, jk, sy*sx.', mx*my, nx*ny);

   else

      dz = max(xd, yd);
      xs = [ xs ones(1,dz-xd) ];
      ys = [ ys ones(1,dz-yd) ];

      v = reshape( reshape( 1:2*dz, dz, 2 )', 1, 2*dz );
      Z = reshape( Y(:)*X(:).', [ ys xs ] );
      Z = reshape( permute( Z, v ), xs.*ys );

   end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function [y,deriv] = linTrans_adaptive(w,map,transmap)
% This is an MV2DF. See MV2DF_API_DEFINITION.readme.
%
% Applies linear transform y = map(w). It needs the transpose of map,
% transmap for computing the gradient. map and transmap are function
% handles.

if nargin==0
    test_this();
    return;
end

if isempty(w)
    y = @(w)linTrans_adaptive(w,map,transmap);
    return;
end

if isa(w,'function_handle')
    outer = linTrans_adaptive([],map,transmap);
    y = compose_mv(outer,w,[]);
    return;
end

y = map(w);
y = y(:);

deriv = @(g2) deriv_this(g2,map,transmap,numel(w));
end

function [g,hess,linear] = deriv_this(g2,map,transmap,wlen)
g = transmap(g2,wlen);
g = g(:);
%linear = false;  % use this to test linearity of map, if in doubt
linear = true;
hess = @(d) hess_this(map,d);
end

function [h,Jd] = hess_this(map,d)
h = [];
if nargout>1
    Jd = map(d);
    Jd = Jd(:);
end

end

function test_this()
first = 2;
len = 3;
map = @(w) w(first:first+len-1);
function w = transmap_test(y,sz)
    w=zeros(sz,1);
    w(first:first+len-1)=y;
end
transmap = @(y,sz) transmap_test(y,sz);
f = linTrans_adaptive([],map,transmap);
test_MV2DF(f,randn(5,1));
end

author: nesl