text:
function [connectivity, pValues, freq] = bst_granger_spectral(X, Y, Fs, order, inputs)
% BST_GRANGER_SPECTRAL  Granger causality at each frequency between any two
%                       signals.
%
% Inputs:
%   X                 - first set of signals, one signal per row
%                       [X: A x N or A x N x T matrix]
%   Y                 - second set of signals, one signal per row
%                       [Y: B x N or B x N x T matrix]
%                       (default: Y = X)
%   Fs                - sampling rate (we assume uniform sampling rate)
%                       [FS: scalar, FS > freq(end)*2]
%   order             - maximum order of autogressive model
%                       [p: integer > 1, default 10]
%   inputs            - structure of parameters:
%   |-freq            - frequencies of interest if desired
%   |-freqResolution  - maximum freq resolution in Hz, to limit NFFT
%   |                   [DF: double, default [] (i.e. no limit)]
%   |-nTrials         - # of trials in concantenated signal
%   |-flagFPE         - if true, optimize order for autoregression
%   |                   if false (default), use same order in autoregression
%   |-standardize     - if true (default), remove mean from each signal
%   |                   if false, assume signal has already been detrended
%
% Outputs:
%   connectivity      - A x B matrix of spectral Granger causalities from
%                       source to sink. For each signal pair (a,b) we calculate
%
%                                        S_{sink} (f)
%                  ------------------------------------------------------------
%                  S_{sink}(f) - |H_{sink, source} (f)|^2 sigma_{source | sink}
%
%                       with S_{sink}(f) as the power spectral density of a @ f
%                            H_{sink, source}(f) as the transfer function @ f
%                            sigma_{source | sink} as the conditional variance
%                             of the residual at b given the residual at a,
%                             calculated using the residual covariance.
%                       By default, GC(a,a) = 0 if Y is empty.
%                       [C: MX x MY x NF matrix]
%   pValues           - parametric p-value for corresponding spectral Granger
%                       causality in mean estimate
%                       [P: MX x MY x NF matrix]
%   freq              - frequencies corresponding to the previous two metrics
%                       [F: NF x 1 vector]
%
% See also BST_GRANGER, BST_COHERENCE_MVAR.
%
% Call:
%   connectivity = bst_granger_spectral(X, Y, 200, 10, inputs); % general call
%   connectivity = bst_granger_spectral(X, [], 200, 30, inputs); % every pair
%   connectivity = bst_granger_spectral(X, [], 200, 30, inputs); % more variance
% Parameter examples:
%   inputs.freq           = 0:0.1:100; % specify desired frequencies
%   inputs.freqResolution = 0.1; % have a high-point FFT
%   inputs.nTrials        = 9; % use trial-average covariances in AR estimation
%   inputs.flagFPE        = true; % use AR model with best information criteria

% Note: for those following Chicharro2012, I have used the equivalence
%                       sigma_{xx}^(xy) |H_{xx}^(xy) (w)|^2
%                                       =
%                S_{xx}(w) - sigma_{yy}^(xy) |H_{xy}^(xy) (w)|^2
% which follows Geweke1982 instead. As Chicharro notes, it may be better to use his
% formulation because it separates out instantaneous causality I think.

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

% reformat to a 2D matrix if necessary, and pull out # of trials
if ndims(X) == 3
  inputs.nTrials = size(X,3);
  X = reshape(X, size(X, 1), []);
elseif ~isfield(inputs, 'nTrials')
  inputs.nTrials = 1;
end

% lengths of things
nSamples = size(X, 2);
nTimes = nSamples / inputs.nTrials;

% standardization: zero mean & unit variance, remove linear & quadratic trends as well
if isfield(inputs, 'standardize') && inputs.standardize
  detrender = [ ...
    ones(1, nTimes); ... % constant trend in data
    linspace(-1, 1, nTimes); ... % linear trend in data
    3/2 * linspace(-1, 1, nTimes).^2 - 1/2 ... % quadratic trend in data
  ];

  % detrend X
  for iTrial = 1:inputs.nTrials
    X(:, (iTrial-1)*nTimes + (1:nTimes)) = X(:, (iTrial-1)*nTimes + (1:nTimes)) - ( X(:, (iTrial-1)*nTimes + (1:nTimes)) / detrender ) * detrender;
    X(:, (iTrial-1)*nTimes + (1:nTimes)) = diag( sqrt(sum(X(:, (iTrial-1)*nTimes + (1:nTimes)).^2, 2)) ) \ X(:, (iTrial-1)*nTimes + (1:nTimes));
  end

  % detrend Y only if it is not empty
  if ~isempty(Y)
    Y = reshape(Y, size(Y, 1), []); % reshape to 2D matrix first
    for iTrial = 1:inputs.nTrials
      Y(:, (iTrial-1)*nTimes + (1:nTimes)) = Y(:, (iTrial-1)*nTimes + (1:nTimes)) - ( Y(:, (iTrial-1)*nTimes + (1:nTimes)) / detrender ) * detrender;
      Y(:, (iTrial-1)*nTimes + (1:nTimes)) = diag( sqrt(sum(Y(:, (iTrial-1)*nTimes + (1:nTimes)).^2, 2)) ) \ Y(:, (iTrial-1)*nTimes + (1:nTimes));
    end
  end
end

% number of FFT bins required
if ~isfield(inputs, 'freq') || isempty(inputs.freq) % frequencies of interest are not defined
  if isfield(inputs, 'freqResolution') && ~isempty(inputs.freqResolution) && (size(X,2) > round(Fs / inputs.freqResolution))
    nFFT = 2^nextpow2( round(Fs / inputs.freqResolution) );
  else % use a default frequency resolution of 1Hz to mirror the standard frequency resolution in bst_coherence_welch.m
    nFFT = 2^nextpow2( round(Fs / 1) );
  end
elseif numel(inputs.freq) == 1
  nFFT = inputs.freq;
else
  nFFT = 2^nextpow2(max( length(inputs.freq)-1, (Fs/2) / min(diff( sort(inputs.freq(:)) )) )) * 2;
end

% default: Order 10 used in BrainStorm
if isempty(order)
  order = 10;
end

% default: Single-trial
if ~isfield(inputs, 'nTrials') || isempty(inputs.nTrials)
  inputs.nTrials = 1;
end

% default: do not optimize model order
if ~isfield(inputs, 'flagFPE') || isempty(inputs.flagFPE)
  inputs.flagFPE = false;
end

%% Differentiate between auto-causality and cross-causality for speed

if isempty(Y) % auto-causality between signals in X, so we can halve the number of models to estimate

  % preallocate the spectral causality matrix
  connectivity = zeros(size(X, 1), size(X, 1), nFFT/2);

  % iterate over all pairs of sinks & sources
  for iX = 1:size(X, 1)
    for iY = iX+1 : size(X, 1) % to avoid auto-causality

        % two-variate model for given source
        [transfers, noiseCovariance, order] = bst_mvar([X(iX, :); X(iY, :)], order, inputs.nTrials, inputs.flagFPE);

        % spectra and power of forward system
        [spectra, freq, forward] = bst_granger_spectral_spectrum(transfers, noiseCovariance, nFFT, Fs);

        % Geweke-Granger spectral causality from source to sink
        unrestricted = abs(spectra(1, 1, :)); restriction = forward(1, 2, :); % S_{sink}(f) and |H_{sink, source} (f)|^2, w/ abs to get rid of 1e-16 imag part
        residualVariance = noiseCovariance(2,2) - noiseCovariance(2, 1) / noiseCovariance(1, 1) * noiseCovariance(1, 2); % partial variance of source
        restricted = abs(unrestricted) - abs(restriction).^2 * residualVariance; % S_{sink} (f) - |H_{sink, source} (f)|^2 sigma_{source | sink}
        connectivity(iX, iY, abs(restricted) > 1e-60) = unrestricted(abs(restricted) > 1e-60) ./ restricted(abs(restricted) > 1e-60) - 1;% Geweke-Granger
        % sigma_{source | sink} = partial covariance which is the formula above (sigma_{source} - rho_{source, sink} / sigma_{sink} * rho_{sink, source})

        % Geweke-Granger spectral causality from sink back to source (to halve the # of MVAR fittings)
        unrestricted = abs(spectra(2, 2, :)); restriction = forward(2, 1, :); % S_{source}(f) and |H_{source, sink} (f)|^2, w/ abs to get rid of 1e-16 imag part
        residualVariance = noiseCovariance(1,1) - noiseCovariance(1, 2) / noiseCovariance(2, 2) * noiseCovariance(2, 1); % partial variance of sink
        restricted = abs(unrestricted) - abs(restriction).^2 * residualVariance; % S_{source} (f) - |H_{source, sink} (f)|^2 sigma_{sink | source}
        connectivity(iY, iX, abs(restricted) > 1e-60) = unrestricted(abs(restricted) > 1e-60) ./ restricted(abs(restricted) > 1e-60) - 1; % Geweke-Granger
        % sigma_{sink | source} = partial covariance which is the formula above (sigma_{sink} - rho_{sink, source} / sigma_{source} * rho_{source, sink})

    end

    % diagonal will equal the maximum of all inflows and outflows for iX, specific to each frequency
    connectivity(iX, iX, :) = max( max(connectivity(iX, :, :), [], 2), max(connectivity(:, iX, :), [], 1) );

  end

else % we have to use all pairs of signals

  % preallocate the spectral causality matrix
  connectivity = zeros(size(X, 1), size(Y, 1), nFFT/2);
  duplicates = zeros(0, 2);

  % iterate over all pairs of sinks & sources
  for iX = 1:size(X, 1)
    for iY = 1:size(Y, 1)

      if max(abs(X(iX, :) - Y(iY, :))) > eps % by default, if X(sink) = Y(source), the causality is 0 everywhere

        % 2-variate model for given source
        [transfers, noiseCovariance, order] = bst_mvar([X(iX, :); Y(iY, :)], order, inputs.nTrials, inputs.flagFPE);

        % spectra and power of forward system
        [spectra, freq, forward] = bst_granger_spectral_spectrum(transfers, noiseCovariance, nFFT, Fs);
        spectra = abs(spectra(1, 1, :)); % limit to the autospectrum of the sink S_{sink} (f) and take absolute value to rid the 1e-16 imaginary part
        forward = forward(1, 2, :); % limit to the cross-transfer from source to sink H_{sink, source} (f)

        % partial covariance of residual
        residualVariance = noiseCovariance(2,2) - noiseCovariance(2,1) / noiseCovariance(1, 1) * noiseCovariance(1, 2);
        % sigma_{source | sink} = partial covariance which is the formula above (sigma_{source} - rho_{source, sink} / sigma_{sink} * rho_{sink, source})

        % Geweke-Granger spectral causality from source to sink
        restricted = spectra - abs(forward).^2 * residualVariance; % S_{sink} (f) - |H_{sink, source} (f)|^2 sigma_{source | sink}
        connectivity(iX, iY, abs(restricted) > 1e-60) = spectra(abs(restricted) > 1e-60) ./ restricted(abs(restricted) > 1e-60) - 1; % Geweke-Granger

      else % save duplicates to modify later

        duplicates(end+1, :) = [iX iY]; %#ok<AGROW>

      end

    end
  end

  % for duplicate indices, set the causality value to the maximum of all inflows for iX and outflows for iY
  for iDuplicate = 1:size(duplicates, 1)
    connectivity(duplicates(iDuplicate, 1), duplicates(iDuplicate, 2), :) = max( ...
      max(connectivity(duplicates(iDuplicate, 1), :, :), [], 2), ...
      max(connectivity(:, duplicates(iDuplicate, 2), :), [], 1) ...
      );
  end

end

%% Interpolate to desired frequencies & perform statistics if desired
if isfield(inputs, 'freq') && ~isempty(inputs.freq) && numel(inputs.freq) > 1
  connectivity = permute(interp1(freq, permute(connectivity, [3 1 2]), inputs.freq), [2 3 1]);
  % pValues = permute(interp1(freq, permute(pValues, [3 1 2]), inputs.freq), [2 3 1]);
  freq = inputs.freq;
end

pValues = NaN; % no parametric p-values for now

end

%% ======================================================== estimation for multivariate autoregression ========================================================
function [spectra, freq, forward] = bst_granger_spectral_spectrum(transfers, noiseCovariance, nFFT, Fs)
% BST_MVAR_SPECTRUM     Calculate the parametric spectra of a bivariate system
%                       with given MVAR coefficients & an estimate of the
%                       covariance matrix of the innovation (i.e. noise)
%                       process.
%
% Inputs:
%   transfers         - transfer matrices in AR process
%                       [A: 2 x 2P matrix, P = order]
%   noiseCovariance   - variance of residuals
%                       [C: 2 x 2 matrix]
%   nFFT              - number of FFT bins to calculate spectra
%                       [NF: positive number, usually power of 2]
%   Fs                - sampling frequency of data
%                       [FS: double, FS > freq(end)*2]
%
% Outputs:
%   spectra           - cross-spectrum between each pair of variables
%                       [S: 2 x 2 x NF/2 matrix]
%   freq              - frequencies used based on # of FFT bins
%                       [F: NF/2 x 1 matrix]
%   forward           - forward transform in frequency from source to sink
%                       [H: 2 x 2 x NF/2 matrix]
%   --> all outputs have length NF/2, ignoring frequenices past Fs/2
%
% Call:
%   spectra = bst_mvar_spectrum(transfers, C, 512, 200); % basic
%   spectra = bst_mvar_spectrum(transfers, C, 2048, 200); % add FFT interp
%   [spectra, freq] = bst_mvar_spectrum(transfers, C, 512, 200); % grab freqs
%   [~, ~, forward] = bst_mvar_spectrum(transfers, C, 64, 200); % causality

% Notes:
% The DTFT we want is
% H = I - sum_p C_p e^{-j 2 pi f/Fs p} = sum_p D_p e^{-j 2 pi f/Fs (p-1)} where D_1 = 1 and D_p = -C_{p-1} for p > 1
% MatLab's FFT provides
% G = sum_n x_n e^{-j 2 pi (k-1)/N (n-1)}
% so the matching is p = n and (k-1)/N = f/Fs, after vectorizing H.

% frequencies to estimate cross-spectra
freq = Fs/2*linspace(0, 1, nFFT/2 + 1);
freq(end) = [];

%% Inverse of transfer function in frequency
% Inverse transfer means the transfer function from the sources to the innovations. We calculate this at each frequency.
% The inverse transfer from source a to innovation b is
%                1 - \sum_{n=1}^N a_{ab} [n] e^{-j2pi * f * n}

inverse = fft(reshape([eye(2) -transfers], 4, [])', nFFT); % reshape so we can do a vector FFT quickly
inverse = inverse(1:nFFT/2, :); % restrict to the first symmetric half of the spectrum

%% Forward transfer in autoregressive model

% pieces of 2x2 inverse
forward = zeros(2, 2, nFFT/2); % an important caveat is that I did not reshape inverse earlier;
forward(1,1,:) = inverse(:,4); forward(1,2,:) = -inverse(:,3); % as a result, we have a column-wise index mapping: 4 = (2,2) and 3 = (1,2)
forward(2,1,:) = -inverse(:,2); forward(2,2,:) = inverse(:,1); % in addition, 2 = (2,1) and 1 = (1,1). then these elements fit the 2x2 matrix inverse
detInverse = inverse(:,1).*inverse(:,4) - inverse(:,3).*inverse(:,2); % the same thing happens here, using the indexing to avoid a reshape() call

% normalization by determinant to get inverse
forward = bst_bsxfun(@rdivide, forward, reshape(detInverse, [1 1 length(freq)])); % complete the inversion by dividing by frequency-dependent determinant

%% Cross-spectrum from forward model
% The forward transfer is the inverse of the inverse transfer at each frequency f. Denoted H(f), the power spectral density is then HH' at each frequency.

% the loop that we won't use
% for idxFreq = 1:nFFT
%   spectra(:, :, idxFreq) = H(:,:,idxFreq) * noiseCovariance * H(:,:,idxFreq)';
% end

% formula for the matrix multiplication
spectra(1,1,:) = ... % 1,1 element is H_11 C_11 H_11^* + H_12 C_12 H_11^* + H_11 C_12 H_12^* + H_12 C_22 H_12^* (and I combine the middle two into 2 Re{.})
  noiseCovariance(1,1) * forward(1,1,:) .* conj(forward(1,1,:)) ...
  + noiseCovariance(1,2) * real(forward(1,2,:) .* conj(forward(1,1,:))) * 2 ...
  + noiseCovariance(2,2) * forward(1,2,:) .* conj(forward(1,2,:));
spectra(1,2,:) = ... % 1,2 element is H_11 C_11 H_21^* + H_12 C_12 H_21^* + H_11 C_12 H_22^* + H_12 C_22 H_22^*
  noiseCovariance(1,1) * forward(2,1,:) .* conj(forward(1,1,:)) ...
  + noiseCovariance(1,2) * forward(2,2,:) .* conj(forward(1,1,:)) ...
  + noiseCovariance(1,2) * forward(2,1,:) .* conj(forward(1,2,:)) ...
  + noiseCovariance(2,2) * forward(1,2,:) .* conj(forward(1,2,:));
spectra(2,1,:) = conj(spectra(1,2,:)); % for speed, force the 2,1 element to be the conjugate of the 1,2 element so we have conjugate symmetry
spectra(2,2,:) = ... % 2,2 element is H_21 C_11 H_21^* + H_22 C_12 H_21^* + H_21 C_12 H_22^* + H_22 C_22 H_22^* (and I combine the middle two into 2 Re{.})
  noiseCovariance(1,1) * forward(2,1,:) .* conj(forward(2,1,:)) ...
  + noiseCovariance(1,2) * real(forward(2,2,:) .* conj(forward(2,1,:))) * 2 ...
  + noiseCovariance(2,2) * forward(2,2,:) .* conj(forward(2,2,:));

%% Normalize for sampling rate, and truncate if necessary

% normalization
forward = forward / sqrt(Fs);
spectra = spectra / Fs;

%% Confidence intervals
% taken from Kay, pp 194-195
% alpha = 1 - ci;
% original = -sqrt(2) * erfcinv( 2 * ( 1 - alpha/2) );
% lower = abs(spectra) * (1 - sqrt(2 * order / nTimes) * original);
% upper = abs(spectra) * (1 + sqrt(2 * order / nTimes) * original);

end %% <== FUNCTION END

author: brainstorm-tools


text:
% @author: Maziar Raissi

function params_list = Fractional()
% quantile(params_list,[0.025 0.25 0.50 0.75 0.975])

clc; close all;

plt = 1;
save_plt = 0;

addpath ..
addpath ../Utilities
addpath ../Kernels/Fractional
addpath ../Utilities/export_fig

function CleanupFun()
    rmpath ..
    rmpath ../Utilities
    rmpath ../Kernels/Fractional
    rmpath ../Utilities/export_fig
end

finishup = onCleanup(@() CleanupFun());

rng('default')

set(0,'defaulttextinterpreter','latex')

%% Load Data
[t_star, x_star, u_star] = gen_data(10^6, 0.01, 5, 100);
% u_star --> 300x5
% t_star --> 5x1
% x_star --> 300x1
N_star = size(x_star,1);
nsteps = size(t_star,1)-1;

if plt ==1
    figure(2);
    plot_surface(t_star, x_star, u_star, '$t$', '$x$', '$u(t,x)$');
    view(3);

    drawnow()
end

%% Setup
noise = 0.0;
u_data = (1 + noise*randn(size(u_star))).*u_star;

N0 = 100;
N1 = 100;
%% Optimize model
params_list = [];
hyp = [log([1.0 0.1]) 0.1 1.2 -4.0];
idx1 = randsample(N_star, N0);
step = 1;
for i = nsteps:step:nsteps
    dt = t_star(i+step) - t_star(i);

    idx0 = idx1;
    x0 = x_star(idx0,:);
    u0 = u_data(idx0,i);

    idx1 = randsample(N_star,N1);
    x1 = x_star(idx1,:);
    u1 = u_data(idx1,i+step);

    model = HPM(x1, u1, x0, u0, dt, hyp);
    model = model.train(500);

    hyp = model.hyp;
    params_list = [params_list; hyp(3:4)];

    [pred_n_star, var_n_star] = model.predict(x_star);
    var_n_star = abs(diag(var_n_star));

    error = norm(pred_n_star - u_star(:,i+step))/norm(u_star(:,i+step));

    fprintf(1,'=========================\n');
    fprintf(1,'Step: %d, Time = %.2f\n\nNLML = %.2f, Error = %.2e\n\n', i, ...
        t_star(i+step), model.NLML, error);

    str = sprintf('%.2f  ', params_list(end,:));
    fprintf('Parameters: %s\n\n', str)

    str = sprintf('%.2f  ', median(params_list,1));
    fprintf('Median: %s\n', str)
    fprintf(1,'=========================\n\n');

    if plt == 1
        if ~exist('fig','var')
            fig = figure(2);
        end
        set(fig,'units','normalized','outerposition',[0 0 1 1])
        clf

        subplot(3,1,1);
        tit = sprintf('Time: %.2f\n%d training points', t_star(i), N0);
        plot_data_1D(x_star, u_star(:,i), x0, u0, '$x$', '$u(t,x)$', tit);

        subplot(3,1,2);
        tit = sprintf('Time: %.2f\n%d training points', t_star(i+step), N1);
        plot_data_1D(x_star, u_star(:,i+step), x1, u1, '$x$', '$u(t,x)$', tit);

        subplot(3,1,3);
        plot_prediction_1D(x_star, u_star(:,i+step), pred_n_star, var_n_star, ...
            '$x$', '$u(t,x)$', tit);

        drawnow;
    end

end

if save_plt == 1
    export_fig ./Figures/Burgers.png -r300
end

end

function [t_star, x_star, u_star] = gen_data(N,dt,m,n)

    pos = cumsum(sqrt(dt)*randn(N,1));

    M = 0;

    P = zeros(length(pos)-m,m);
    for i = 1:length(pos)-m
        y = pos(i+1:i+m) - pos(i);
        M = max([M, max(abs(y))]);
        P(i,:) = y;
    end

    M = M/2;

    bins = linspace(-M,M,n+1);
    x_star = linspace(M*(1/n-1), M*(1-1/n), n)';
    t_star = dt*(1:1:m)';
    u_star = zeros(n,m);

    figure(1)
    clf
    hold
    for i = 1:m
        h = histogram(P(:,i), bins, 'Normalization', 'pdf');
        u_star(:,i) = h.Values;
    end

end

author: maziarraissi

repo: HPM

save_path: github-repos/MATLAB/maziarraissi-HPM


text:
classdef J0Kernel < Kernel
    % Object of type kernel but optimized for J0 so that the computation of the
    % radial quadrature goes faster
    % Radial quadrature is only made of 1 component !

    properties
        R=1, C=1;
        % Such that G(r) = C*besselj(0,R*r)
    end

    methods
        function[this] = J0Kernel(RR,CC)
            if nargin == 0
                RR = 1;
            end
            if nargin <= 1
                CC = 1;
            end
            this@Kernel(@(x)(CC*besselj(0,RR*x)),@(x)(-RR*CC*besselj(1,RR*x)))
            this.C = CC;
            this.R = RR;
        end
        function[out] = dilatation(this,lambda)
            out = J0Kernel(lambda*this.R,this.C);
            % No changes in gamma_est, scalFunc nor normFunc
        end
        function[out] = mtimes(this,mu)
            if isa(this,'Kernel')
                out = J0Kernel(this.R,mu*this.C);
            else
                out = mtimes(mu,this);
            end
        end
        function[out] = radialQuadKernel(this,a,tol,varargin)
            Cmem = this.C;
            this = (1/Cmem)*this;
            rq = RadialQuadrature(a,this,tol,varargin{:});
            out = Cmem*rq;
        end
    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function varargout=airProp(T, prop)
%---------------------------------------------
% Interpolates thermodynamic air properties
% Temp. range: 250 - 1200 K
% According to W. C. Reynolds, cp. Heywood,
% Internal Combustion Engine Fundamentals,
% p. 127 and p. 912.
%
% Values are in SI-units:
%
% 	col-#	prop. 	units
%	------------------------
%	1		T		K
%	2		h		kJ/kg
%	3		u		kJ/kg
%	4		Psi		kJ/(kgK)
%	5		Fi		kJ/(kgK)
%	6		pr		-
%	7		vr		-
%	8		cp		kJ/(kgK)
%	9		cv		kJ/(kgK)
%
% Example 1:	out=airProp(304, 'Fi')
% Example 2:
% 		[h,u]=airProp([333 999],{'h' 'u'})
%---------------------------------------------
% (c)2004 by Stefan Billig
%---------------------------------------------
% Last Change:	19-May-2004
%---------------------------------------------

% check # of input arguments
if ~isequal(nargin,2)
	error('airProp requires 2 input arguments!')
	return
% check temperature request
elseif find(T<250) | find(T>1200) | ~isnumeric(T)
	error('Valid temperature range: 250 >= T[K] >= 1200')
	return
end
% get table
load propTabAir
% if multi property request
if iscell(prop)
	% scan along cells
	for idx=1:length(prop)
		% identify property column
		col=find(strcmp(propInfo,prop(idx)));
		if isempty(col)
			disp(['Property "' char(prop(idx)) '" not recognized!'])
		else
			% create output
			varargout{idx}=interp1(airTab(:,1),airTab(:,col),T);
		end
	end
% single property request
else
	% identify property column
	col=find(strcmp(propInfo,prop));
	if isempty(col)
		disp(['Property "' prop '" not recognized!'])
	else
		% create output
		varargout{1}=interp1(airTab(:,1),airTab(:,col),T);
	end
end

author: Sable

repo: mcbench-benchmarks


text:
function [c,s] = wavefast(x,n,varargin)
%WAVEFAST Computes the FWT of a '3-D extended' 2-D array.
%   [C, L] = WAVEFAST(X, N, LP, HP, EXTMODE) computes 'PAGES' 2D N-level
%   FWTs of a 'ROWS x COLUMNS x PAGES' matrix X with respect to
%   decomposition filters LP and HP and boundary extension mode EXTMODE.
%
%   [C, L] = WAVEFAST(X, N, WNAME, EXTMODE) performs the same operation
%   but fetches filters LP and HP for wavelet WNAME using WAVEFILTER.
%
%   Scale parameter N must be less than or equal to log2 of the maximum
%   image dimension.  Filters LP and HP must be even. If EXTMODE = 'SYM'
%   (the default), X is symmetrically extended; if X = [c1 c2 c3 ... cn]
%   (in 1D), then its symmetric extension would be [... c3 c2 c1 c1 c2
%   c3 ... cn cn cn-1 cn-2 ...]. If EXTMODE = 'PER', X is periodically
%   extended to [... cn-1 cn c1 c2 c3 ... cn c1 c2 ...].
%
%   OUTPUTS:
%     Vector C is a coefficient decomposition vector:
%
%      C = [ a1(n)...ak(n) h1(n)...hk(n) v1(n)...vk(n)
%            d1(n)...dk(n) h1(n-1)... d1(1)...dk(1) ]
%
%     where ai, hi, vi, and di for i = 0,1,...k are columnwise vectors
%     containing approximation, horizontal, vertical, and diagonal
%     coefficient matrices, respectively, and k is the number of pages
%     in the 3-D extended array X. C has 3n + 1 sections where n is the
%     number of wavelet decompositions.
%
%     Matrix S is an [(n+2) x 2] bookkeeping matrix if k = 1; else it is
%     [(n+2) x 3]:
%
%      S = [ sa(n, :); sd(n, :); sd(n-1, :); ... ; sd(1, :); sx ]
%
%     where sa and sd are approximation and detail size entries.
%
%   See also WAVEBACK and WAVEFILTER.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

% Check the input arguments for reasonableness.
extmode = 'SYM';
if ischar(varargin{1})
   [lp,hp] = wavefilter(varargin{1},'d');
   if nargin > 3
      extmode = varargin{2};
   end
else
   lp = varargin{1};   hp = varargin{2};
   if nargin > 4
      extmode = varargin{3};
   end
end

% Get the filter length, 'lp', input array size, 'sx', and number of
% pages, 'pages', in extended 2-D array x.
fl = length(lp);       sx = size(x);        pages = size(x,3);

if ((~ismatrix(x)) && (ndims(x) ~= 3)) || (min(sx) < 2) ...
      || ~isreal(x) || ~isnumeric(x)
   error('X must be a real, numeric 2-D or 3-D matrix.');
end

if (~ismatrix(lp)) || ~isreal(lp) || ~isnumeric(lp) ...
      || (~ismatrix(hp)) || ~isreal(hp) || ~isnumeric(hp) ...
      || (fl ~= length(hp)) || rem(fl,2) ~= 0
   error(['LP and HP must be even and equal length real, ' ...
      'numeric filter vectors.']);
end

if ~isreal(n) || ~isnumeric(n) || (n < 1) || (n > log2(max(sx)))
   error(['N must be a real scalar between 1 and ' ...
      'log2(max(size((X))).']);
end

% Init the starting output data structures and initial approximation.
c = [];        s = sx(1:2);
app = cell(pages, 1);
for i = 1:pages
   app{i} = double(x(:,:,i));
end

% For each decomposition ...
for i = 1:n
   % Extend the approximation.
   [app,keep] = extend(app,fl,pages,extmode);

   % Convolve rows with HP and downsample. Then convolve columns
   % with HP and LP to get the diagonal and vertical coefficients.
   rows = convolve(app,hp,'row',fl,keep,pages,extmode);
   coefs = convolve(rows,hp,'col',fl,keep,pages,extmode);
   c = addcoefs(c,coefs,pages);
   s = [size(coefs{1}); s];
   coefs = convolve(rows,lp,'col',fl,keep,pages,extmode);
   c = addcoefs(c,coefs,pages);

   % Convolve rows with LP and downsample. Then convolve columns
   % with HP and LP to get the horizontal and next approximation
   % coeffcients.
   rows = convolve(app,lp,'row',fl,keep,pages,extmode);
   coefs = convolve(rows,hp,'col',fl,keep,pages,extmode);
   c = addcoefs(c,coefs,pages);
   app = convolve(rows,lp,'col',fl,keep,pages,extmode);
end

% Append the final approximation structures.
c = addcoefs(c,app,pages);
s = [size(app{1}); s];
if ~ismatrix(x)
   s(:,3) = size(x,3);
end

%----------------------------------------------------------------------%
function nc = addcoefs(c,x,pages)
% Add 'pages' array coefficients to the wavelet decomposition vector.

nc = c;
for i = pages:-1:1
   nc = [x{i}(:)' nc];
end

%----------------------------------------------------------------------%
function [y,keep] = extend(x,fl,pages,extmode)
% Extend the 'pages' arrays of x in both dimensions and return
% 'keep' to determine the number of coefficients to keep after
% convolution and downsampling.

y = cell(pages,1);
for i = 1:pages
   if strcmpi(extmode,'SYM')
      keep = floor((fl + size(x{i}) - 1) / 2);
      y{i} = padarray(x{i},[(fl - 1) (fl - 1)],'symmetric','both');
   elseif strcmpi(extmode,'PER')
      keep = size(x{i});
      y{i} = padarray(x{i},[fl/2 fl/2],'circular','both');
   else
      error('Invalid extension mode!');
   end
end

%----------------------------------------------------------------------%
function y = convolve(x,h,type,fl,keep,pages,extmode)
% For the 'pages' 2-D arrays in x, convolve the rows or columns with
% h, downsample, and extract the section defined by 'keep'.

y = cell(pages,1);
for i = 1:pages
   if strcmp(type,'row')
      if strcmpi(extmode, 'SYM')
         y{i} = conv2(x{i},h);
         y{i} = y{i}(:,1:2:end);
         y{i} = y{i}(:,fl / 2 + 1:fl / 2 + keep(2));
      else
         y{i} = conv2(x{i},h,'valid');
         y{i} = y{i}(:, 2:2:2 * ceil(keep(2) / 2));
      end
   else
      if strcmpi(extmode,'SYM')
         y{i} = conv2(x{i},h');
         y{i} = y{i}(1:2:end,:);
         y{i} = y{i}(fl / 2 + 1:fl / 2 + keep(1), :);
      else
         y{i} = conv2(x{i},h','valid');
         y{i} = y{i}(2:2:2 * ceil(keep(1) / 2),:);
      end
   end
end

author: dipum

repo: dipum-toolbox

save_path: github-repos/MATLAB/dipum-dipum-toolbox


text:
function [r] = issubfield(s, f)

% ISSUBFIELD tests for the presence of a field in a structure just like the standard
% Matlab ISFIELD function, except that you can also specify nested fields
% using a '.' in the fieldname. The nesting can be arbitrary deep.
%
% Use as
%   f = issubfield(s, 'fieldname')
% or as
%   f = issubfield(s, 'fieldname.subfieldname')
%
% This function returns true if the field is present and false if the field
% is not present.
%
% See also ISFIELD, GETSUBFIELD, SETSUBFIELD

% Copyright (C) 2005, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% 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: issubfield.m 2885 2011-02-16 09:41:58Z roboos $

try
  getsubfield(s, f);    % if this works, then the subfield must be present
  r = true;
catch
  r = false;                % apparently the subfield is not present
end

author: PatternRecognition

repo: OpenBMI


text:
function test_stepsize_alg_demo()
% demonstration file for original stepsize algorithm.
%
% This file illustrates how to set user's own stepsize algorithm in case of linear
% regression problem. This demonstrates SGD and SVRG algorithms.
%
% This file is part of SGDLibrary.
%
% Created by H.Kasai on Sep. 25, 2017

    clc;
    clear;
    close all;

    %% generate synthetic data
    % set number of dimensions
    d = 10;
    % set number of samples
    n = 1000;
    % generate data
    data = logistic_regression_data_generator(n, d);

    %% define problem definitions
    problem = logistic_regression(data.x_train, data.y_train, data.x_test, data.y_test);

%     %% define user-defined stepsize algorithm
%     function step = my_stepalg(iter, options)
%         step = options.step_init / (10 + iter*0.5);
%     end
%

    %% perform algorithms SGD and SVRG
    options.w_init = data.w_init;
    options.step_init = 0.01;
    options.verbose = 2;

    options.step_alg = 'fix';
    [w_sgd_fix, info_sgd_fix] = sgd(problem, options);

    options.step_alg = 'decay';
    [w_sgd_decay, info_sgd_decay] = sgd(problem, options);

    options.step_alg = 'decay-2';
    [w_sgd_decay2, info_sgd_decay2] = sgd(problem, options);

    options.stepsizefun = @my_stepalg;  % set my_stepalg (user-defined stepsize algorithm)
    [w_sgd_my, info_sgd_my] = sgd(problem, options);

    %% display cost/optimality gap vs number of gradient evaluations
    display_graph('grad_calc_count','cost', {'SGD (fix)','SGD (decay)', 'SGD (decay-2)', 'SGD (My stepsize algorithm)'}, ...
            {w_sgd_fix, w_sgd_decay w_sgd_decay2, w_sgd_my}, {info_sgd_fix, info_sgd_decay, info_sgd_decay2, info_sgd_my});

end

    %% define user-defined stepsize algorithm
    function step = my_stepalg(iter, options)
        step = options.step_init / (10 + iter*0.5);
    end

author: hiroyuki-kasai

repo: SGDLibrary


text:
function state = is_octave()
% is_octave: checks if the script is running on Octave or MATLAB.
%
% INPUT
%   None.
%
% OUTPUT
%	  state, TRUE or FALSE.
%
%   Copyright (C) 2014, Rodrigo Gonzalez, all rights reserved.
%
%   This file is part of NaveGo, an open-source MATLAB toolbox for
%   simulation of integrated navigation systems.
%
%   NaveGo is free software: you can redistribute it and/or modify
%   it under the terms of the GNU Lesser General Public License (LGPL)
%   version 3 as published by the Free Software Foundation.
%
%   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 Lesser General Public License for more details.
%
%   You should have received a copy of the GNU Lesser General Public
%   License along with this program. If not, see
%   <http://www.gnu.org/licenses/>.
%
% References:
%
%
% Version: 001
% Date:    2021/12/09
% Author:  Rodrigo Gonzalez <rodralez@frm.utn.edu.ar>
% URL:     https://github.com/rodralez/navego

  state = exist('OCTAVE_VERSION', 'builtin') ~= 0;
end

author: rodralez

repo: NaveGo

save_path: github-repos/MATLAB/rodralez-NaveGo


text:
function val=numMPartitions(n,m,getUpToM)
%%NUMMPARTITIONS Get the number of ways of partitioning the integer n into
%                m parts, where order does not matter (unlike compositions,
%                where order matters). That is the number of ways of
%                getting m nonzero integers that sum to n. The amount of
%                memory for this implementation is proportional to n*m and
%                thus the solution becomes slow for very large n and m
%                values.
%
%INPUTS: n The positive n>0 integer to be partitioned.
%        m The number of parts into which the integer is to be partitioned.
% getUpToM This parameter changes what is returned. If false (the default
%          if omitted or an empty matrix is passed), the number of
%          ways of partitioning n into m parts is returned. If true, then
%          the total number of ways of partitioning n into m or fewer parts
%          is returned.
%
%OUTPUTS: val The number of ways of partitioning the integer n into m
%             parts or if getUpToM=true, then this is the total number of
%             ways of partitioning n into m or fewer parts.
%
%The counting method is based on the recurrence relation in Chapter 7.2.1.4
%of [1]. The algorithm could be implemented to use less memory. However,
%that would require shifting rows in a table up (to reuse memory) and would
%make the code a lot more complicated.
%
%REFERENCES:
%[1] D. E. Knuth, The Art of Computer Programming. Vol. 4, Fascicle 3:
%    Generating all Combinations and Partitions, Upper Saddle River, NJ:
%    Addison-Wesley, 2009.
%
%October 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<3||isempty(getUpToM))
    getUpToM=false;
end

P=eye(n+1,m+1);
P(:,1+1)=1;

for curN=2:n
    for curK=1:min(curN,m)
        P(curN+1,curK+1)=P(curN-1+1,curK-1+1)+P(curN-curK+1,curK+1);
    end
end
if(getUpToM)
    val=sum(P(end,:));
else
    val=P(end,end);
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function chararray_out=textart(image_data, varargin)
%TEXTART - Conert digital image to ASCII text image
%   TEXTART(IMAGE) returns a character array based on IMAGE, which may be
%     either a filename of an image file, or the pixel data of an image
%     (2-D or 3-D)
%   TEXTART(...,PROPERTIES) allows the caller to specifiy PROPERTIES, pair-wise:
%     'zoom' - enlarges (or reduces) image (default 1.0)
%     'aspect' - Aspect correction ratio (default 1.0)
%     'gamma' - Skews the intensity levels (default 1.0)
%     'file' - Specifies (text) file to which to write character data (default none)
%     'bg' - Specifies background color of original image, typically 0. (default none)
%     'interp' - Interpolate per 1969 I.D.G. MacLead algorithm
%
%   Example:
%     chardata = textart('c:\matlab6p5\bin\win32\matlab.ico','zoom',8,'bg',0);

% 02/05/09 jaj Add 'interp' function per 1969 I.D.G. MacLead algorithm,
%   "Pictorial Output with a Line Printer", IEEE Transactions on Computers,
%   Volume 19 ,  Issue 2  (February 1970), Pages 160-162

% Default properties
properties = struct('zoom',1.0,'aspect',1.0,'gamma',1.0,'file',[],'bg',[],'interp',false);

if nargin == 0
  error('Missing image file name to convert')
end

properties = getopt(properties,varargin{:});

% Load/convert the image data
if ischar(image_data)
  image_data = imread(image_data);
end

image_data = double(image_data); % Some images are uint8, etc.

% If 3-D image (color) convert to black/white
if ndims(image_data) == 3
  image_data = mean(image_data,3);
end

if ~isempty(properties.bg)
  image_data(image_data == properties.bg) = 255;
end

% Process then input arguments

% Define the set of ASCII characters to use, and compute their denisties
% (intensities) by drawing them, capuring the images, and counting the
% number of pixels lit up.
% (Need do this only once.)
[density_table,h_char_pix,w_char_pix] = get_density_data;

A = 1-image_data/255; % Convert image to B/W, negate, normalize

% Make the size of a divisibale by the number of characters
h_im_pix = size(A,1); % height of image, in pixels
w_im_pix = size(A,2); % width of image, in pixels

h_im_char = round(properties.aspect*properties.zoom*h_im_pix/h_char_pix); % height of image, in characters
w_im_char = round(properties.zoom*w_im_pix/w_char_pix); % width of image, in characters

h_step_pix = h_im_pix/h_im_char;
w_step_pix = w_im_pix/w_im_char;

% Gamma adjustment
A = A.^properties.gamma;
% hist(A(:),100) % View inensity distribution

h_step_idx_list = round(h_step_pix*(1:h_im_char));
w_step_idx_list = round(w_step_pix*(1:w_im_char));

h_step_idx_list = max(h_step_idx_list,1);
w_step_idx_list = max(w_step_idx_list,1);

% This filtering technique is just a convenient means to average
% the intensity of pixels in sucessive blocks.
filt_len_h = floor(h_step_pix);
filt_len_w = floor(w_step_pix);

filtcoef_h = [ones(filt_len_h,1);h_step_pix-filt_len_h];
filtcoef_w = [ones(filt_len_w,1);w_step_pix-filt_len_w];

filtcoef_h = filtcoef_h ./ sum(filtcoef_h); % Normalize
filtcoef_w = filtcoef_w ./ sum(filtcoef_w); % Normalize

A1 = filter(filtcoef_h, 1, A, [], 2);
A1 = filter(filtcoef_w, 1, A1, [], 1);
A1 = A1(h_step_idx_list,w_step_idx_list); % decimate

% % This section provides (a very slow) preview
% for h=1:h_im_char
%   for w=1:w_im_char
%     c = 1-A1(h,w);
%     try
%       plot(w,h,'.','color',[c c c]); hold on
%     catch
%       fprintf('uh oh\n');
%     end
%   end
% end
% hold off;axis equal;axis tight;set(gca,'ydir','reverse');
% set(gca,'DataAspectRatio',[h_char_pix w_char_pix 1])

density = [density_table.density];
ch_list = {density_table.ch_list};
nchar   = [density_table.nchar];
% A1 is the density of each block in the picture. No compare against the
% character densities, using clever(?) vector processing.
d_im = repmat(A1,[1 1 length(ch_list)]); % Now a 3-d matrix
d_ch = repmat(reshape(density,[1 1 length(ch_list)]),[size(A1),1]);

% 03/05/09 jaj Use I.D.G. MacLeod algorithm per 1969 paper "Pictorial
% output with a line printer"

if ~properties.interp
  % Original implementation
  [dif,ind] = min(abs(d_im-d_ch),[],3); % Find closest density match: block to character
else
  dist = d_im-d_ch;
  [junk,ind] = min(abs(dist),[],3); % Find closest density match: block to character
  dif = A1-density(ind);
  ind0 = ind - (dif <= 0);
  ind1 = ind + (dif > 0);
  ind1(ind0 == 0) = ind1(ind0 == 0)+1;
  ind0(ind0 == 0) = ind0(ind0 == 0)+1;
  % Linear interpolation
  d0 = A1-density(ind0);
  d1 = density(ind1)-A1;
  p = d0./(d0+d1);
  r = rand(size(A1));
  sel = r <= p;
  ind(sel) = ind0(sel);
  ind(~sel) = ind1(~sel);
end
chararray = repmat(' ',size(A1));
for ii=1:numel(A1)
  ix = ind(ii);
  chararray(ii) = ch_list{ix}(ceil(nchar(ix)*rand(1)));
end
%chararray=ch_list(ind);

% Return the character data only under certain conditions
if nargout > 0 || isempty(properties.file)
  chararray_out = chararray;
end

% Write character data to file, if requested
if ~isempty(properties.file)
  fid=fopen(properties.file,'wt');
  if fid <= 0
    error(['Unable to open file ',properties.file]);
  end
  for y=1:size(chararray,1);
    fprintf(fid,'%s\n',chararray(y,:));
  end
  fclose(fid);
end

%---------------------------------------------------------------
function [density_table_out,h_char_pix_out,w_char_pix_out] = get_density_data

persistent h_char_pix w_char_pix density_table
fontname = 'courier'; % Must be fixed width font
if isempty(density_table)
  ch_list = char(32:126);
  figure; subplot(1,1,1); cla; reset(gca)
  for ii=1:length(ch_list);
    ch = ch_list(ii);
    h=text(0.5,0.5,ch);
    set(h,'Interp','none','FontName',fontname,'FontUnits','pixels',...
      'FontSize',10,'units','pixels','margin',1);
    set(gca,'xlim',[0 1],'ylim',[0 1]);
    rect = get(h,'extent');
    F = getframe(gca,rect+[0 2 -2 0]);
    %   x=[char(183),'o'];
    bw = ~F.cdata(:,:,1); % Black(1) and white(0) pixel array
    chardata{ii} = bw;
    density(ii) = sum(sum(bw));
    delete(h)
  end
  close(gcf)
  density = density./max(density);
  ch_list_save = ch_list;
  fontname_save = fontname;
  % 03/05/09 sort by density
  [density,order] = sort(density);
  ch_list = ch_list(order);
  chardata = chardata(order);
  %03/05/09 Create a table, based on unique densities
  density_table = struct(...
    'density',0,... % Density, 0..1
    'nchar',0,... % number of characters having this density
    'ch_list',''); % List of characters (of length nchar)
  [u,m,n] = unique(density);
  density_table = repmat(density_table,1,length(u));
  for ii=1:length(m)
    f = find(density == u(ii));
    density_table(ii).density = u(ii);
    density_table(ii).nchar = length(f);
    density_table(ii).ch_list = ch_list(f);
  end
  h_char_pix = size(chardata{1},1); % Height of a character, in pixels
  w_char_pix = size(chardata{1},2); % Width of a character, in pixels

end
density_table_out = density_table;
h_char_pix_out = h_char_pix;
w_char_pix_out = w_char_pix;

function properties = getopt(properties,varargin)
%GETOPT - Process paired optional arguments as 'prop1',val1,'prop2',val2,...
%
%   getopt(properties,varargin) returns a modified properties structure,
%   given an initial properties structure, and a list of paired arguments.
%   Each argumnet pair should be of the form property_name,val where
%   property_name is the name of one of the field in properties, and val is
%   the value to be assigned to that structure field.
%
%   No validation of the values is performed.
%
% EXAMPLE:
%   properties = struct('zoom',1.0,'aspect',1.0,'gamma',1.0,'file',[],'bg',[]);
%   properties = getopt(properties,'aspect',0.76,'file','mydata.dat')
% would return:
%   properties =
%         zoom: 1
%       aspect: 0.7600
%        gamma: 1
%         file: 'mydata.dat'
%           bg: []
%
% Typical usage in a function:
%   properties = getopt(properties,varargin{:})

% Process the properties (optional input arguments)
prop_names = fieldnames(properties);
TargetField = [];
for ii=1:length(varargin)
  arg = varargin{ii};
  if isempty(TargetField)
    if ~ischar(arg)
      error('Propery names must be character strings');
    end
    f = find(strcmp(prop_names, arg));
    if length(f) == 0
      error('%s ',['invalid property ''',arg,'''; must be one of:'],prop_names{:});
    end
    TargetField = arg;
  else
    % properties.(TargetField) = arg; % Ver 6.5 and later only
    properties = setfield(properties, TargetField, arg); % Ver 6.1 friendly
    TargetField = '';
  end
end
if ~isempty(TargetField)
  error('Property names and values must be specified in pairs.');
end


text:
function coordinates = createCoordinates(bounds, step)
    bounds = abs(bounds);
    coordinates = -bounds + step/2:step: bounds - step/2;
    %due to numerical instability of floating point numbers, somtimes series
    %is missing the last point.
    if(abs(coordinates(end) - (bounds - step/2)) > 1e-10)
        coordinates(end + 1) = bounds;
    end
end

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function [x,fval,exitflag,output,lambda,grad,hessian] = fmincon_fix(fix,fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options,varargin)
% Wrapper for fmincon that allows the user to specify a number of model
% parameters that remain fixed to the initial settings.
%
% fix is a binary array.  Zero indicates that the parameter in the
% corresponding position in x0 varies during fitting; one indicates that
% the value remains fixed at the alue in x0.
%
% The other parameters are all as for fmincon.
%
% author: Daniel C Alexander (d.alexander@ucl.ac.uk)
%

% Construct the parameter vector with fixed values removed.
x0f = x0(find(fix==0));

% Run the opt.
[x,fval,exitflag,output,lambda,grad,hessian] = fmincon(fun,x0f,A,b,Aeq,beq,lb,ub,nonlcon,options,varargin{:}, fix, x0);

% Reconstruct the full fitted parameter list including the fixed values.
xf = fix.*x0;
xf(find(fix==0)) = x;
x = xf;

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function K = covPeriodic(hyp, x, z, i)

% Stationary covariance function for a smooth periodic function, with period p
% in 1d (see covPERiso and covPERard for multivariate data):
%
% k(x,z) = sf^2 * exp( -2*sin^2( pi*||x-z||/p )/ell^2 )
%
% where the hyperparameters are:
%
% hyp = [ log(ell)
%         log(p)
%         log(sf) ]
%
% Copyright (c) by Carl Edward Rasmussen and Hannes Nickisch, 2011-01-05.
%
% See also COVFUNCTIONS.M.

if nargin<2, K = '3'; return; end                  % report number of parameters
if nargin<3, z = []; end                                   % make sure, z exists
xeqz = isempty(z); dg = strcmp(z,'diag');                       % determine mode

[n,D] = size(x);
if D>1, error('Covariance is defined for 1d data only.'), end
ell = exp(hyp(1));
p   = exp(hyp(2));
sf2 = exp(2*hyp(3));

% precompute distances
if dg                                                               % vector kxx
  K = zeros(size(x,1),1);
else
  if xeqz                                                 % symmetric matrix Kxx
    K = sqrt(sq_dist(x'));
  else                                                   % cross covariances Kxz
    K = sqrt(sq_dist(x',z'));
  end
end

K = pi*K/p;
if nargin<4                                                        % covariances
    K = sin(K)/ell; K = K.*K; K =   sf2*exp(-2*K);
else                                                               % derivatives
  if i==1
    K = sin(K)/ell; K = K.*K; K = 4*sf2*exp(-2*K).*K;
  elseif i==2
    R = sin(K)/ell; K = 4*sf2/ell*exp(-2*R.*R).*R.*cos(K).*K;
  elseif i==3
    K = sin(K)/ell; K = K.*K; K = 2*sf2*exp(-2*K);
  else
    error('Unknown hyperparameter')
  end
end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function dat = index_loss(algo,dat)

  [x y]=get_xy(dat);
  res = sum((x-y).^2,2);
  lss = length(find(res > 0))/max(size(res));

  dat=data([get_name(dat) ' -> index_loss=' num2str(lss,4) ],[],lss);

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function y = fx(x)
i = find((x > -pi) & (x <= 0));
y = 2*x;
y(i) = 0;
y(1) = y(end);


text:
function pos = ApproxMedianFilter_RB_LED(file)

% This code tests an approximate median filter for separating out the red
% and blue leds mounted on top of the subjects head
% I was using the DataMax system, which was synced using a red LED. Hence,
% I am tracking the x,y position of the animal as well as when the sync
% light comes on.

%file       = 'Mouse12-120808-01.mpg';
%[~,name,~] = fileparts(file);
% Create Text file for writing out LED locations and sync trigger
%fid        = fopen(sprintf('%s.pos',name),'w+');
% Creater readerobj of file
readerobj  = VideoReader(file);
width      = readerobj.Width;
height     = readerobj.Height;
threshF    = 110;    % threshold on foreground pixel intensity
% Initial frame
Fint       = 1;

% Initialize grid for locating centroid of LED
[X,Y] = meshgrid(1:width,1:height);

% Initialize background as a grayscale image of the first frame
bg_bw     = rgb2gray(read(readerobj,Fint));
pos
% Initialize foreground image and difference image
fg          = zeros(size(bg_bw));
fr_diff     = zeros(size(bg_bw));

% Initialize color mask
mask  = zeros(height,width,3,'uint8');

% Initialize fr in case first frame reading returns error
fr  = zeros(height,width,3,'uint8');

% Initialize pos matrix
pos = zeros(readerobj.NumberOfFrames,4);

for i = Fint:readerobj.NumberOfFrames
    %fprintf('%i',i);

    % Access frame of interest - if error (mostly for the last frames,
    % don't know why), will analyze previous frame...
    try
        fr    = read(readerobj,i);
    end
    % convert frame to grayscale
    fr_bw = rgb2gray(fr);
%     keyboard
    %%% Label Color Mask
    label         = repmat(logical(fr_bw>threshF),[1 1 3]);
    mask(label)   = fr(label);
    mask(~label)  = 0;

    %%% Find centroid of remaining pixels %%%
    bw_mask = rgb2gray(mask);
    [L,Nl]  = bwlabel(double(bw_mask));
    ids     = [];

    % Initialize red and blue LEDS as missing (=> -1)
    Rr = [-1 -1];
    Br = [-1 -1];

    % Update centroid if connected components are found
    if Nl > 0
        a = regionprops(L,'PixelList');

        for ii = 1:Nl
            % Skip single pixels
            if size(a(ii).PixelList,1)<2
                continue
            end

            % Access color information and calculate means
            R  = mask(a(ii).PixelList(:,2),a(ii).PixelList(:,1),1);
            mR = mean(R(R>0));
            B  = mask(a(ii).PixelList(:,2),a(ii).PixelList(:,1),3);
            mB = mean(B(B>0));

            if mR > mB && mR > 200
                ids = [ids; a(ii).PixelList];
                Rr   = round(mean(a(ii).PixelList,1));
            elseif mR < mB && mB > 200
                ids = [ids; a(ii).PixelList];
                Br   = round(mean(a(ii).PixelList,1));
            else
                continue
            end

        end
    end

    pos(i,:) = [Rr(1),Rr(2),Br(1),Br(2)];

    % End processing time
    if mod(i,100)==0 & i<1000
        h = waitbar(i/readerobj.NumberOfFrames);
        if 1
            ixr = pos(i-99:i,1);
            ixr = ixr>-1;
            ixb = pos(i-99:i,3);
            ixb = ixb>-1;
            figure(1),clf,
            subplot(3,1,1);imshow(uint8(fr));title('If looks wrong, Ctrl+C and change threshold \!!')
            subplot(3,1,2),imshow(uint8(bw_mask));title('Mask')
            subplot(3,1,3)
                plot(pos(ixr,1),pos(ixr,2),'r')
                hold on
                plot(pos(ixb,3),pos(ixb,4),'b')
            fprintf('Detection of Red LED failed %i times (%i times for the blue LED) in the last 100 frames\n',sum(~ixr),sum(~ixb));
        end
    elseif mod(i,1000)==0
        h = waitbar(i/readerobj.NumberOfFrames);
    end

end
close(h)

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function result = analyzePathDeviationDetection(debug, PerformanceSpec)

% Author: Alexander Wischnewski     Date: 01-07-2019
%
% Description:
%   scans the debug data and ground truth whether the path deviation
%   detection works as specified
%
% Input:
%   debug:          Debug structure of the vehicle
%   PerformanceSpec Structure with fields:
%       HighDeviationLat_m
%       VeryHighDeviationLat_m
%       HighDeviationHeading_rad
%       VeryHighDeviationHeading_rad
%       TimeSpec_s
%
% Outputs:
%   result:       1 if valid according to spec, 0 otherwise

result = 1;

disp(' ');
disp('------------------------------');
disp('Check High Path Deviation Detection: ');
disp('------------------------------');
idxFailure = find(abs(debug.debug_mvdc_path_matching_debug_PathPos_d_m.Data) >= PerformanceSpec.HighDeviationLat_m | ...
    abs(debug.debug_mvdc_path_matching_debug_PathPos_psi_rad.Data) >= PerformanceSpec.HighDeviationHeading_rad, 1);
if isempty(idxFailure)
    disp('The vehicle did not deviate enough from the path for high path deviation detection.');
else
    % load time where failure occured
    tFailure = debug.debug_mloc_statemachine_debug_TUMVehicleState.Time(idxFailure);
    % find time where vehicle has gone to emergency
    idxDetect = find(uint16(debug.debug_mloc_statemachine_debug_TUMVehicleState.Data) == 50, 1);
    tFailureDetected_s = debug.debug_mloc_statemachine_debug_TUMVehicleState.Time(idxDetect);
    disp(['The failure was detected within ' num2str(tFailureDetected_s-tFailure) 's.']);
    if(tFailureDetected_s - tFailure > PerformanceSpec.TimeSpec_s)
        disp('The failure was not detected within the specified time window');
        result = 0;
        return
    end
    disp('The high path deviation failure was detected within the specified time window');
end

disp(' ');
disp('------------------------------');
disp('Check Very High Path Deviation Detection: ');
disp('------------------------------');
idxFailure = find(abs(debug.debug_mvdc_path_matching_debug_PathPos_d_m.Data) >= PerformanceSpec.VeryHighDeviationLat_m | ...
    abs(debug.debug_mvdc_path_matching_debug_PathPos_psi_rad.Data) >= PerformanceSpec.VeryHighDeviationHeading_rad, 1);
if isempty(idxFailure)
    disp('The vehicle did not deviate enough from the path for very high path deviation detection.');
else
    % load time where failure occured
    tFailure = debug.debug_mloc_statemachine_debug_TUMVehicleState.Time(idxFailure);
    % find time where vehicle has gone to emergency
    idxDetect = find(uint16(debug.debug_mloc_statemachine_debug_TUMVehicleState.Data) == 60, 1);
    tFailureDetected_s = debug.debug_mloc_statemachine_debug_TUMVehicleState.Time(idxDetect);
    disp(['The failure was detected within ' num2str(tFailureDetected_s-tFailure) 's.']);
    if((tFailureDetected_s - tFailure) > PerformanceSpec.TimeSpec_s)
        disp('The failure was not detected within the specified time window');
        result = 0;
        return
    end
    disp('The very high path deviation failure was detected within the specified time window');
end

disp(' ');
disp('All checks passed.');

author: TUMFTM


text:
function [airfoildata] = PerformRegression(data)
% Relate the coefficients of the drag polar to Re as parabolic functions.
% Drag polar is Cd = a0 + a1*alpha + a2*alpha^2
% The contents of airfoildata are the coefficients resulting from the
% polynomial regression with respect to Reynolds number:
%'NACA 4 or 5 digit' | 3 coefficients for a0 |
% 3 coefficients for a1 | 3 coefficients for a2 |
% 3 coefficients for alpha_stall
% so that:
% a0 = b1 + b2*log10(Re) + b3*log10(Re)^2
% a1 = c1 + c2*log10(Re) + c3*log10(Re)^2
% a2 = d1 + d2*log10(Re) + d3*log10(Re)^2
% and:
% alphastall = k1 + k2*log10(Re) + k3*log10(Re)^2

for i = 1:length(data)
    if length(data(i).Re) >= 2
        airfoildata{1}{i} = data(i).airfoil;
        airfoildata{2}{i} = [ones(size(data(i).Re)) log10(data(i).Re) log10(data(i).Re).^2]\[data(i).a0];
        airfoildata{3}{i} = [ones(size(data(i).Re)) log10(data(i).Re) log10(data(i).Re).^2]\[data(i).a1];
        airfoildata{4}{i} = [ones(size(data(i).Re)) log10(data(i).Re) log10(data(i).Re).^2]\[data(i).a2];
        airfoildata{5}{i} = [ones(size(data(i).Re)) log10(data(i).Re) log10(data(i).Re).^2]\[data(i).alphastall];
    end
end


text:
%A 2D uniform PDF function.
%xrange and yrange are inclusive limits of region
%z is a 2D sample point.
function val = unifpdf_2d(xrange, yrange, z)
    minX = xrange(1);
    maxX = xrange(2);
    minY = yrange(1);
    maxY = yrange(2);
    evalX = z(1);
    evalY = z(2);
    if(evalX < minX)
        val = 0;
        return;
    elseif(evalX > maxX)
        val = 0;
        return
    elseif (evalY < minY)
        val = 0;
        return;
    elseif(evalY > maxY)
        val = 0;
        return;
    else
        val = 1 / ((maxX - minX) * (maxY - minY));
    end

end


text:
function sc = census(I,J)

% CENSUS  Census Correlation coefficient
%   CENSUS(I,J) computes the census score of matrices I and J:
%
%     CENSUS = 1/N * ~xor( I > Io , J > Jo )
%
%   where N  = prod(size(I))
%         Io = I(central pixel)
%         Jo = J(central pixel)
%   and   size(I) = size(J)
%
%   See also PATCHCORR, ZNCC, SSD

% (c) 2005 Joan Sola

if size(I) ~= size(J)
    error ('Matrices must be the same size.')
else
    s = size(I); % patch size
    if any(iseven(s))
        error('Matrix sizes must be odd.')
    else
        c = (s+1)/2; % patch center
        i = I(c(1),c(2)); % central pixel
        j = J(c(1),c(2)); % central pixel

        V = I > i;
        W = J > j;

        SC = ~xor(V,W);
        sc = sum(sum(SC))/numel(I);
    end
end
% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

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

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

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function Bnew = pollGPS2N(B,x,gpstruct,LB,UB,optimState,options)
%POLLGPS2N Poll 2N fixed basis vectors (Generalized Pattern Search).

nvars = numel(x);

if isempty(B)
    % Axis-oriented poll vector set
    Bnew = [eye(nvars); -eye(nvars)];
else
    Bnew = [];
end

end

author: acerbilab

repo: bads

save_path: github-repos/MATLAB/acerbilab-bads


text:
for i=1:21
    index = i;
    disp(num2str(index));

    fuse_path1 = ['./fused_infrared/fused',num2str(index),'_vgg19_zca1.png'];
    fuse_path2 = ['./fused_infrared/fused',num2str(index),'_vgg19_zca2.png'];
    fuse_path3 = ['./fused_infrared/fused',num2str(index),'_vgg19_zca3.png'];
    fuse_path4 = ['./fused_infrared/fused',num2str(index),'_vgg19_zca4.png'];
%     fuse_path5 = ['./fused_infrared/fused',num2str(index),'_vgg19_zca5.png'];

    image1 = imread(fuse_path1);
    image2 = imread(fuse_path2);
    image3 = imread(fuse_path3);
    image4 = imread(fuse_path4);
%     image5 = imread(fuse_path5);

    imwrite(image1,['./fused_infrared/fused',num2str(index),'_vgg19_nu_zca1.png'],'png');
    imwrite(image2,['./fused_infrared/fused',num2str(index),'_vgg19_nu_zca2.png'],'png');
    imwrite(image3,['./fused_infrared/fused',num2str(index),'_vgg19_nu_zca3.png'],'png');
    imwrite(image4,['./fused_infrared/fused',num2str(index),'_vgg19_nu_zca4.png'],'png');
%     imwrite(image5,['./fused_infrared/fused',num2str(index),'_vgg19_nu_zca5.png'],'png');

    delete(fuse_path1);
    delete(fuse_path2);
    delete(fuse_path3);
    delete(fuse_path4);
%     delete(fuse_path5);
end

author: thfylsty


text:
function [GAmp,GTime]=ADCRadial(p)

global VCtl;
global VObj;

tMiddle=p.tMiddle;

% 2D radial encoding
FOV = VCtl.FOVFreq; % choose FOVFreq as real FOV
Res = VCtl.ResFreq; % choose ResFreq as real resolution
ERes = VCtl.R_SampPerSpoke; % choose R_SampPerSpoke as effective radial resolution

GxAmp=(1/FOV)/((VObj.Gyro/(2*pi))*(1/VCtl.BandWidth));
tHalf=1/(2*(VObj.Gyro/(2*pi))*GxAmp*(FOV/Res));
[GAmp1,GTime1]=StdTrap(tMiddle+VCtl.TEAnchorTime-tHalf-VCtl.MinUpdRate, ...
                       tMiddle+VCtl.TEAnchorTime+tHalf+VCtl.MinUpdRate, ...
                       tMiddle+VCtl.TEAnchorTime-tHalf,               ...
                       tMiddle+VCtl.TEAnchorTime+tHalf,               ...
                       1,2,ERes,2);

GAmp=[GAmp1];
GTime=[GTime1];
[GTime,m,n]=unique(GTime);
GAmp=GAmp(m);

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function  [n, adchans] = plx_ad_chanmap(filename)
% plx_ad_chanmap(filename): return map of raw continuous channel numbers for each channel
%            for the specified .plx or .pl2 file
%
% [n, adchans] = plx_ad_chanmap(filename)
%
% INPUT:
%   filename - if empty string, will use File Open dialog
%
% OUTPUT:
%   n - number of continuous channels
%   adchans - 1 x n array of continuous channel numbers
%
% Normally, there is one channel entry in the .plx for for each raw continuous channel,
% so the mapping is trivial adchans[i] = i-1 (because continuous channels start at 0).
% However, for certain .plx files saved in some ways from OFS (notably after
% loading data files from other vendors), the mapping can be more complex.
% E.g. there may be only 2 non-empty channels in a .plx file, but those channels
% correspond to raw channel numbers 7 and 34. So in this case NChans = 2,
% and adchans[1] = 7, adchans[2] = 34.
% The plx_ routines that return arrays always return arrays of size NChans. However,
% routines that take channels numbers as arguments always expect the raw
% channel number.  So in the above example, to get the data from
% the second channel, use
%   [adfreq, n, ts, fn, ad] = plx_ad(filename, adchans[2])

n = 0;
adchans = -1;

if nargin ~= 1
    error 'Expected 1 input argument';
end

[ filename, isPl2 ] = internalPL2ResolveFilenamePlx( filename );
if isPl2 == 1
    pl2 = PL2GetFileIndex(filename);
    n = numel(pl2.AnalogChannels);
    adchans = 0:(n-1);
    return;
end

[n, adchans] = mexPlex(27,filename);

end


text:
function x = p00_root ( prob, i )

%*****************************************************************************80
%
%% P00_ROOT returns a known root for any problem.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 January 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer PROB, the number of the problem.
%
%    Input, integer I, the index of the root to return.
%
%    Output, real X, the I-th root.
%
  if ( prob == 1 )
    x = p01_root ( i );
  elseif ( prob == 2 )
    x = p02_root ( i );
  elseif ( prob == 3 )
    x = p03_root ( i );
  elseif ( prob == 4 )
    x = p04_root ( i );
  elseif ( prob == 5 )
    x = p05_root ( i );
  elseif ( prob == 6 )
    x = p06_root ( i );
  elseif ( prob == 7 )
    x = p07_root ( i );
  elseif ( prob == 8 )
    x = p08_root ( i );
  elseif ( prob == 9 )
    x = p09_root ( i );
  elseif ( prob == 10 )
    x = p10_root ( i );
  elseif ( prob == 11 )
    x = p11_root ( i );
  elseif ( prob == 12 )
    x = p12_root ( i );
  elseif ( prob == 13 )
    x = p13_root ( i );
  elseif ( prob == 14 )
    x = p14_root ( i );
  elseif ( prob == 15 )
    x = p15_root ( i );
  elseif ( prob == 16 )
    x = p16_root ( i );
  elseif ( prob == 17 )
    x = p17_root ( i );
  elseif ( prob == 18 )
    x = p18_root ( i );
  elseif ( prob == 19 )
    x = p19_root ( i );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P00_ROOT - Fatal error!\n' );
    fprintf ( 1, '  Illegal problem number = %d\n', prob );
    error ( 'P00_ROOT - Fatal error!' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [area] = radius2area(radius)
% r2A returns the area of a circle of radius r.
% Chad Greene 2012
area = pi.*radius.^2;

author: Sable


text:
%% RRT*FN Toolbox Examples
%
%
RRT Examples

%

%
	2D Mobile point-like robot

%
	6 DOF Planar Redundant Manipulator

%
	3 DOF Mobile rectangular robot

%

%

%
RRT* Examples

%

%
	2D Mobile point-like robot

%
	6 DOF Planar Redundant Manipulator

%
	3 DOF Mobile rectangular robot

%

%

%
RRT*FN Examples

%

%
	2D Mobile point-like robot

%
	6 DOF Planar Redundant Manipulator

%
	3 DOF Mobile rectangular robot

%

%

author: olzhas


text:
function [Du,area,Dlambda] = gradu(node,elem,u,Dlambda)
%% GRADU gradient of a finite element function.
%
% Du = GRADU(node,elem,u) compute the gradient of a finite element function
% u on a mesh representing by (node,elem).
%
% [Du,area,Dlambda] = GRADU(node,elem,u) also outputs area and Dlambda
% which is the gradient of P1 basis.
%
% Du = GRADU(node,elem,u,Dlambda) compute the gradient with Dlambda. It
% will save time when Dlambda is available. See recovery.m
%
% See also gradu3, gradbasis
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

if ~exist('Dlambda','var')
    [Dlambda,area] = gradbasis(node,elem);
end
dudx =  u(elem(:,1)).*Dlambda(:,1,1) + u(elem(:,2)).*Dlambda(:,1,2) ...
      + u(elem(:,3)).*Dlambda(:,1,3);
dudy =  u(elem(:,1)).*Dlambda(:,2,1) + u(elem(:,2)).*Dlambda(:,2,2) ...
      + u(elem(:,3)).*Dlambda(:,2,3);
Du = [dudx, dudy];

repo: ifem